Problems with Runge-Kutta methods for different functions

Question

Here is a great code of Runge-Kutta method proposed by @Alex Trounev in question Solving ODE with Runge-Kutta method :

n = 150;h = 5/n;

y = Table[0, n + 1]; z = Table[0, n + 1];xxx = Table[0, n + 1];
xxx[[1]] = 0; y[[1]] = 1; z[[1]] = -1;

fy[xe_, ye_, ze_] := ze;
fz[x_, yt_, zt_] := -3 yt - 4 zt;


Do[
 yy = Solve[
k1y == fy[xxx[[i]] + h*1/3, y[[i]] + 1/3*h*k1y, z[[i]]] && 
k2y == fy[xxx[[i]] + h*1, y[[i]] + 1/3*h*k1y + 2/3*h*k2y,  z[[i]]] && 
k3y == fy[xxx[[i]] + h*1, y[[i]] + 0 + 0 + 1*h*k3y, z[[i]]], {k1y,k2y, k3y}];
k1y0 = k1y /. yy[[1]]; k2y0 = k2y /. yy[[1]]; k3y0 = k3y /. yy[[1]];
 zz = Solve[
k1z == fz[xxx[[i]] + h*1/3, y[[i]], z[[i]] + 1/3*h*k1z] && 
k2z == fz[xxx[[i]] + h*1, y[[i]], z[[i]] + 1/3*h*k1z + 2/3*h*k2z] && 
k3z == fz[xxx[[i]] + h*1, y[[i]], z[[i]] + 0 + 0 + 1*h*k3z], {k1z,k2z, k3z}];

k1z0 = k1z /. zz[[1]]; k2z0 = k2z /. zz[[1]]; k3z0 = k3z /. zz[[1]];
 y[[i + 1]] = y[[i]] + h (3/4 k1y0 + 3/4 k2y0 - 1/2 k3y0);
 z[[i + 1]] = z[[i]] + h (3/4 k1z0 + 3/4 k2z0 - 1/2 k3z0);
 xxx[[i + 1]] = xxx[[i]] + h;, {i, 1, n}]

Show[Plot[Exp[-x], {x, 0, 5}, PlotRange -> All], 
 ListPlot[Transpose[{xxx, y}], PlotRange -> All, PlotStyle -> Red]] 

The problem is that it is not working for fz[x_, yt_, zt_] := 11yt +10 zt; and fz[x_, yt_, zt_] := x Sqrt[yt];... That is very sad because it is not possible to use it for real problems...

Is there a way to fix it?

Answer 1

Actually this is code @YuraHolubeu, not mine. I just made it working. Also it is not so good as it announced. Let us consider next equation with initial condition and 6 methods to solve it with NDSolve

eq = y''[t] - 10 y'[t] - 11 y[t] == 0;
ic = {y[0] == 1, y'[0] == -1};
m = {Automatic, "Adams", "BDF", "ExplicitRungeKutta", 
   "ImplicitRungeKutta", "StiffnessSwitching"};

Table[sol[i] = 
   NDSolveValue[{eq, ic}, y, {t, 0, 4}, Method -> m[[i]]], {i, 
   Length[m]}];
Table[Plot[{Exp[-t], sol[i][t]}, {t, 0, 4}, PlotLabel -> m[[i]], 
  PlotRange -> All], {i, Length[m]}]

We see that all methods failed since numerical solutions strong deviate from the exact solution $y=e^{-t}$. Not, there are no any messages from the system with warnings. Perhaps Mathematica developers should take this into account for the next version. We easily can find right solution with the next modification of this code

eq = y''[t] - 10 y'[t] - 11 y[t] == 0;
ic = {y[0] == 1, y'[0] == -1};
m = {Automatic, "Adams", "BDF", "ExplicitRungeKutta", 
  "ImplicitRungeKutta", "StiffnessSwitching"}; p = 
 Table[32, {Length[m]}];

Table[sol[i] = 
   NDSolveValue[{eq, ic}, y, {t, 0, 4}, Method -> m[[i]], 
    WorkingPrecision -> p[[i]]], {i, Length[m]}];
Table[Plot[{Exp[-t], sol[i][t]}, {t, 0, 4}, PlotLabel -> m[[i]], 
  PlotRange -> All], {i, Length[m]}] 

Now we testing method "ExplicitRungeKutta" with order from 2 to 7. We get the right numerical solution for all tests

 methods = 
  Table[{"ExplicitRungeKutta", "DifferenceOrder" -> i, 
    "StiffnessTest" -> False}, {i, 2, 7}];

Table[sol6[i] = 
   NDSolveValue[{eq, ic}, y, {t, 0, tm}, Method -> methods[[i]], 
     WorkingPrecision -> 32, MaxSteps -> 10^6] // AbsoluteTiming, {i, 
   2, 7}];


Table[Plot[{Exp[-t], sol6[i][[2]][t]}, {t, 0, tm}, 
  PlotLabel -> {i, sol6[i][[1]]}, PlotRange -> {0, 1}], {i, 2, 7}] 

Therefore we can build explicit RK method to solve this problem even for order 2. But we should calculate with high precision. I take this piece of code for RK4 implementation @Henrik Schumacher. It solves problem as it is

nsteps = 500; nsys = 2; \[Tau] = 0.01; F = 
 X \[Function] {Indexed[X, 2], 11 Indexed[X, 1] + 10 Indexed[X, 2]};
cFlow = Block[{YY, Y, k1, k2, k3, k4, \[Tau], Ylist, j}, 
   YY = Table[Compile`GetElement[Ylist, j, i], {i, 1, nsys}];
   k1 = \[Tau] F[YY];
   k2 = \[Tau] F[0.5 k1 + YY];
   k3 = \[Tau] F[0.5 k2 + YY];
   k4 = \[Tau] F[k3 + YY];
   With[{code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]], 
     code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]}, 
    Compile[{{Y0, _Real, 1}, {\[Tau], _Real}, {n, _Integer}}, 
     Block[{Ylist}, Ylist = Table[0., {n + 1}, {Length[Y0]}];
      Ylist[[1]] = Y0;
      Do[Ylist[[j + 1, 1]] = code1;
       Ylist[[j + 1, 2]] = code2;, {j, 1, n}];
      Ylist], CompilationTarget -> "WVM", RuntimeOptions -> "Speed"]]];
Ylist2 = cFlow[{-1., 1.}, \[Tau], nsteps];  

Visualization and comparison with exact solution

Show[Plot[Exp[-x], {x, 0, 5}], 
 ListPlot[Table[{i \[Tau] , Ylist2[[i, 2]]}, {i, Length[Ylist2]}], 
  PlotStyle -> Orange]]

We see that this solution exhibits the same problem as NDSolve[] with "ExplicitRungeKutta" method, because we using Compile function and it calculates with MachinePrecision only. So we take code @Szabolcs from here. With this code we get right numerical solution

ClearAll[RK4step]
RK4step[f_, h_][{t_, y_}] := Module[{k1, k2, k3, k4}, k1 = f[t, y];
  k2 = f[t + h/2, y + h k1/2];
  k3 = f[t + h/2, y + h k2/2];
  k4 = f[t + h, y + h k3];
  {t + h, y + h/6*(k1 + 2 k2 + 2 k3 + k4)}]
f[t_, {x_, v_}] := {v, 11 x + 10 v}

res = NestList[RK4step[f, 1/20], {0, {1, -1}}, 100];

Show[Plot[Exp[-x], {x, 0, 5}], 
 ListPlot[Transpose[{res[[All, 1]], res[[All, 2, 1]]}], 
  PlotStyle -> Orange]]