Need help solving a differential equation numerically

[MrOmar]

Hello all, I am currently trying to solve a differential equation numerically. The equation is as follows:

dv/dt = (u*q)/(m₀-q*t) - g - (c_d*ρ*A*0.5)*v²

If you haven't already guessed, it's the rocket equation with added gravity and drag. Now, I'm not even sure if that's what it's supposed to look like. Anyway, u is the speed of the exhaustion gasses of the rocket, q is the amount of gas spewed out every second in kg/s, m₀ is the total mass of the rocket, t is time, g is gravity, c_d is the drag coefficient, ρ is the density of air, A is the reference area of the rocket that is affected by the air, and finally v² is, well, velocity squared.

Glad that's out of the way.

For a way to make this simpler for myself, I've assumed that the air density is constant. So now my equation would look like this:

dv/dt = (u*q)/(m₀-q*t) - g - k*v²

Where k is a constant. I have all the data for the various constants, and now I can solve this numerically, right? Wrong.

To solve this numerically I need a difference quotient, and I have just that (with some help from my teacher):

(v(n+1)-v(n))/(t(n+1)-t(n))=k1/t(n)-k2-k3*v(n)²

where t(n+1)-t(n)=h is the stepsize.

Isolating v(n+1) I get v(n+1)=(k1/t(n)-k2-k3*v(n)²)*h+v(n)

Now, m₀-q*t is the mass of the rocket as a function of time. My teacher calls it t(n) in the difference quotient, while I also have a t₀, the time. I don't get it. But I tried solving it anyway, using Eulers method, and maybe eventually probably using RK4 at some point. I have never done this before, mind you, and I have just learned about this, by myself. My teacher is not answering the emails I've sent him, and he won't (read: He's not allowed to) help me anyway, as I am writing an assignment where I have to delve into the math and physics by myself without much help.

So I try to solve this using excel, and set some start values: y(10)=97 where t₀ is 10 and v₀ is 97. v(n) in my difference quotient, in the very first step, is 97, and I have to use that to find the next one hence the reason I isolated v(n+1).

But when I insert all the values into the difference quotient, and try to calculate my step, my velocity goes from 97 m/s to -3600. I tried using a stepsize of 1, and 0.5 but the end result was always the same; it would jump drastically down below 0

I really, really, don't know what I'm doing wrong, I don't know if my difference quotient is correct, and at this point, I'm not even sure if my rocket equation is correct. I've messaged my teacher but he's not answering.
Any help would be appreciated. Sorry for the long post, hope it was coherent enough to understand.

 

[rikblok]

It looks like you're on the right track but I think you're using t(n) to mean two different things?

t(n+1)-t(n)=h

v(n+1)=(k1/t(n)-k2-k3*v(n)2)*h+v(n)

It might be better to define the mass of the rocket as $m(n)=m_0-q t(n)$.

Then things look good. You've got v(n+1) defined strictly in terms of the state of the system at step n, which is just what you need for Euler's method.

 

[MrOmar]

It looks like you're on the right track but I think you're using t(n) to mean two different things?It might be better to define the mass of the rocket as $m(n)=m_0-q t(n)$.

Then things look good. You've got v(n+1) defined strictly in terms of the state of the system at step n, which is just what you need for Euler's method.

Thanks for your input, but I am afraid it is still not working. I think I may be doing something wrong, but I'm not sure. I'll keep trying though, see if I maybe stumble upon a solution by accident..

 

[Ray Vickson]

Hello all, I am currently trying to solve a differential equation numerically. The equation is as follows:

dv/dt = (u*q)/(m₀-q*t) - g - (c_d*ρ*A*0.5)*v²

If you haven't already guessed, it's the rocket equation with added gravity and drag. Now, I'm not even sure if that's what it's supposed to look like. Anyway, u is the speed of the exhaustion gasses of the rocket, q is the amount of gas spewed out every second in kg/s, m₀ is the total mass of the rocket, t is time, g is gravity, c_d is the drag coefficient, ρ is the density of air, A is the reference area of the rocket that is affected by the air, and finally v² is, well, velocity squared.

Glad that's out of the way.

For a way to make this simpler for myself, I've assumed that the air density is constant. So now my equation would look like this:

dv/dt = (u*q)/(m₀-q*t) - g - k*v²

Where k is a constant. I have all the data for the various constants, and now I can solve this numerically, right? Wrong.

To solve this numerically I need a difference quotient, and I have just that (with some help from my teacher):

(v(n+1)-v(n))/(t(n+1)-t(n))=k1/t(n)-k2-k3*v(n)²

where t(n+1)-t(n)=h is the stepsize.

Isolating v(n+1) I get v(n+1)=(k1/t(n)-k2-k3*v(n)²)*h+v(n)

Now, m₀-q*t is the mass of the rocket as a function of time. My teacher calls it t(n) in the difference quotient, while I also have a t₀, the time. I don't get it. But I tried solving it anyway, using Eulers method, and maybe eventually probably using RK4 at some point. I have never done this before, mind you, and I have just learned about this, by myself. My teacher is not answering the emails I've sent him, and he won't (read: He's not allowed to) help me anyway, as I am writing an assignment where I have to delve into the math and physics by myself without much help.

