Differentialgleichungen Demonstrationsübung II

Adaptive Schrittweiten: Eingebettete Runge-Kutta-Verfahren

Numerische Analysis & Computer-Algebra mit Mathematica	
Prof. Dr. Bernhard Zgraggen, HSR

Mini-Einführung Mathematica

Dieses Dokument ist ein "selbsterklärendes" Mathematica-Notebook, das interaktiv bearbeitet werden kann ("Hands-On"). Input-Zellen können mittels der Tastenkombination ⇧⌤ (in Worten: Shift-Enter) ausgeführt werden. Falls Textzellen geschlossen sind, so können sie durch Klick auf das Icon ▼ oder auf die blaue Zellklammer rechterhand geöffnet werden.

Wenn Sie in Mathematica auf Befehle oder Begriffe stossen, die Sie nicht verstehen, so steht Ihnen die mächtige Hilfe/Doc Center zur Verfügung; diese ist interaktiv.

Allfällige Aufgaben/Übungen befinden sich am Ende des Dokuments.

Generating Butcher tableaus

Needs["DifferentialEquations`NDSolveUtilities`"];

Read out abc-numbers and embedded coefficients (e) used in Mathematica:

abce=NDSolve`EmbeddedExplicitRungeKuttaCoefficients[3,Infinity];

Grid representation utility function:

butcherTableau[abce_List]:=Join[
{{0}},
Table[Flatten@{abce[[3]][[i]],abce[[1]][[i]]},{i,1,Length[abce[[3]]]}],
{PrependTo[abce[[2]],""],PrependTo[abce[[4]],""]}
];

butcherTableau[abce]//Grid//TraditionalForm//Quiet

Simulation/Coding of Step-size formulas

BSR3(2)

{a,b}={0,8};
f[x_,y_]:=1.(((1.-1./4 Cos[y]))^2); (* r.h.s. function *)

Clear[res];nn=1000;
Module[{x,y,err,h,hneu,k1,k2,k3,k4,s1,s2,n,(* res={} *),print,ag,pg},

(* Initial data *)
x=0.;y=0.;
s1=1.; s2=1.;
h=0.001;
n=1;
ag=4;pg=4;

res={{Style["x",Bold],Style["y",Bold],Style["h",Bold],Style["(Subscript[e, k])",Bold],Style["(|\!\(\*FractionBox[SubscriptBox[\(e\), \(k\)], \(\*SubscriptBox[\(ϵ\), \(a\)] + y\\ \*SubscriptBox[\(e\), \(r\)]\)]\)|)^(-1/3)",Bold],Style["(Subscript[h, neu])",Bold],Style["Status",Bold],Style["k1",Bold],Style["k2",Bold],Style["k3",Bold],Style["k4",Bold]}};

(* ... *)
While[x+h≤b&&n<nn,
k1=f[x,y];
k2=f[x+(1/2)h,y+h(1/2) k1];
k3=f[x+(3/4)h,y+h(3/4) k2];
k4=f[x+(1)h,y+h(2/9) k1+h(1/3) k2+h(4/9) k3];

(* Heun-Euler 2(1) error formula *)
err=h((2/9-(-5/72))k1+(1/3-1/12)k2+(4/9-1/9)k3+(0-(-1/8))k4);

hneu=h s1 (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3));

(* Schrittwiederholung: wenn err > tolerance *)
If[(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3))≥1,
print="Proceed";
res=Append[res,{x,y,h,err,(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3)),hneu,print,k1,k2,k3,k4}] ;
x=x+h;
y=y+h(2/9) k1+h(1/3) k2+h(4/9) k3,
print="Reject";
res=Append[res,{x,y,h,err,(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3)),hneu,print,k1,k2,k3,k4}] ;
];
h=hneu;
n++;

]
]

res//Length
TableForm[Take[SetPrecision[res,10],Min[Length[res],nn,100]]]

ListPlot[DeleteDuplicates@Transpose@(Transpose[res][[{1,2}]]),PlotStyle->Red]

ListLinePlot[Transpose[res][[3]],PlotRange->All]

SS3(2)

{a,b}={0,8};
f[x_,y_]:=1.(((1.-1./4 Cos[y]))^2); (* 10Cos[x+y] *); (* r.h.s. function *)

Clear[res];nn=1000;
Module[{x,y,err,h,hneu,k1,k2,k3,k4,s1,s2,n,(* res={} *),print,ag,pg},

(* Initial data *)

x=0.;y=0. (* Pi/2. *);
s1=1.; s2=1.;
h=0.001;
n=1;
ag=4;pg=4;

res={{Style["x",Bold],Style["y",Bold],Style["h",Bold],Style["(Subscript[e, k])",Bold],Style["(|\!\(\*FractionBox[SubscriptBox[\(e\), \(k\)], \(\*SubscriptBox[\(ϵ\), \(a\)] + y\\ \*SubscriptBox[\(e\), \(r\)]\)]\)|)^(-1/3)",Bold],Style["(Subscript[h, neu])",Bold],Style["Status",Bold],Style["k1",Bold],Style["k2",Bold],Style["k3",Bold],Style["k4",Bold]}};

(* ... *)
While[x+h≤b&&n<nn,
k1=f[x,y];
k2=f[x+(1/2)h,y+h(1/2) k1];
k3=f[x+(1)h,y+-h k1+h(2) k2];
k4=f[x+(1)h,y+h(1/6) k1+h(2/3) k2+h(1/6) k3];

(* ... *)
err=h((1/6-(1/72) (-10+(Sqrt[82])))k1+(2/3-(1/36) (10-(Sqrt[82])))k2+(1/6-(1/144) (28-(Sqrt[82])))k3+(0-(1/48) (-16+(Sqrt[82])))k4);

hneu=h s1 (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3));

(* Ev. Schrittwiederholung: wenn hneu < h *)
If[ (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3))≥1,
print="Proceed";
res=Append[res,{x,y,h,err,(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3)),hneu,print,k1,k2,k3,k4}] ;
x=x+h;
y=y+h(1/6) k1+h(2/3) k2+h(1/9) k3,
print="Reject";
res=Append[res,{x,y,h,err,(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/3)),hneu,print,k1,k2,k3,k4}] ;
];
h=hneu;
n++;

]
]

res//Length
TableForm[Take[SetPrecision[res,10],Min[Length[res],nn,100]]]

ListPlot[DeleteDuplicates@Transpose@(Transpose[res][[{1,2}]]),PlotStyle->Red]

ListLinePlot[Transpose[res][[3]],PlotRange->All]

HE2(1)

{a,b}={0.,8};
f[x_,y_]:=1.(((1.-1./4 Cos[y]))^2); (* r.h.s. function *)

Clear[res];nn=200;
Module[{x,y,err,h,hneu,k1,k2,s1,s2,n,(* res={} *),print,ag,pg},
(* Initial data *)
x=0.;y=0.;
s1=1.; s2=1.;
h=0.001;
n=1;
ag=4;pg=4;

res={{Style["x",Bold],Style["y",Bold],Style["h",Bold],Style["(Subscript[e, k])",Bold],Style["(|\!\(\*FractionBox[SubscriptBox[\(e\), \(k\)], \(\*SubscriptBox[\(ϵ\), \(a\)] + y\\ \*SubscriptBox[\(e\), \(r\)]\)]\)|)^(-1/2)",Bold],Style["(Subscript[h, neu])",Bold],Style["Status",Bold],Style["k1",Bold],Style["k2",Bold]}};

(* ... *)
While[x≤b+h &&n<nn,
k1=f[x,y];
k2=f[x+h,y+h k1];

(* Heun-Euler 2(1) error formula *)
err=h((1/2-1)k1+(1/2-0)k2);

hneu=h s1 (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/2));

(* Ev. Schrittwiederholung: wenn err > s2*tol *)

If[ (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/2))>=1.,
print="Proceed";
res=Append[res,{x,y,h,err,(((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/2)),hneu,print,k1,k2}] ;
x=x+h;
y=y+h(1/2) k1+h(1/2) k2,

print="Reject";
res=Append[res,{x,y,h,err, (((Abs[((err/(s2((10^(-ag))+ (10^(-pg))Abs[y]))))]))^(-1/2)),hneu,print,k1,k2}] ;
];

h=hneu;
n++;
]
]

res//Length
TableForm[Take[SetPrecision[res,10],Min[Length[res],nn,25]]]

ListPlot[DeleteDuplicates@Transpose@(Transpose[res][[{1,2}]]),PlotStyle->Red]

ListLinePlot[Transpose[res][[3]],PlotRange->All]

HE2(1) multi-variate

{a,b}={0,100};
f[x_,{y_,v_}]:={v,0.2(1-y^2)v-y}; (* r.h.s. function van der Pol *)

Clear[res]; nn=1000;
Module[{x,y,err,h,hneu,k1,k2,s1,s2,n,(* res={} *),print,ag,pg},
(* Initial data *)
x=0.;y={1,-1};
s1=1.; s2=1.;
h=0.001;
n=1;
res={};
ag=4;pg=4;
res={{Style["x",Bold],Style["y",Bold],Style["h",Bold],Style["(Subscript[e, k])",Bold],Style["\[LeftDoubleBracketingBar]\!\(\*FractionBox[SubscriptBox[\(e\), \(k\)], \(\*SubscriptBox[\(ϵ\), \(a\)] + y\\ \*SubscriptBox[\(e\), \(r\)]\)]\)\!\(\*SuperscriptBox[\(\[RightDoubleBracketingBar]\), \(\(\(-1\)/2\)\(\\n\)
\)]\)",Bold],Style["(Subscript[h, neu])",Bold],Style["Status",Bold],Style["k1",Bold],Style["k2",Bold]}};

(* ... *)
While[x≤b&&n<nn,
k1=f[x,y];
k2=f[x+h,y+h k1];

(* Heun-Euler 2(1) error formula *)
err=h((1/2-1)k1+(1/2-0)k2);

hneu=h s1 (((Norm[{((err[[1]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[2]]])))},Infinity]))^(-1/2));

(* Ev. Schrittwiederholung: wenn err > s2*tol *)

If[ (((Norm[{((err[[1]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[2]]])))},Infinity]))^(-1/2))>=1.,
print="Proceed";
 res=Append[res,{x,y (* StringForm["(``)",Column[y,", "]] *),h,StringForm["(``)",Column[err,", "]],(((Norm[{((err[[1]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[2]]])))},Infinity]))^(-1/2)),hneu,print,StringForm["(``)",Column[k1,", "]],StringForm["(``)",Column[k2,", "]]}] ; 
(* res=Append[res,{x,y,h,StringForm["(``)",Row[err,", "]],(((Norm[{((err[[1]])/(s2((10^(-1))+ (10^(-2))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-1))+ (10^(-2))Abs[y[[2]]])))},Infinity]))^(-1/2)),hneu,print}] ; *)

x=x+h;
y=y+h(1/2) k1+h(1/2) k2,
print="Reject";
 res=Append[res,{x,y (* StringForm["(``)",Column[y,", "]] *),h,StringForm["(``)",Column[err,", "]],(((Norm[{((err[[1]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-ag))+ (10^(-pg))Abs[y[[2]]])))},Infinity]))^(-1/2)),hneu,print,StringForm["(``)",Column[k1,", "]],StringForm["(``)",Column[k2,", "]]}] ; 
(* res=Append[res,{x,y,h,StringForm["(``)",Row[err,", "]],(((Norm[{((err[[1]])/(s2((10^(-1))+ (10^(-2))Abs[y[[1]]]))),((err[[2]])/(s2((10^(-1))+ (10^(-2))Abs[y[[2]]])))},Infinity]))^(-1/2)),hneu,print}] *)
]; 

h=hneu;
n++;
]
]

Length[res]
TableForm[Take[SetPrecision[res,10],Min[Length[res],100]]]

ListPlot[DeleteDuplicates@Transpose@(Transpose[res][[{2}]]),PlotStyle->Red]

ListLinePlot[Transpose[res][[3]],PlotRange->All]

NDSolve: Illustrations for embedded pair Heun-Euler 2(1) -> Version 7

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system,sol]]

Heun-Euler coefficients (a-b-c-e-numbers):

HEamat={{1}};
HEbvec={1/2,1/2};
HEcvec = {1};
HEevec={1,0};
HECoefficients[2, p_]:=
N[{HEamat,HEbvec,HEcvec,HEevec}, p];

HeunCoefficients[2, p_]:=
N[{HEamat,HEbvec,HEcvec}, p];

Solution method for adaptive Heun-Euler 2(1):

methodHE={"ExplicitRungeKutta","DifferenceOrder"->2,"EmbeddedDifferenceOrder"->1,"Coefficients"->HECoefficients,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

solHE=NDSolve[system,StartingStepSize->0.001 ,Method->methodHE,AccuracyGoal->4,PrecisionGoal->4 ,MaxSteps->Infinity,MaxStepSize->Infinity,MaxStepFraction->Infinity]

Solution method for non-adaptive classical Heun method:

methodH={"FixedStep",Method->{"ExplicitRungeKutta","DifferenceOrder"->2,"Coefficients"->HeunCoefficients,"StiffnessTest"->False}};

Solution with classical explicit Heun method (order 2) and fixed step size:

hstart=0.1;
solH=NDSolve[system,StartingStepSize->hstart, Method->methodH]

Readout of coordinate data for different methods and generate plots:

coordsHE = First[InterpolatingFunctionCoordinates[First[y/.solHE]]];
gsolHE=
ListPlot[Transpose[{coordsHE, First[y/.solHE][coordsHE]}],PlotStyle->Gray];


coordsH = First[InterpolatingFunctionCoordinates[First[y/.solH]]];
gsolH=ListPlot[Transpose[{coordsH, First[y/.solH][coordsH]}],PlotStyle->Red];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],Plot[y[x]/.solHE,{x,a,b},PlotStyle->{Dotted,Darker@Gray}],gsolHE,gsolH]

Read out coordinate data for different methods and generate tabular data:

coords=coordsHE;
n={1,20}; (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[{coords, First[y/.solHE][coords],First[y/.solH][coords],First[y/.sol][coords]}],n],TableHeadings->{None,{"(Subscript[t, k])","Heun-Euler 2(1)","Heun","Reference"}},
TableSpacing->{0.2,2}]

Plot the step-sizes as a function:

StepDataPlot[solHE]

Comparison of methods (utility)

Compare the methods (maybe long running time). For this the system unknown y has to be replaced by y[x]:

{a,b}={0,8}; (* Interval *)
system1=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y[x]},{x,a,b},{},{},{}}]

sol = NDSolve[system1,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system1,sol]]

data=CompareMethods[system1,refsol,{methodHE,methodH},MaxStepSize->Infinity,MaxStepFraction->Infinity,MaxSteps->Infinity];

TabulateResults[{"Method","Steps","Costs","Error"}, {"HE","Heun"}, data]

NDSolve: Illustrations for embedded pair Bogacki-Shampine 3(2) -> Version 7

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system,sol]]

Bogacki-Shampine-Ralston BSR coefficients (a-b-c-e-numbers):

BSR32amat={{1/2},{0,3/4},{2/9,1/3,4/9}};
BSR32bvec={2/9,1/3,4/9,0};
BSR32cvec = {1/2,3/4,1};
BSR32evec={-5/72,1/12,1/9,-1/8};
BSR32Coefficients[3, p_]:=
N[{BSR32amat,BSR32bvec,BSR32cvec,BSR32evec}, p];

Solution method for adaptive BSR32:

methodBSR32={"ExplicitRungeKutta","DifferenceOrder"->3,"EmbeddedDifferenceOrder"->2,"Coefficients"->BSR32Coefficients,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

solBSR32=NDSolve[system,StartingStepSize->0.001,AccuracyGoal->4,PrecisionGoal->4,MaxSteps->Infinity,Method->methodBSR32]

Readout of coordinate data for different methods and generate plots:

coordsBSR32 = First[InterpolatingFunctionCoordinates[First[y/.solBSR32]]];
gsolBSR32=
ListPlot[Transpose[{coordsBSR32, First[y/.solBSR32][coordsBSR32]}],PlotStyle->Red];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],Plot[y[x]/.solBSR32,{x,a,b},PlotStyle->{Dotted,Darker@Gray}],gsolBSR32]

Read out coordinate data for different methods and generate tabular data:

coords=coordsBSR32;
n={1,Length[coords]} (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[{coords, First[y/.solBSR32][coords],First[y/.sol][coords]}],n],TableHeadings->{None,{"(Subscript[t, k])","Bogacki-Shampine-Ralston 3(2)","Reference"}},
TableSpacing->{0.2,2}]

Plot the step-sizes as a function:

StepDataPlot[solBSR32]
StepDataPlot[sol]

NDSolve: Illustrations for default embedded pair Sofroniou-Spaletta 3(2) and comparison -> Version 7

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];

Sofroniou-Spalette 3(2) coefficients (a-b-c-e-numbers):

NDSolve`EmbeddedExplicitRungeKuttaCoefficients[3,Infinity]

