# Problems with Runge-Kutta methods for different functions

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[] = 0; y[] = 1; z[] = -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[]; k2y0 = k2y /. yy[]; k3y0 = k3y /. yy[];
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[]; k2z0 = k2z /. zz[]; k3z0 = k3z /. zz[];
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?

## 1 Answer

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 == 1, y' == -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 == 1, y' == -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][][t]}, {t, 0, tm},
PlotLabel -> {i, sol6[i][]}, 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[CompileGetElement[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)[],
code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[]},
Compile[{{Y0, _Real, 1}, {\[Tau], _Real}, {n, _Integer}},
Block[{Ylist}, Ylist = Table[0., {n + 1}, {Length[Y0]}];
Ylist[] = 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]] • cool! But what to do wit Sqrt, does it work for you? – Yura Holubeu May 14 '20 at 19:34
• @YuraHolubeu Please publish equation and initial condition as well for Sqrt[]. – Alex Trounev May 14 '20 at 19:41
• Given that any starting values for $y,\,y'$ corresponds to a linear combination $a\,e^{-t}+b\,e^{11t}$ and that at each step there is rounding error leading to a nonzero value for $b$, eventually one should expect the $e^{11t}$ component to dominate. I think that is the issue with the OP's direct use of an ERK method. Concerning another point, some methods let you adjust the sensitivity of the stiffness test: NDSolveValue[{eq, ic}, y, {t, 0, 4}, Method -> {"ExplicitRungeKutta", "StiffnessTest" -> {True, "SafetyFactor" -> 0.5} }]. But this is not a very stiff system to begin with. – Michael E2 May 14 '20 at 20:27
• @MichaelE2 You are right, we have the problem with stiffness here. But OP actually want handmade code to solve this problem as it reported in the previous post on mathematica.stackexchange.com/questions/221811/… – Alex Trounev May 14 '20 at 21:12
• Re Szablolcs' infinite-precision method: The RK4 method has the same eigenvectors as the exact autonomous linear system. Hence if you start on an eigensolution $y = a\,e^{-t}$, the numerical error in the RK4 method bumps you onto a nearby eigensolution $y = a\,e^{-t+t_0}$ with a constant phase shift $t_0$ at each step (for a fixed step size). Check: res[[All, 2, 1]]/Exp[-res[[All, 1]]] // Log // Differences // Differences // Simplify` – Michael E2 May 15 '20 at 3:21