# How to write more simple variable-step Runge Kutta algorithm to solve this differential equation?

Here is a case of Runge Kutta algorithm with variable-step size:

ClearAll["Global*"]
x[0] = 0;
y[0] = 1.;
ε = 0.01;(*Maximum allowable error*)
L = 0.1;(*Initial step size*)
f[x_, y_] := x^2
SingleRK[u_, v_, L_ : 0.1] :=
Module[{x = u, y = v, h = L, K1, K2, K3, K4, ynext, xnext},
K1 = f[x, y];
K2 = f[x + h/2, y + h*K1/2];
K3 = f[x + h/2, y + h*K2/2];
K4 = f[x + h, y + h*K3];
ynext = y + h/6 (K1 + 2 K2 + 2 K3 + K4);
xnext = x + h;
Return[{xnext, ynext}]]
LocalRK[{u_, v_}] :=
Module[{xnext, xinitial, ynext, yinitial, x, y, h = L}, x = u; y = v;
{xnext, ynext} = SingleRK[x, y, h];
If[Abs[ynext - y] < ε,
While[Abs[ynext - y] < ε, h = 2 h;
ynext = SingleRK[x, y, h] // Last;],
While[Abs[ynext - y] > ε, h = h/2;
ynext = SingleRK[x, y, h] // Last;]];
{x + h, ynext}]
data = NestList[LocalRK, {x[0], y[0]}, 35];
Abs /@ Subtract @@@
Partition[data[[;; , 1]], 2, 1](*Step-size information*)
Show[Plot[1/3 (x^3 + 3), {x, 0, 1}, PlotStyle -> {Opacity[0.5], Red}],
ListPlot[data, AxesOrigin -> {0, 0},
PlotStyle -> {Black, PointSize[0.01]}]]
data


But I found that the error of the first step was obviously more than 0.01, but it did not automatically adjust to the appropriate step size. How to improve this problem(I need the step size of each step to ensure that the error is exactly less than 0.01 in the maximum step size)?

Then I need more concise examples of variable-step algorithm for solving definite integral or differential equations.

I hope to get more adaptive step size algorithms to solve the following differential equation:

DSolve[{y'[x] == x^2, y[0] == 1}, y[x], x] // FullSimplify

• Why not use NDSolve? Commented Oct 6, 2020 at 13:13
• @ChrisK I want to use a specific numerical algorithm to solve this problem. This is an exercise in numerical analysis. Commented Oct 7, 2020 at 0:35
• It is also possible to define user methods for NDSolve, see tutorial/NDSolvePlugIns
– I.M.
Commented Oct 7, 2020 at 0:58
• @I.M. Thank you for your comments. I need to use specific numerical algorithms to solve differential equations. This is required for homework. Commented Oct 7, 2020 at 2:02

In order to make the error of the first step less than 0.01, we can solve this problem with the help of the intermediate variable sol:

ClearAll["Global*"]
x[0] = 0;
y[0] = 1.;
ε = 0.01;
L = 0.1;
f[x_, y_] := x^2
SingleRK[u_, v_, L_ : 0.1] :=
Module[{x = u, y = v, h = L, K1, K2, K3, K4, ynext, xnext},
K1 = f[x, y];
K2 = f[x + h/2, y + h*K1/2];
K3 = f[x + h/2, y + h*K2/2];
K4 = f[x + h, y + h*K3];
ynext = y + h/6 (K1 + 2 K2 + 2 K3 + K4);
xnext = x + h;
Return[{xnext, ynext}]]
LocalRK[{u_, v_}] :=
Module[{xnext, xinitial, ynext, yinitial, x, y, h = L, sol}, x = u; y = v;
{xnext, ynext} = SingleRK[x, y, h];
If[Abs[ynext - y] < ε,
While[Abs[ynext - y] < ε, h = 2 h;
ynext = SingleRK[x, y, h] // Last; sol = SingleRK[x, y, h/2];],
While[Abs[ynext - y] > ε, h = h/2;
ynext = SingleRK[x, y, h] // Last; sol = SingleRK[x, y, h];]];
sol]
data = NestList[LocalRK, {x[0], y[0]}, 35];
Abs /@ Subtract @@@
Partition[data[[;; , 1]], 2, 1](*Step-size information*)
Show[Plot[1/3 (x^3 + 3), {x, 0, 1}, PlotStyle -> {Opacity[0.5], Red}],
ListPlot[data, AxesOrigin -> {0, 0},
PlotStyle -> {Black, PointSize[0.01]}]]
data


But the logic of this algorithm is not very clear. I hope community members can provide more and more sophisticated algorithms.

Runge-Kutta-Felberg variable-step method:

RKF[f_, x0_, y0_, xmax_, hmax_, hmin_, eps_] :=
Module[{x = x0, y = y0, pts = {{x0, y0}}, K1, K2, K3, K4, K5, K6,
delta, h = hmax, flag = 1, R},
While[flag == 1, K1 = h*f[x, y]; K2 = h*f[x + h/4, y + 1/4*K1];
K3 = h*f[x + 3 h/8, y + 3/32*K1 + 9/32*K2];
K4 = h*f[x + 12 h/13,
y + 1932/2197*K1 - 7200/2197*K2 + 7296/2197*K3];
K5 = h*f[x + h, y + 439/216*K1 - 8*K2 + 3680/513 K3 - 845/4104 K4];
K6 = h*f[x + h/2,
y - 8/27*K1 + 2*K2 - 3544/2565 K3 + 1859/4104 K4 - 11/40*K5];
R = Abs[K1/360 - 128/4275 K3 - 2197/75240 K4 + K5/50 + 2/55 K6]/h;
If[R <= eps,
y = y + (16/135 K1 + 6656/12825 K3 + 28561/56430 K4 - 9/50 K5 +
2/55 K6); x = x + h; AppendTo[pts, {x, y}]];
delta = 0.84 (eps/R)^(1/4);
Which[delta <= 0.1, h = 0.1 h, delta >= 4, h = 4 h, True,
h = h*delta]; If[h > hmax, h = hmax];
If[x >= xmax, flag = 0,
If[x + h > xmax, h = xmax - x,
If[h < hmin, flag = 0; Print["minimum h exceeded"]]]]]; pts]

f[x_, y_] := y^2
RKF[f, 0, 1, 2, 0.1, 0.01, 0.01]
Differences[%[[All, 1]]](*Step-size information*)

• All this is doable in NDSolve as @I.M. pointed out, though IIRC the step control depends on the order when it's not an embedded method. In fact, for RK methods, you don't even need the plug-in framework. Commented Oct 8, 2020 at 2:49