Solution method for adaptive SS32:

methodSS32={"ExplicitRungeKutta","DifferenceOrder"->3,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

solSS32=NDSolve[system,StartingStepSize->Automatic,AccuracyGoal->4,PrecisionGoal->4,MaxSteps->Infinity,Method->methodSS32]

Readout of coordinate data for different methods and generate plots:

coordsSS32 = First[InterpolatingFunctionCoordinates[First[y/.solSS32]]];
gsolSS32=
ListPlot[Transpose[{coordsSS32, First[y/.solSS32][coordsSS32]}],PlotStyle->Red];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],Plot[y[x]/.solSS32,{x,a,b},PlotStyle->{Dotted,Darker@Gray}],gsolSS32]

Read out coordinate data for different methods and generate tabular data:

coords=coordsSS32;
n={1,Length[coords]} (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[{coords, First[y/.solSS32][coords],First[y/.sol][coords]}],n],
TableHeadings->{None,{"(Subscript[t, k])","Sofroniou-Spaletta 3(2)","Reference"}},
TableSpacing->{0.2,2}]

Plot the step-sizes as a function:

StepDataPlot[solSS32]

Comparison of methods (utility function)

Compare methods (maybe long running time). Therefore the system has to be (re-)defined with y[x] as unknown.

{a,b}={0,8}; (* Interval *)
system1=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y[x]},{x,a,b},{},{},{}}]

