# Numerical method for first-order non-linear differential equation

So I'm modeling a cycler on a 1000 m race track using the equation $$\frac{dv}{dt}=\frac{P}{m}\left(\frac{1}{v}\right)-\frac{k}{m}\left(v^2\right)-ug$$ where $$P$$ = power, $$m$$ = mass, $$v$$ = velocity of cycler, $$k$$ = drag coefficient, $$u$$ = kinetic friction and $$g$$ = gravity.

Now logical, the cycler will start at $$v = 0$$ for $$t = 0$$, i.e. $$v(0)=0$$ (initial condition) but when I use Euler's Method, it says that the initial condition is not part of the domain.

Now my questions

1. Then what should I use as the initial condition
2. How can I calculate the next $$v$$, i.e $$v(t+h$$) where $$h$$ is the step size
3. How can I simulate the race and plot the velocity-time graph?
4. How can I plot the distance-time graph?

What would be a more accurate model? Thank you

• Welcome to this site. You should use MathJax for showing formulas. Oct 5, 2019 at 11:21
• Where did you get this model? Oct 5, 2019 at 11:46
• We're advised 2 develop a more realistic model. This is just 2 stimulate sum race with different parameters. What would a more realistic model then? Oct 6, 2019 at 20:40
• @Schematic There are several models in the literature for simulating cyclists on a track. Including one of mine. Others have also developed models. Look for papers by Olds or by Martin (they are cited in both of the linked articles) for the original cycling models. Oct 8, 2019 at 0:31

Assuming that the model is indeed correct, there is obviously a problem with using $$v(t=0)=0$$ due to the division by the velocity. There are a few ways I can think of to fix this:

1. Use a non-zero (but small) starting velocity (e.g., $$v(t=0)=10^{-10}$$).
2. Ignore the $$1/v$$ term for the first step
3. Add a very small value (e.g., machine-precision) to the denominator so that it's $$1/(v+\epsilon)$$ instead of $$1/v$$
4. Determine a velocity at $$t=-\mathrm{d}t/2$$ and integrate forwards from there (such that the first period "skips over" the $$t=0$$ issue)
5. Use implicit methods so that it's $$1/v(t=\mathrm{d}t)$$ instead of $$1/v(t=0)$$ you're doing now.
• Alternatively, you could even just use a root finder to find the $$v(t=\mathrm{d}t)$$ given $$v=0$$ & $$\mathrm{d}t$$ value & then integrate using explicit method

The last two are probably the most difficult but possibly the most accurate of them. The first is probably the easiest to implement. The second can be easily done as well, but if you think the power term is significant when accelerating, it's probably not good to skip out on it for the first step.1
The third can be dangerous because when $$v=0$$, then $$1/\epsilon$$ is actually very big and can lead to erroneous calculations for all $$t>0$$. Using large values of $$\epsilon$$ can also mean the $$\epsilon$$ is still (falsely) contributing to the future values.

Once you choose one of these options, the remaining questions can easily be discovered.

As an example as why one shouldn't use $$v+\epsilon$$ with too large a $$\epsilon$$, consider the following diagram.2 The purple line uses the bullet under item 5 to determine the first $$\mathrm{d}t$$ (i.e., finds the $$v(t=\mathrm{d}{t})$$ given the $$\mathrm{d}t$$) while the the green line uses $$\epsilon=10^{-3}$$. Other parameters are: $$\mathrm{d}t=10^{-5}$$, $$P=10^3$$, $$m=10$$, $$u=0.2$$ and $$k=100$$. As you can clearly see, adding the artificial value of $$\epsilon$$ to avoid a numerical issue introduces a jump in the data and then converges to the wrong solution.
Hence, using such a means to avoid divide-by-zero cases should be discouraged.

1. This can possibly be mitigated to some degree by using a very small first time step, e.g., $$\mathrm{d}t\sim10^{-8}$$, before continuing with a more practical $$10^{-3}$$ or whatever it needs to be.
2. Code that generated this is posted to my Github page. I wrote it in Rust, but should be easily transcribed into other, familiar languages.

• I suspect that the model is based on $P = F v$ and either (a) assuming constant power rather than constant force (which is probably a bad assumption) or (b) assumption some functional form for $P(v)$ (which the OP has missed out on in the problem statement.) Oct 6, 2019 at 22:10
• @Seifert Assumption made was that the power is constant. We are advised to develop a more accurate model. Oct 7, 2019 at 20:44
• @Michael I was assuming the first of those as well. Oct 7, 2019 at 23:06

So when velocity is small, you can estimate acceleration with $$\dot{v} = \frac{P}{m v}$$ (ignores small terms) which you can directly solve with \begin{aligned} x(v) &= \int \frac{v}{\dot{v}}\,{\rm d}v = \frac{v^3}{3 \tfrac{P}{m} } \\ t(v) &= \int \frac{1}{\dot{v}}\,{\rm d}v = \frac{v^2}{2 \tfrac{P}{m} } \end{aligned} \tag{1}

So initially you have, \begin{aligned} v(t) &= \sqrt{ \frac{2 P t}{m} } \\ x(t) &= \sqrt{ \frac{8 P t^3}{9 m} } \end{aligned} \tag{2}

You can use the above for the first step when $$t>0$$.

But, generally, I do not recommend using the Euler method here. Even, the midpoint method would be a huge improvement in stability and accuracy.

BTW. The overall equation $$\dot{v} = \frac{a}{v} - \beta v^2 - g v$$ does have an analytical solution, which I cannot write here because it is too long, but any CAS system can produce rather easily.

Use the equations in (1) to solve directly if you wish.

As a side note, I do question the contribution of gravity as a velocity related term. The way I see it, the equation of motion should be

$$\dot{v} = \frac{P}{m v} - \beta v^2 - g \sin\theta$$ where $$\theta$$ is the incline angle, with positive being uphil, and negative downhill.