Runge-Kutta of order 4 to solve m first order IVPs using Mathematica

Question

I have the following m first-order IVPs

d(y1)/dx=f1(x,y1,…,ym), y1(x0)=y10,
d(y2)/dx=f2(x,y1,…,ym), y2(x0)=y20,

⋮

d(ym)/dx =fm(x,y1,…,ym), ym(x0)=ym0.

The above system can be written using vector notations as

dY/dx =F(x,Y ), Y(x0) = Y0,

where Y =(y1,…,ym), F =(f1,…,fm), and Y0=(y10,…,ym0). By applying the RK4 method

k1=h*F(xn,Y(xn)), 
k2=h*F(xn + 0.5h,Y(xn) + 0.5k1),
k3=h*F(xn + 0.5h, Y(xn) + 0.5k2),
k4=h*F(xn + h, Y(xn) + k3)

the solutions are then given byY(xn+1) = Y(xn)+1/6(k1+2k2+2k3+k4).

Example :

x1’=x2, x1(0)=1,
x2=x3, x2(0)=2,
x3’=x4, x3(0)=3,
x4’= −8*x1+sin(t)*x2−3*x3 + t^2, x4(0)=4.

How to write a code to implement the Runge-Kutta 4th order to solve the above m first order IVPs using Mathematica.

https://mathematica.stackexchange.com/questions/23516/solving-a-system-of-odes-with-the-runge-kutta-method?rq=1, https://mathematica.stackexchange.com/questions/194409/numerical-solution-using-rk4-to-solve-nonlinear-ode?rq=1 — Moo, Sep 30 '22 at 16:11
Once again: if you just want to solve ODEs numerically, forget about RK4, use NDSolve. — xzczd, Sep 30 '22 at 16:34
Thanks Dr. Moo and Dr. xzcxd. I know NDSolve does this, but I want to do my own. — Osama, Oct 01 '22 at 08:30
Can you please run the mentioned code for the following example? x1’=x2, x1(0)=1, x2=x3, x2(0)=2, x3’=x4, x3(0)=3, x4’= −8x1+sin(t)x2−3*x3 + t^2, x4(0)=4. — Osama, Oct 01 '22 at 08:33
You need to add @xzczd in your answer or I won't get the reminder. The package in this answer of mine should be straightforward enough to use, try it out. — xzczd, Oct 01 '22 at 11:24

score 0 · Answer 1 · answered Sep 30 '22 at 16:16

Here is a simple example of a 3D system of homogeneous linear DG:

n = 3;
SeedRandom[5];
m = RandomReal[{-1, 1}, {n, n}];
f[x_, y_] := m . y;
y0 = RandomReal[{-1, 1}, {n}];
rk[x_, y_, h_] := Module[{k1, k2, k3, k4},
  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];
  y + 1/6 (k1 + 2 k2 + 2 k3 + k4)
  ]
dt = 0.1;
tmax = 1;
y = y0;
res = Reap[
    Sow[y0];
    Do[
     Sow[y = rk[x, y, dt]];
     , {x, 0, tmax, dt}]
    ][[2, 1]];
ListLinePlot[Transpose[{Range[0, tmax, dt], #}] & /@ Transpose[res]]

Runge-Kutta of order 4 to solve m first order IVPs using Mathematica

1 Answers1