4th order Runge-Kutta for $y' = y$

Question

My question is quite simple, but the more I look at it, the less content I am. My question is how to do a RK4 method for $y'=y$. At first I would assume the following:

$$k_1=y_n$$

$$k_2=y_n+\frac{h}{2}k_1$$

and right here you can actually start to see part of the issue that I have, namely that it seems in the second equation we are adding apples and oranges since $h$ is not unitless. Nonetheless, in continuing onward:

$$k_3=y_n+\frac{h}{2}k_2$$

$$k_4=y_n+h\ k_3$$

and now put all together,

$$y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$$.

So for all of these it seems as though we are adding apples and oranges (my case is exactly this pretty much, finding position by taking the derivative of the velocity). Furthermore, in this last equation to solve for the next step, $y'$ does not appear in it at all, which to me seems really weird (not right at all). To clarify why I am confused on this point, is because normally the first step is an Euler Step, in which case the whole thing put together would look like:

$$y_{n+1}=y_n+h\ y'_{n}$$

Part of the reason I bring this up is because I am given $y'$ and need to now map it's derivative (just like the above Euler Method). I have used RK4 many times on more complicated equations, but for some reason this one just doesn't seem quite right. Thanks!

Answer 1

The answer is quite simple. You are already comparing apples and oranges in the first equation. Garbage in, garbage out. The equation $y'=y$ if written properly is $$dy/dx=y.$$ Do you see it now? To correct it, simply write: $dy/dx=ay,$ where $a$ is a constant and in our example, $a=1$ in units of $1/x$.

Answer 2

You have to keep in mind that you are integrating something. An integral usually looks like this

$$F(x) = \int f(x) \mathop{}\!dx$$

Here $dx$ is also not unitless, so in principle you are multiplying an apple ($f$) by something ($dx$) to obtain an orange ($F$).

The issue that $y'$ does not show up is unclear. Perhaps it becomes clear when I restate the Runge-Kutta scheme here

Given $y' = f(t,y)$ with initial value $y(t_0) = y_0$ you can calculate the discrete time series $y_n = y(t_n)$ with $$ y_{n+1} = y_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k4) \;,\quad t_{n+1} = t_n + h$$ where $$\begin{aligned} k_1 &= f(t,y_n) \\ k_2 &= f(t+h/2 ,y_n + h k_1/2) \\ k_3 &= f(t+h/2 ,y_n + h k_2/2) \\ k_4 &= f(t+h ,y_n + h k_3) \\ \end{aligned}$$

---

Here a quick RK4 implementation in Julia to show that there is nothing wrong.

using Plots

# Function definition y' = y
f(t,y) = y

# Boundary condition y(0) = 1
y = 1

# Time step
dt = 0.01

# Arrays for time an solution
time = collect(0:dt:1)
sol = similar(time)

# integration loop
for (n,t) in enumerate(time)
    sol[n] = y
    k1 = f(t     ,y        )
    k2 = f(t+dt/2,y+dt/2*k1)
    k3 = f(t+dt/2,y+dt/2*k2)
    k4 = f(t+dt  ,y+dt  *k3)
    y  = y + dt/6*(k1 + 2*k2 + 2*k3 + k4)
end

plot(time,sol)
plot!(time,exp(time))

Answer 3

This is a tricky question because mathematicians usually ignore units. I actually ran into this myself when making DifferentialEquations.jl compatible with units. In this notebook, I use Unitful.jl's unitful-arithmetic. Essentially, you just use Julia's string macros to make u"s" put units in seconds and implicit multiplication to make 2u"s" mean "two seconds. Knowing this, notice that

using DifferentialEquations
f = (t,y) -> 0.5*y
u0 = 1.5u"N"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5())

that throws an error, and is very close to your equation. Why? The error is very clear:

DimensionError: 1.5 N and 0.015 N s are not dimensionally compatible.

This is exactly what you're seeing. It calculates $hk_1$ as having units of $Ns$ to generate the intermediate estimated value, but this value needs to have units of $N$. Is it wrong? No, I was wrong! Remember that $f(t,u) = y\prime = \frac{dy}{dt}$. It should be a rate, so it needs to be $y$ per second. When you do that:

f = (t,y) -> 0.5*y/3.0u"s"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5()) 

# the printout

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Array{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},1}:
      0.0 s
 0.143116 s
      1.0 s
u: 3-element Array{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},1}:
     1.5 N
 1.53621 N
 1.77204 N

you see that it works. So yes, $k$ is a rate, so it needs to be in units per time, and then when you multiply by $h$, you now get a value in units again.

What's actually going on? $hk_1$ is just an estimate of $y(t+h)$ (notice the units line up). That step is just Euler. You then use $y(t)$ and the estimate of $y(t+h)$ to make an estimate of $y(t+h/2)$ (some kind of midpoint). You then use the current $y(t+h/2)$ to correct the estimate at $y(t+h/2)$ (a kind of Picard iteration). You then use the corrected $y(t+h/2)$ estimate to take a half step and get a $y(t+h)$ estimate. Then you average this $y(t+h)$ estimate with the first one (remember $hk_1$?), and double the changes to the midpoint to get two more estimates for the final step, and take a weighted average of all of these together. That is RK4. And if you check the units on this, they all line up.

So $y\prime = y$ is bad shorthand if you care about units.