0
$\begingroup$

I want to use Runge-Kutta 4 to solve the IVP $$ y' = y, \ y(0) = 33 $$ numerically and I want to find $y(2^k)$ where $k = 1,2,3,4...$ with an accuracy of two correct decimals.

I've created a python program for this where I'm using the RK4 method iterating from
$t=0$ to $t=2^k$

However I can't seem to get this to work for $k \gt 4$, for $k=5$ I get an error of 12.5 if I use $dt = 10^{-4}$.
If I change the step length to $dt = 10^{-5}$ I end up with an error of 42.5 for $k = 5$
For $dt = 10^{-6}$ it gets even larger.

I would like to be able to at least run this for $k = 10$ with an accuracy of two correct decimals if possible.

I would greatly appreciate any help I can get, thanks in advance!
I'll paste my code below:

import numpy as np

# D.E.: y' = y
def f(y):
    return y

# Runge-Kutta 4
def rk4(y0,t,dt):
    for _ in range(len(t)):
        k1 = (f(y0))
        k2 = (f(y0 + dt*k1/2))
        k3 = (f(y0 + dt*k2/2))
        k4 = (f(y0 + dt*k3)) 
        k = dt*(k1 + 2*(k2 + k3) + k4)/6
        y0 = y0 + k
    return y0

# Exact function
def yfunc(t):
    return 33*np.exp(t)

# Initial values
t0 = 0
y0 = 33
dt = 10**-4
n = 5

# Loop
for i in range(n):

    tn = 2**(i+1) # tn = 2^k where k = 1,2,3,4...
    t = np.arange(t0,tn,dt)
    y = rk4(y0,t,dt)
    error = abs(y-yfunc(tn))

    print(f"--------------2^{i+1}---------------")
    print(f"Exact value: {yfunc(tn)}")
    print(f"Found value: {y}")
    print(f"Error: {error}")
    print(f"--------------------------------")

When I run this code I get the following result:

--------------2^1---------------
Exact value: 243.83885126471145
Found value: 243.83885126471452
Error: 3.069544618483633e-12
--------------------------------
--------------2^2---------------
Exact value: 1801.7389510937599
Found value: 1801.7389510937614
Error: 1.5916157281026244e-12
--------------------------------
--------------2^3---------------
Exact value: 98371.61357237703
Found value: 98371.61357237639
Error: 6.402842700481415e-10
--------------------------------
--------------2^4---------------
Exact value: 293241647.1767598
Found value: 293241647.1767584
Error: 1.3709068298339844e-06
--------------------------------
--------------2^5---------------
Exact value: 2605777686028462.5
Found value: 2605777686028450.0
Error: 12.5
--------------------------------
$\endgroup$

1 Answer 1

3
$\begingroup$

You got 14 correct digits in the last case, that is slightly better than to be expected. For better results you need to change from 64bit floating point to a multi-precision format, for instance mpmath. You get a little more if you apply a compensated or Kahan summation to the step accumulation, which is an ad-hoc solution for a little bit of multi-precision.

If you add $y+h·k$ (in your code, $k$ already contains the factor $h$), with $y\sim 1$, $k\sim y$ and $h=10^{-4}$, then the shift of the second operand will lose 4 trailing digits from the update. These digits are irrecoverably lost, the information supplied to the sum is only good for 12 digits, so the error in the result is expected to have a similar accuracy.

To avoid this, capture the lost digits with some extra operations

# Runge-Kutta 4
def rk4(y0,t,dt):
    dy = 0
    for _ in range(len(t)):
        k1 = (f(y0))
        k2 = (f(y0 + dt*k1/2))
        k3 = (f(y0 + dt*k2/2))
        k4 = (f(y0 + dt*k3))
        dy += dt*(k1 + 2*(k2 + k3) + k4)/6  # add the increment to the previous overflow
        y1, y0 = y0, y0 + dy                # perform the lossy addition
        dy -= y0-y1                         # remove the actual increment from the update to get the new overflow
    return y0
   

With that I get almost 16 digits accuracy

--------------2^1---------------
Exact value: 243.83885126471145
Found value: 243.83885126471148
Error: 2.842170943040401e-14
rel.Err.: 1.1655939684340703e-16
--------------------------------
--------------2^2---------------
Exact value: 1801.7389510937599
Found value: 1801.7389510937603
Error: 4.547473508864641e-13
rel.Err.: 2.523935837710596e-16
--------------------------------
--------------2^3---------------
Exact value: 98371.61357237703
Found value: 98371.61357237707
Error: 4.3655745685100555e-11
rel.Err.: 4.437839748656839e-16
--------------------------------
--------------2^4---------------
Exact value: 293241647.1767598
Found value: 293241647.17676
Error: 2.384185791015625e-07
rel.Err.: 8.130447410761157e-16
--------------------------------
--------------2^5---------------
Exact value: 2605777686028462.5
Found value: 2605777686028467.0
Error: 4.5
rel.Err.: 1.7269316657855668e-15
--------------------------------
$\endgroup$
2
  • $\begingroup$ Thanks a lot for the help, I've never seen that you can capture the lost digits with that kind of extra operations - very helpful. I did add this to my code and I also imported "Decimal" in order to get bigger numbers since at 2^10 the result gets too big to handle otherwise - it looks much better now. Once again - thanks! $\endgroup$
    – Minea
    Oct 11, 2022 at 16:13
  • 1
    $\begingroup$ It got a bit separated in writing up the two ideas, but the summation trick is, as mentioned, Kahan summation (usually introduced to increase the accuracy of series summation). I expounded on the background for this in ODE solvers in scicomp.stackexchange.com/q/37425/…. $\endgroup$ Oct 11, 2022 at 16:34

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .