Using Runge-Kutta 4 to approximate y(t)

Question

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
--------------------------------

Answer 1

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
--------------------------------