Runge-Kutta: Maintaining Units in Numerical Methods

[Roboto]

When conducting numerical methods using 4th Order Runge-Kutta do the physical units have to be maintained?

This never occurred to me until I was writing out all the steps in detail when showing someone I work with the method using a simple projectile motion with drag. It had 4th Order time terms this making the units wrong in the iteration process.

So when using numerical methods, do the units have to be conserved/maintained?

 

[scottdave]

I don't see how any problem can be successfully solved if the units are not preserved. Can you give your example, for clarification with what happened in your situation?

 

[Roboto]

Exactly, that was my thoughts. The first sanity check of something is to check the units.

But does that have to be maintained when using numerical methods to solve something as opposed to algebraic.

I will get my document cleaned up and figure out how to post it here. It several pages.

 

[scottdave]

If it is a missile launch, then I guess you will have derivatives with respect to time. So your steps will be Δt, which should be in the correct unit.

 

[Roboto]

if we take a simple particle velocity decay in the air, mass*acceleration = K*Velocity² where K is drag and geometry coefficents of the particle.
If we let u₁=position, y, and du₁/dt=velocity, and then we let u₂=velocity and du₂/dt=acceleration

With Runge Kutta,
dy/dt=f(u(t),t)
y_n+1=y_n+(h/6)*(k₁+2*k₂+2*k₃+k₄)
t_n+1=t_n+h
k₁=f(u_n,t)
k₂=f(u_n+(h/2)*k₁,t+h/2)
k₃=f(u_n+(h/2)*k₂,t+h/2)
k₄=f(u_n+h*k₃,t+h)

for the first equation f(u_n,t) = u₂

Then substituting u₂ into the above set of steps we get the time iteration for position as follows:
y_n+1=y_n+(u₂*h/6)*(6+3h+h²+h³/4)

The units of u₁ is the same as y_n which are meters. The units of u₂ is meters/sec. The unit of h is time, seconds

Note the term (u₂*h/6) is velocity (m/s) times time (sec) which is distance, meters. So from a units perspective the units are for y_n+1, y_n, and (u₂*h/6) for meters just as one would expect.

But the kicker are the (6+3h+h²+h³/4) terms. The 6 has no units, the 3h has units of time, and h² has units of time², and the h³/4 has units of time³. The units in this last group are not preserved.

As you can see, the units are not preserved when applying 4th order Runge -Kutta to numerically solve the ODE

This was a surprise to me in that I thought that the units must be conserved. But the classical method to numerically solve ODE's clearly shows units are not preserved. Hence the question, when using numerical methods, do units have to be conserved?

 

[scottdave]

How did you get to the yn+1=y_n+(u₂*h/6)*(6+3h+h²+h³/4) ?
I am just not seeing that jump. I have done a Euler's Method with 2 variables, but not Runge Kutta. I did some search and found this, article, which sets it up as a linear algebra style problem.
https://www.myphysicslab.com/explain/runge-kutta-en.html
There is no combining to get squared or cubed h terms in that.

 

[Roboto]

if k₁=u₂
then k₂=u₂+h/2*k₁=u₂+h*u₂/2
Then k₃=u₂+h/2*k₂=u₂+(h/2)*(u₂+h*u₂/2)=u₂+h*u₂/2+h²*u₂/4
then k₄=u₂+h*k₃=u₂+h(u₂+h*u₂/2+h²*u₂/4)=u₂+h*u₂+h²u₂/2+h³u₂/4

thus
y_n+1=y_n+(h/6)*(k₁+2k₂+2k₃+k₄) = y_n+1=y_n+(u₂*h/6)*(6+3h+h²+h³/4)

See how this happens? Crazy right? Units are not preserved, and yet this is one of the most popular methods for numerically solving ODE's

 

[scottdave]

Except that k2 is a function of (u_n+(h/2)*k1,t+h/2), not equivalent to u_n+(h/2)*k1, correct?

 

[Roboto]

In this example, it happens to be that way. The function can have numerous parameters passed to it. In this example, time is not a variable that is being passed.

if we want to get really formal, dy/dt=f(x, dy/dt, dy²/d²t,t). since the original equation was mass*d²y/dt²=K*(dy/dy)²

For Runge Kutta, we convert the 2nd order differential equation into a set of 1st order equations
thus letting u₁=y, du₁/dt = dy/dt Then we let u₂=dy/dt and then du2/dt = d²y/dt²

du₁/dt=u₂ which is the first function = f₁(u₂)
du₂/dt=(K/mass)(u₂)² which is the second function = f₂(u₂)

u_1,n+1=u_1,n+(h/6)*(k11+2k12+2k13+k14)=y_n+1

k11, k12, k13, k14 are for evaluating the first function
k21, k22, k23, k24 are for evaluating the second function

