# 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}]? Feb 8, 2020 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}]; Feb 8, 2020 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? Feb 8, 2020 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}] Feb 8, 2020 at 14:18
• Dear Michael, RK4 with fixed step size 0.1, I think all other steps like residual, plot...depend on it, right? Feb 8, 2020 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, Feb 8, 2020 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. Feb 7, 2020 at 21:01