Given the example Differential equation:






y
′

(
t
)
=
t
×


y
(
t
)




{\displaystyle y'(t)=t\times {\sqrt {y(t)}}}


With initial condition:






t

0


=
0


{\displaystyle t_{0}=0}

 and 




y

0


=
y
(

t

0


)
=
y
(
0
)
=
1


{\displaystyle y_{0}=y(t_{0})=y(0)=1}


This equation has an exact solution:





y
(
t
)
=



1
16



(

t

2


+
4

)

2




{\displaystyle y(t)={\tfrac {1}{16}}(t^{2}+4)^{2}}




Task
Demonstrate the commonly used explicit   fourth-order Runge–Kutta method   to solve the above differential equation.

 Solve the given differential equation over the range 



t
=
0
…
10


{\displaystyle t=0\ldots 10}

 with a step value of 



δ
t
=
0.1


{\displaystyle \delta t=0.1}

 (101 total points, the first being given)
 Print the calculated values of 



y


{\displaystyle y}

 at whole numbered 



t


{\displaystyle t}

's (



0.0
,
1.0
,
…
10.0


{\displaystyle 0.0,1.0,\ldots 10.0}

) along with error as compared to the exact solution.


Method summary
Starting with a given 




y

n




{\displaystyle y_{n}}

 and 




t

n




{\displaystyle t_{n}}

 calculate:





δ

y

1


=
δ
t
×

y
′

(

t

n


,

y

n


)



{\displaystyle \delta y_{1}=\delta t\times y'(t_{n},y_{n})\quad }






δ

y

2


=
δ
t
×

y
′

(

t

n


+



1
2



δ
t
,

y

n


+



1
2



δ

y

1


)


{\displaystyle \delta y_{2}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{1})}






δ

y

3


=
δ
t
×

y
′

(

t

n


+



1
2



δ
t
,

y

n


+



1
2



δ

y

2


)


{\displaystyle \delta y_{3}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{2})}






δ

y

4


=
δ
t
×

y
′

(

t

n


+
δ
t
,

y

n


+
δ

y

3


)



{\displaystyle \delta y_{4}=\delta t\times y'(t_{n}+\delta t,y_{n}+\delta y_{3})\quad }


then:






y

n
+
1


=

y

n


+



1
6



(
δ

y

1


+
2
δ

y

2


+
2
δ

y

3


+
δ

y

4


)


{\displaystyle y_{n+1}=y_{n}+{\tfrac {1}{6}}(\delta y_{1}+2\delta y_{2}+2\delta y_{3}+\delta y_{4})}







t

n
+
1


=

t

n


+
δ
t



{\displaystyle t_{n+1}=t_{n}+\delta t\quad }




from math import sqrt
 
def rk4(f, x0, y0, x1, n):
    vx = [0] * (n + 1)
    vy = [0] * (n + 1)
    h = (x1 - x0) / float(n)
    vx[0] = x = x0
    vy[0] = y = y0
    for i in range(1, n + 1):
        k1 = h * f(x, y)
        k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
        k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
        k4 = h * f(x + h, y + k3)
        vx[i] = x = x0 + i * h
        vy[i] = y = y + (k1 + k2 + k2 + k3 + k3 + k4) / 6
    return vx, vy
 
def f(x, y):
    return x * sqrt(y)
 
vx, vy = rk4(f, 0, 1, 10, 100)
for x, y in list(zip(vx, vy))[::10]:
    print("%4.1f %10.5f %+12.4e" % (x, y, y - (4 + x * x)**2 / 16))
 
 0.0    1.00000  +0.0000e+00
 1.0    1.56250  -1.4572e-07
 2.0    4.00000  -9.1948e-07
 3.0   10.56250  -2.9096e-06
 4.0   24.99999  -6.2349e-06
 5.0   52.56249  -1.0820e-05
 6.0   99.99998  -1.6595e-05
 7.0  175.56248  -2.3518e-05
 8.0  288.99997  -3.1565e-05
 9.0  451.56246  -4.0723e-05
10.0  675.99995  -5.0983e-05