1
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Maybe LogPlot[Evaluate[Abs[{residual[x] /. y -> yf}]], {x, 0, 1}, PlotStyle -> {GrayLevel[0]}, AxesOrigin -> {0, 0}]? $\endgroup$ – Michael E2 Feb 8 at 3:34
  • $\begingroup$ 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}]; $\endgroup$ – user62716 Feb 8 at 7:00
  • $\begingroup$ 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? $\endgroup$ – Michael E2 Feb 8 at 14:10
  • $\begingroup$ 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[0] == 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}] $\endgroup$ – user62716 Feb 8 at 14:18
  • $\begingroup$ Dear Michael, RK4 with fixed step size 0.1, I think all other steps like residual, plot...depend on it, right? $\endgroup$ – user62716 Feb 8 at 14:34
4
$\begingroup$

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
 ]

enter image description here

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}
 ]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Many thanks Michael, $\endgroup$ – user62716 Feb 8 at 20:33
0
$\begingroup$

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}]

Log plot of residuals for selected points.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Many thanks for your response, please I modified my question. $\endgroup$ – user62716 Feb 7 at 21:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.