Reference solution:

sol = NDSolve[system1,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system1,sol]]

data=CompareMethods[system1,refsol,{methodBSR32,methodSS32},AccuracyGoal->4,PrecisionGoal->4];

TabulateResults[{"Method","Steps","Costs","Error"}, {"BSR32","SS32"}, data]

data=CompareMethods[system1,refsol,{methodBSR32,methodSS32},AccuracyGoal->8,PrecisionGoal->8];
TabulateResults[{"Method","Steps","Costs","Error"}, {"BSR32","SS32"}, data]

data=CompareMethods[system1,refsol,{methodBSR32,methodSS32},AccuracyGoal->16,PrecisionGoal->16];
TabulateResults[{"Method","Steps","Costs","Error"}, {"BSR32","SS32"}, data]

NDSolve: Illustrations for embedded pair Dormand-Prince 5(4) -> Version 7

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system,sol]]

Dormand-Prince DP54 coefficients (a-b-c-e-numbers):

DP54amat={
{1/5},
{3/40,9/40},
{44/45,-56/15,32/9},
{19372/6561,-25360/2187,64448/6561,-212/729},
{9017/3168,-355/33,46732/5247,49/176,-5103/18656},
{35/384,0,500/1113,125/192,-2187/6784,11/84}};
DP54bvec={35/384,0,500/1113,125/192,-2187/6784,11/84, 0};
DP54cvec = {1/5,3/10,4/5,8/9,1, 1};
DP54evec={71/57600,0,-71/16695,71/1920,-17253/339200,22/525,-1/40};
DP54Coefficients[5, p_]:=
N[{DP54amat,DP54bvec,DP54cvec,DP54evec}, p];

