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 }




with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;
procedure RungeKutta is
   type Floaty is digits 15;
   type Floaty_Array is array (Natural range <>) of Floaty;
   package FIO is new Ada.Text_IO.Float_IO(Floaty); use FIO;
   type Derivative is access function(t, y : Floaty) return Floaty;
   package Math is new Ada.Numerics.Generic_Elementary_Functions (Floaty);
   function calc_err (t, calc : Floaty) return Floaty;
 
   procedure Runge (yp_func : Derivative; t, y : in out Floaty_Array;
                    dt : Floaty) is
      dy1, dy2, dy3, dy4 : Floaty;
   begin
      for n in t'First .. t'Last-1 loop
         dy1 := dt * yp_func(t(n), y(n));
         dy2 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy1 / 2.0);
         dy3 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy2 / 2.0);
         dy4 := dt * yp_func(t(n) + dt, y(n) + dy3);
         t(n+1) := t(n) + dt;
         y(n+1) := y(n) + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
      end loop;
   end Runge;
 
   procedure Print (t, y : Floaty_Array; modnum : Positive) is begin
      for i in t'Range loop
         if i mod modnum = 0 then
            Put("y(");   Put (t(i), Exp=>0, Fore=>0, Aft=>1);
            Put(") = "); Put (y(i), Exp=>0, Fore=>0, Aft=>8);
            Put(" Error:"); Put (calc_err(t(i),y(i)), Aft=>5);
            New_Line;
         end if;
      end loop;
   end Print;
 
   function yprime (t, y : Floaty) return Floaty is begin
      return t * Math.Sqrt (y);
   end yprime;
   function calc_err (t, calc : Floaty) return Floaty is
      actual : constant Floaty := (t**2 + 4.0)**2 / 16.0;
   begin return abs(actual-calc);
   end calc_err;   
 
   dt : constant Floaty := 0.10;
   N : constant Positive := 100;
   t_arr, y_arr : Floaty_Array(0 .. N);
begin
   t_arr(0) := 0.0;
   y_arr(0) := 1.0;
   Runge (yprime'Access, t_arr, y_arr, dt);
   Print (t_arr, y_arr, 10);
end RungeKutta;