# RK4 residual error for ODE

Below I am giving the code that I am using

ClassicalRungeKuttaCoefficients[4, prec_] :=
With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}},
N[{amat, bvec, cvec}, prec]]

{yf} = {y} /.
First@NDSolve[{y'[x] == 1/(2 y[x]), y == 1/10}, {y}, {x, 0, 1},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients},
StartingStepSize -> 0.1];

residual[x_] = y'[x] - 1/(2 y[x]);

er1 = Evaluate[Abs[{residual[x] /. y -> yf}]];

LogPlot[er1, {x, 0, 1}]


I need to get the residual and plot for only at the values of x between 0 1 with step 0.1, i.e. x = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1};

Thanks.

• Maybe LogPlot[Evaluate[Abs[{residual[x] /. y -> yf}]], {x, 0, 1}, PlotStyle -> {GrayLevel}, AxesOrigin -> {0, 0}]? – Michael E2 Feb 8 '20 at 3:34
• Dear Michael, many thanks for your very useful comment, please can we make the plot for just the values of x between [0,1] step 0.1?x = Table[i, {i, 0, 1, 0.1}]; – user62716 Feb 8 '20 at 7:00
• Did you want the step size taken by NDSolve to be fixed at 0.1 as well, or do you want NDSolve to adapt the step size to reduce the error? – Michael E2 Feb 8 '20 at 14:10
• Dear Michael, I need the new code using Runge-kutta of order 4, provides me similar output and plot like the previous codes. lowsol = NDSolve[{y'[x] == 1/(2 y[x]), y == 1/10}, y, {x, 0, 1}, InterpolationOrder -> All]; residual[x_] = y'[x] - 1/(2 y[x]); x = Table[i, {i, 0, 1, 0.1}]; er1 = Abs[residual[x] /. lowsol]; z1 = Max[er1] ListLogPlot[er1, Joined -> True, PlotRange -> All, Frame -> True, Axes -> True, PlotMarkers -> {Automatic, 15}] – user62716 Feb 8 '20 at 14:18
• Dear Michael, RK4 with fixed step size 0.1, I think all other steps like residual, plot...depend on it, right? – user62716 Feb 8 '20 at 14:34

My understanding is that the OP wants the "FixedStep" method with the RK4 integration method. I usually add the option MaxStepFraction -> 1, even though it's not strictly necessary in this particular problem. But sometimes you want to alter the interval or step size -- and how many remember the default setting, anyway? -- and it's more convenient to turn off this check when you start.

ClassicalRungeKuttaCoefficients[4, prec_] :=
With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}},
N[{amat, bvec, cvec}, prec]]

{yf} = {y} /.
First@NDSolve[{y'[x] == 1/(2 y[x]), y == 1/10}, {y}, {x, 0, 1},
Method -> {"FixedStep",
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients}},
StartingStepSize -> 0.1, MaxStepFraction -> 1,
InterpolationOrder -> All];


Check the steps taken (the "Grid"), which are as desired:

yf@"Grid"
(*  {{0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6}, {0.7}, {0.8}, {0.9}, {1.}}  *)


Replace y[x] and y'[x] by their values at the steps:

rr = residual[x] /. {y[x] -> yf@"ValuesOnGrid",
y'[x] -> yf'@"ValuesOnGrid"}
(*
{-8.88178*10^-16, 0., 2.22045*10^-16, 0., 0.,
2.22045*10^-16, -3.33067*10^-16, 2.22045*10^-16,
2.22045*10^-16,  1.11022*10^-16, 0.}
*)


The point to be observed is that the residuals are roughly zero up to round-off error. This is because the residual equation is used to compute the value of y'[x]; so there is no approximation/truncation error, only machine round-off error.

Remark: The first error rr[] is not zero because of a 1 ulp mistake in the initial condition:

yf - 0.1
(*  -2.77556*10^-17  *)


I guess that's a bug. You can get the same error with something like (0.1 + 0.5 - 0.5) - 0.1, but I don't see a need to shift the initial condition.

A plot, especially a log plot, is going to be disappointing, but since it was asked for several times, here it is:

ListLogPlot[
yf@"Coordinates" ~Append~ Abs@rr // Transpose
] Some points are missing because Log[0.] is Indeterminate.

If you want to estimate the error, measure the difference between two approximations, either with different step sizes $$h$$ and $$h/2$$ or with different orders. One problem we run into with NDSolve is that round-off error means that it is possible that step size $$h/2$$ does not have exactly twice as many steps as the number with step size $$h$$. If you pick a step size $$h=2^{-n}$$, then there won't be such a problem, but in the OP's example, the step size is $$h=0.1$$. This causes some complications in the code, if we want to illustrate the convergence of the OP's problem as $$h \rightarrow 0$$:

yseq = Table[
NDSolveValue[
{y'[x] == 1/(2 y[x]), y == 1/10}, y, {x, 0, 1},
Method -> {"FixedStep",
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients}},
StartingStepSize -> 1/10*2^-n, MaxStepFraction -> 1,
InterpolationOrder -> All],
{n, 0, 10}];

ListLogPlot[
With[{mostparts = Range[1, Length@#[]@"Grid" - 2, 2]},
With[{parts1 = Range@Length@mostparts~Append~-1,
parts2 = mostparts~Append~-1},
Transpose@{
Flatten[#[]["Grid"][[parts1]]],
Abs[
#[]["ValuesOnGrid"][[parts2]] -
#[]["ValuesOnGrid"][[parts1]]
]
}
]] & /@ Partition[yseq, 2, 1],
Frame -> True, GridLines -> {None, Automatic}
] • Many thanks Michael, – user62716 Feb 8 '20 at 20:33

Is this what you want?

lowsol = NDSolve[{y'[x] == 1/(2 y[x]), y == 1/10}, y, {x, 0, 1},
Method -> "ExplicitRungeKutta"]


pts = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1};

ListLogPlot[Transpose[{pts, Abs[residual@pts]}], PlotStyle -> Black,
AxesOrigin -> {0, 0}, PlotRange -> {-0.02, 0.03}] • Many thanks for your response, please I modified my question. – user62716 Feb 7 '20 at 21:01