Solution method for adaptive BSR32:

methodDP54={"ExplicitRungeKutta","DifferenceOrder"->5,"Coefficients"->DP54Coefficients,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

solDP54=NDSolve[system,StartingStepSize->Automatic,AccuracyGoal->8,PrecisionGoal->8,MaxSteps->Infinity,Method->methodDP54]

Readout of coordinate data for different methods and generate plots:

coordsDP54 = First[InterpolatingFunctionCoordinates[First[y/.solDP54]]];
gsolDP54=
ListPlot[Transpose[{coordsDP54, First[y/.solDP54][coordsDP54]}],PlotStyle->Red];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],Plot[y[x]/.solDP54,{x,a,b},PlotStyle->{Dotted,Darker@Gray}],gsolDP54]

Read out coordinate data for different methods and generate tabular data:

coords=coordsDP54;
n={1,Length[coords]} (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[SetPrecision[{coords, First[y/.solDP54][coords],First[y/.sol][coords]},8+2]],n],TableHeadings->{None,{"(Subscript[t, k])","Dormand-Prince 5(4)","Reference"}},
TableSpacing->{0.2,2}]

Plot the step-sizes as a function:

StepDataPlot[solDP54]
StepDataPlot[sol]

NDSolve: Illustrations for default embedded pair Bogacki-Shampine 5(4) and comparison -> Version 7

