Differentialgleichungen Demonstrationsübung IIIb

Stiffness & Stiffness Testing

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.

Eine allgemeine Hands-On-Einführung in Mathematica mit Übungen ist im englischen Notebook 01Introduction.nb vorhanden. 

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.

Stiff equation y' = -100y+99(e^(-x))

The system y' = -100y+99(e^(-x)) with inital condition y(0) = 0 is stiff with solution y = (e^(-x))-(e^(-100x)).

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

Special Definition of differential equation:

{a,b}={1,20}; (* Interval *)
system=NDSolveProblem[{{y'[x]==-100y[x]+99Exp[-x]},{y[1]==Exp[-1]-Exp[-100]},{y},{x,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

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

Exact analytic solution:

solex = DSolve[{y'[x]==-100y[x]+99Exp[-x],y[1]==Exp[-1]-Exp[-100]},y,x]

NDSolve Illustration for classical explicit Runge-Kutta

Classical explicit Runge-Kutta coefficients (a-b-c-numbers):

ClassicalRungeKuttaCoefficients[4,p_]:=N[{{{1/2},{0,1/2},{0,0,1}},{1/6,1/3,1/3,1/6},{1/2,1/2,1}},p];

h=0.1; (* Starting step size *)

Solution with classical explicit Runge - Kutta and fixed step size:

sol1a=NDSolve[system,StartingStepSize->h, Method->{"FixedStep",Method->{"ExplicitRungeKutta","DifferenceOrder"->4,"Coefficients"->ClassicalRungeKuttaCoefficients}}]

Solution with classical explicit Euler method and fixed step size:

sol1b=NDSolve[system,StartingStepSize->h, Method->{"FixedStep",Method->{"ExplicitEuler"}}]

Readout of coordinate data for different methods and generate plots:

coords1a = First[InterpolatingFunctionCoordinates[First[y/.sol1a]]];
coords1b= First[InterpolatingFunctionCoordinates[First[y/.sol1b]]];

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

coords=coords1a;
n=20; (* Number of steps = rows in tabular data ≤ Length of interval*)
TableForm[
Take[Transpose[{coords, First[y/.sol1a][coords],First[y/.sol1b][coords],First[y/.sol][coords],First[y/.solex][coords]}],Min[(b-a)/h+1,n]],TableHeadings->{None,{"(Subscript[x, k])","Classical Runge-Kutta","ExplicitEuler","Reference","Exact"}},
TableSpacing->{0.2,2}]

NDSolve Comparisons of BSR3(2) vs SS3(2) -> Version 7

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,"Coefficients"->BSR32Coefficients,"StiffnessTest"->False,"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

Solution method for adaptive SS32 (default in Mathematica) without stiffness testing:

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

Solution method for adaptive SS32 with strict stiffness testing:

methodSS32stiff={"ExplicitRungeKutta","DifferenceOrder"->3,"StiffnessTest"->{True, "MaxRepetitions" -> {1,1}, "SafetyFactor"->1},"StepSizeControlParameters"->{1,0},"StepSizeSafetyFactors"->{1,1},"StepSizeRatioBounds"->{0,Infinity}};

NDSolve[system,Method->methodSS32stiff,AccuracyGoal->8,PrecisionGoal->8]

StepDataPlot[NDSolve[system,Method->methodSS32,MaxSteps->Infinity,AccuracyGoal->8,PrecisionGoal->8]]

StepDataPlot[NDSolve[system,Method->methodBSR32,MaxSteps->Infinity,AccuracyGoal->8,PrecisionGoal->8]]

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

{a,b}={1,20}; (* Interval *)
system1=NDSolveProblem[{{y'[x]==-100y[x]+99Exp[-x]},{y[1]==Exp[-1]-Exp[-100]},{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 Comparisons of DP5(4) vs BS5(4) -> Version 7

Dorman-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 DP54:

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

Solution method for adaptive BS54 (default in Mathematica):

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

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

{a,b}={1,20}; (* Interval *)
system1=NDSolveProblem[{{y'[x]==-100y[x]+99Exp[-x]},{y[1]==Exp[-1]-Exp[-100]},{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"}, {"DP54","BS54"}, data]

data=CompareMethods[system1,refsol,{methodDP54,methodBS54},AccuracyGoal->16,PrecisionGoal->16];
TabulateResults[{"Method","Steps","Costs","Error"}, {"DP54","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]];