So I try to solve this using excel, and set some start values: y(10)=97 where t₀ is 10 and v₀ is 97. v(n) in my difference quotient, in the very first step, is 97, and I have to use that to find the next one hence the reason I isolated v(n+1).

But when I insert all the values into the difference quotient, and try to calculate my step, my velocity goes from 97 m/s to -3600. I tried using a stepsize of 1, and 0.5 but the end result was always the same; it would jump drastically down below 0

I really, really, don't know what I'm doing wrong, I don't know if my difference quotient is correct, and at this point, I'm not even sure if my rocket equation is correct. I've messaged my teacher but he's not answering.
Any help would be appreciated. Sorry for the long post, hope it was coherent enough to understand.

There are two issues:
(1) Does the differential equation make sense?
(2) How to solve the DE numerically and reliably.

As to (1): it looks to me like the answer is no. Assuming constant air density, you can write the DE as
$$v'(t) = \frac{a}{b-t} - g - c v(t)^2$$
where all the parameters are $> 0$. As $t \uparrow b$ we have $v(t) \to +\infty$, so there the acceleration goes to infinity! That does not look reasonable, especially in the presence of air resistance. So, IMHO there is something very wrong with your DE.

(2) Assuming your DE is OK, the singularity at $t = b$ will cause any numerical method to fail: no numerical method will be able to produce a reliable solution getting close to $t = b$.

In general, the (forward) Euler method you presented is just about the worst possible way to try to solve a DE numerically. The reason that so many fancy alternative methods have been developed is from pure necessity: getting reliable solutions for a variety of sometimes stiff, sometimes unstable DEs is tricky. Some methods, such as Runge-Kutta, are pretty good overall, but some other DEs need even more sophisticated methods. You would be well advised to go to the library and check out a book on numerical solutions of DEs, or failing that, look up several web pages that deal with the issues.

 

[MrOmar]

There are two issues:
(1) Does the differential equation make sense?
(2) How to solve the DE numerically and reliably.

As to (1): it looks to me like the answer is no. Assuming constant air density, you can write the DE as
$$v'(t) = \frac{a}{b-t} - g - c v(t)^2$$
where all the parameters are $> 0$. As $t \uparrow b$ we have $v(t) \to +\infty$, so there the acceleration goes to infinity! That does not look reasonable, especially in the presence of air resistance. So, IMHO there is something very wrong with your DE.

(2) Assuming your DE is OK, the singularity at $t = b$ will cause any numerical method to fail: no numerical method will be able to produce a reliable solution getting close to $t = b$.

In general, the (forward) Euler method you presented is just about the worst possible way to try to solve a DE numerically. The reason that so many fancy alternative methods have been developed is from pure necessity: getting reliable solutions for a variety of sometimes stiff, sometimes unstable DEs is tricky. Some methods, such as Runge-Kutta, are pretty good overall, but some other DEs need even more sophisticated methods. You would be well advised to go to the library and check out a book on numerical solutions of DEs, or failing that, look up several web pages that deal with the issues.

Thank you! Atleast I now know what's wrong with the whole thing. I have some questions though. I am not familiar with the up arrow in this quote -> $t \uparrow b$, would you mind explaining what it means in that context?
And it looks like I need to rewrite my equation, but I am way too tired to see what's wrong with it right now. Are there any big problems where I would have to start from scratch, or is it a matter of not having the $v(t) \to +\infty$?

 

[Ray Vickson]

Thank you! Atleast I now know what's wrong with the whole thing. I have some questions though. I am not familiar with the up arrow in this quote -> $t \uparrow b$, would you mind explaining what it means in that context?
And it looks like I need to rewrite my equation, but I am way too tired to see what's wrong with it right now. Are there any big problems where I would have to start from scratch, or is it a matter of not having the $v(t) \to +\infty$?

The notation $t \uparrow b$ means $t \to b$ from below; that is, it is a left-hand-limit at $t = b$. Similarly, the notation $t \downarrow b$ would mean $t \to b$ from above, so would be a right-hand-limit.

I cannot deal with your second question; you need to consult the literature on "rocket equations" in the presence of drag, if you can find it. I am sure that the problem has been studied---for obvious reasons---but it may not be available freely on-line. You might need to go to the library. There may also be quite a lot of such work that is still "Classified", but I would be willing to bet that at least some of it is in the public domain.

 

[MrOmar]

I cannot deal with your second question; you need to consult the literature on "rocket equations" in the presence of drag, if you can find it. I am sure that the problem has been studied---for obvious reasons---but it may not be available freely on-line. You might need to go to the library. There may also be quite a lot of such work that is still "Classified", but I would be willing to bet that at least some of it is in the public domain.

No, I can't find anything on the subject that I can use, I have been looking for days now. Anyway, thanks for your help.