Thus this leads back to the previously posted details.

 

[mpresic3]

It is usually smart to put the problem in dimensionless form before doing the Runge-Kutta integration. For example, if you integrate the equations of motion for a missile launch, scale the distances by the radius of the Earth, scale the velocity by the orbital velocity at the surface, scale the time by one unit of v divided by a unit of r. etc. On the other hand, if you are working with the Hydrogen atom, scale the length by the Bohr radius etc, maybe scale the velocity by the velocity of light. It is usually possible to scale the problem so you lose only a little precision with the calculation. You wouldn't want to scale the position vector of a missile relative to the center of the Earth as multiples of the Bohr radius.

 

[mpresic3]

P.S I wrote many Runge Kutta routines without putting the problem in dimensionless terms from the beginning. Sometimes they worked, and sometimes I ran into trouble. After several years (my earlier advise is not obvious) I always put it in dimensionless form first.

 

[kimbyd]

BTW, if you use Latex, your equations will be clearer (and easier to write).

y_n+1=y_n+(h/6)*(k₁+2*k₂+2*k₃+k₄)

The units of u₁ is the same as y_n which are meters. The units of u₂ is meters/sec. The unit of h is time, seconds

The equation I quoted requires the units of $k_n$ to be the same as the units of $y_n / h$, or meters per second. That should make your units work out.

 

[scottdave]

The k1 through k4 are estimates of y' at point a point then slightly away from that point (see the +h/2 twice then the +h). These different estimates are then weight averaged, to provide a better estimate for y' than the Euler's Method (1st order Runge Kutta). This may explain it for you better than I am - http://lpsa.swarthmore.edu/NumInt/NumIntFourth.html

 

[atyy]

Fourth order Runge-Kutta basically uses the Taylor series expansion. The Taylor series expansion should have the right units.

 

[BvU]

y_n+1=y_n+(h/6)*(k₁+2*k₂+2*k₃+k₄)

is dimensionally correct, so if you substitute correctly, the result should also be correct.

Then substituting u₂ into the above set of steps we get the time iteration for position as follows:
y_n+1=y_n+(u₂*h/6)*(6+3h+h²+h³/4)

conversely, this result can not be right.

As you can see, the units are not preserved when applying 4th order Runge -Kutta to numerically solve the ODE

nonsense

This was a surprise to me in that I thought that the units must be conserved. But the classical method to numerically solve ODE's clearly shows units are not preserved. Hence the question, when using numerical methods, do units have to be conserved?

kudos for your gut feelings! units are conserved. anything else means there is an error.but there is something else being overlooked here: your u has two components, one in meters and one in m/s. likewise your f has two components, one in m/s and one in m/s². idem all k

 

[BvU]

if we want to get really formal, dy/dt=f(x, dy/dt, dy²/d²t,t). since the original equation was mass*d²y/dt²=K*(dy/dy)²

getting formal means you don’t mix notations and work accurately.
if the original equation is really $$m\ddot y = k\dot y^2$$ then there is no $x$ in the problem and $\dot y$ is a not a function of any of those variables you mention directly, but instead the item you want to find through integration (together with $y$ itself)

and the set of first order equations is $${dy\over dt}=\dot y\\\mathstrut \\ {d\dot y\over dt} ={k\over m}\dot y^2$$

For Runge Kutta, we convert the 2nd order differential equation into a set of 1st order equations
thus letting u₁=y, du₁/dt = dy/dt Then we let u₂=dy/dt and then du2/dt = d²y/dt²

du₁/dt=u₂ which is the first function = f₁(u₂)
du₂/dt=(K/mass)(u₂)² which is the second function = f₂(u₂)

u_1,n+1=u_1,n+(h/6)*(k11+2k12+2k13+k14)=y_n+1

k11, k12, k13, k14 are for evaluating the first function
k21, k22, k23, k24 are for evaluating the second function

Thus this leads back to the previously posted details.

looks a lot better than post #1 where $u$ and $k$ only had one index
and of course you mean ‘follow from evaluating’ instead of ‘are for evaluating’

just to reassure us you now know what to do: please re-write post#1 in an absolutely immaculate form...

 

[scottdave]

Just to recap, dimensions (or units if you like) should match. You cannot say 3 miles per hour + 2 miles, for example. This is a powerful tool to see if your answer makes sense. This means if you come up with something like you did (6+3h+h²+h³/4) where h carries dimension (time in this case), then they cannot be added as @BvU stated that you did not think it was correct.
You should suspect that something is not right, but don't jump to the conclusion that the Mathematicians, who created the formulas, are wrong. The more likely scenario is that there is a mistake in your workings, somewhere along the line.

Anybody is capable of a mistake. So yes, it is possible that a textbook or test question has an error on it.

Were you able to figure out where your error was, to get the dimensions to match?