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 }




INCLUDE "D2:PRINTF.ACT" ;from the Action! Tool Kit
INCLUDE "H6:REALMATH.ACT"
 
DEFINE PTR="CARD"
 
REAL one,two,four,six
 
PROC Init()
  IntToReal(1,one)
  IntToReal(2,two)
  IntToReal(4,four)
  IntToReal(6,six)
RETURN
 
PROC Fun=*(REAL POINTER x,y,res)
  DEFINE JSR="$20"
  DEFINE RTS="$60"
  [JSR $00 $00 ;JSR to address set by SetFun
   RTS]
 
PROC SetFun(PTR p)
  PTR addr
 
  addr=Fun+1 ;location of address of JSR
  PokeC(addr,p)
RETURN
 
PROC Rate(REAL POINTER x,y,res)
  REAL tmp
 
  Sqrt(y,tmp)         ;tmp=sqrt(y)
  RealMult(x,tmp,res) ;res=x*sqrt(y)
RETURN
 
PROC RK4(PTR f REAL POINTER dx,x,y,res)
  REAL k1,k2,k3,k4,dx2,k12,k22,tmp1,tmp2,tmp3
 
  SetFun(f)
  Fun(x,y,tmp1)        ;tmp1=f(x,y)
  RealMult(dx,tmp1,k1) ;k1=dx*f(x,y)
 
  RealDiv(dx,two,dx2)  ;dx2=dx/2
  RealDiv(k1,two,k12)  ;k12=k1/2
  RealAdd(x,dx2,tmp1)  ;tmp1=x+dx/2
  RealAdd(y,k12,tmp2)  ;tmp2=y+k1/2
  Fun(tmp1,tmp2,tmp3)  ;tmp3=f(x+dx/2,y+k1/2)
  RealMult(dx,tmp3,k2) ;k2=dx*f(x+dx/2,y+k1/2)
 
  RealDiv(k2,two,k22)  ;k22=k2/2
  RealAdd(y,k22,tmp2)  ;tmp2=y+k2/2
  Fun(tmp1,tmp2,tmp3)  ;tmp3=f(x+dx/2,y+k2/2)
  RealMult(dx,tmp3,k3) ;k3=dx*f(x+dx/2,y+k2/2)
 
  RealAdd(x,dx,tmp1)   ;tmp1=x+dx
  RealAdd(y,k3,tmp2)   ;tmp2=y+k3
  Fun(tmp1,tmp2,tmp3)  ;tmp3=f(x+dx,y+k3)
  RealMult(dx,tmp3,k4) ;k4=dx*f(x+dx,y+k3)
 
  RealAdd(k2,k3,tmp1)     ;tmp1=k2+k3
  RealMult(two,tmp1,tmp2) ;tmp2=2*k2+2*k3
  RealAdd(k1,tmp2,tmp1)   ;tmp3=k1+2*k2+2*k3
  RealAdd(tmp1,k4,tmp2)   ;tmp2=k1+2*k2+2*k3+k4
  RealDiv(tmp2,six,tmp1)  ;tmp1=(k1+2*k2+2*k3+k4)/6
  RealAdd(y,tmp1,res)     ;res=y+(k1+2*k2+2*k3+k4)/6
RETURN
 
PROC Calc(REAL POINTER x,res)
  REAL tmp1,tmp2
 
  RealMult(x,x,tmp1)      ;tmp1=x*x
  RealDiv(tmp1,four,tmp2) ;tmp2=x*x/4
  RealAdd(tmp2,one,tmp1)  ;tmp1=x*x/4+1
  Power(tmp1,two,res)     ;res=(x*x/4+1)^2
RETURN
 
PROC RelError(REAL POINTER a,b,res)
  REAL tmp
 
  RealDiv(a,b,tmp)     ;tmp=a/b
  RealSub(tmp,one,res) ;res=a/b-1
RETURN
 
PROC Main()
  REAL x0,x1,x,dx,y,y2,err,tmp1,tmp2
  CHAR ARRAY s(20)
  INT i,n
 
  Put(125) PutE() ;clear the screen
  MathInit()
  Init()
  PrintF("%-2S %-11S %-8S%E","x","y","rel err")
 
  IntToReal(0,x0)
  IntToReal(10,x1)
  ValR("0.1",dx)
 
  RealSub(x1,x0,tmp1)    ;tmp1=x1-x0
  RealDiv(tmp1,dx,tmp2)  ;tmp2=(x1-x0)/dx
  n=RealToInt(tmp2)      ;n=(x1-x0)/dx
  i=0
  IntToReal(1,y)
  DO
    IntToReal(i,tmp1)      ;tmp1=i
    RealMult(dx,tmp1,tmp2) ;tmp2=i*dx
    RealAdd(x0,tmp2,x)     ;x=x0+i*dx
 
    IF i MOD 10=0 THEN
      Calc(x,y2)
      RelError(y,y2,err)
      StrR(x,s) PrintF("%-2S ",s)
      StrR(y,s) PrintF("%-11S ",s)
      StrR(err,s) PrintF("%-8S%E",s)
    FI
 
    i==+1
    IF i>n THEN EXIT FI
 
    RK4(rate,dx,x,y,tmp1)  ;tmp1=rk4(rate,dx,x0+dx*(i-1),y)
    RealAssign(tmp1,y)     ;y=rk4(rate,dx,x0+dx*(i-1),y)
  OD
RETURN