Special Definition of differential equation:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Definition of differential equation:

{a,b}={0,8}; (* Interval *)
system=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

sol = NDSolve[system,Method->"StiffnessSwitching",WorkingPrecision->32];

Bogacki-Shampine 5(4) coefficients (a-b-c-e-numbers):

NDSolve`EmbeddedExplicitRungeKuttaCoefficients[5,Infinity]

Solution method for adaptive SS32:

methodBS54={"ExplicitRungeKutta","DifferenceOrder"->5,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

solBS54=NDSolve[system,StartingStepSize->Automatic,AccuracyGoal->8,PrecisionGoal->8,MaxSteps->Infinity,Method->methodBS54]

Readout of coordinate data for different methods and generate plots:

coordsBS54 = First[InterpolatingFunctionCoordinates[First[y/.solBS54]]];
gsolBS54=
ListPlot[Transpose[{coordsBS54, First[y/.solBS54][coordsBS54]}],PlotStyle->Red];

Show the plots inluding reference solution:

Show[Plot[y[x]/.sol,{x,a,b},PlotStyle->{Black,Dashed}],Plot[y[x]/.solBS54,{x,a,b},PlotStyle->{Dotted,Darker@Gray}],gsolBS54]

Read out coordinate data for different methods and generate tabular data:

coords=coordsBS54;
n={1,Length[coords]} (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[SetPrecision[{coords, First[y/.solBS54][coords],First[y/.sol][coords]},8+2]],n],
TableHeadings->{None,{"(Subscript[t, k])","Bogacki-Shampine 5(4)","Reference"}},
TableSpacing->{0.2,2}]

Plot the step-sizes as a function:

StepDataPlot[solBS54]

Comparison of methods (utility function)

Compare methods (maybe long running time). Therefore the system has to be (re-)defined with y[x] as unknown.

{a,b}={0,8}; (* Interval *)
system1=NDSolveProblem[{{y'[x]==1(((1-1/4 Cos[y[x]]))^2)},{y[0]==0},{y[x]},{x,a,b},{},{},{}}]

Reference solution:

sol = NDSolve[system1,Method->"StiffnessSwitching",WorkingPrecision->32];
refsol = First[FinalSolutions[system1,sol]]

data=CompareMethods[system1,refsol,{methodDP54,methodBS54},AccuracyGoal->8,PrecisionGoal->8];

TabulateResults[{"Method","Steps","Costs","Error"}, {"DP5432","BS54"}, data]

$MachinePrecision

data=CompareMethods[system1,refsol,{methodDP54,methodBS54},AccuracyGoal->16,PrecisionGoal->16];
TabulateResults[{"Method","Steps","Costs","Error"}, {"DP5432","BS54"}, data]

Anhang: Demo-Codes / Codes

TabulateResults[labels_List,names_List,data_List]:=
DisplayForm[
FrameBox[
GridBox[
Apply[{labels,##}&,MapThread[Prepend,{data,names}]],
RowLines->True,ColumnLines->True
]
]
] /; SameQ[Length[names], Length[data]];