RK4 residual error for ODE

Question

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[0] == 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.

Answer 1

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[0] == 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[[1]] is not zero because of a 1 ulp mistake in the initial condition:

yf[0] - 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[0] == 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@#[[2]]@"Grid" - 2, 2]},
    With[{parts1 = Range@Length@mostparts~Append~-1,
      parts2 = mostparts~Append~-1},
     Transpose@{
       Flatten[#[[1]]["Grid"][[parts1]]],
       Abs[
        #[[2]]["ValuesOnGrid"][[parts2]] -
         #[[1]]["ValuesOnGrid"][[parts1]]
        ]
       }
     ]] & /@ Partition[yseq, 2, 1],
 Frame -> True, GridLines -> {None, Automatic}
 ]

Answer 2

Is this what you want?

lowsol = NDSolve[{y'[x] == 1/(2 y[x]), y[0] == 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}]