Differentialgleichungen Demonstrationsübung III

Lineare Stabilitätsfunktionen, Absolute Stabilitätsbereiche

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.

Absolute stability regions for classical Runge-Kutta 4

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

Needs["DifferentialEquations`NDSolveUtilities`"];

RKamat={{1/2},{0,1/2},{0,0,1}};
RKbvec={1/6,1/3,1/3,1/6};
RKcvec = {1/2,1/2,1};
RKevec={1/6,1/3,1/3,1/6};
RKCoefficients[2, p_]:=
N[{RKamat,RKbvec,RKcvec,RKevec}, p];

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[RKamat,RKbvec,z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-3,3},{y,-3,3},Axes->True,PlotPoints->25]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]

Absolute stability regions for embedded pair Heun-Euler 2(1)

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

Needs["DifferentialEquations`NDSolveUtilities`"];

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

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[HEamat,HEbvec,z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-3,3},{y,-3,3},Axes->True,PlotPoints->25]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]

Absolute stability regions for embedded pair Bogacki-Shampine 3(2)

BS32 coefficients (a-b-c-e-numbers):

Needs["DifferentialEquations`NDSolveUtilities`"];

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

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[BSR32amat,BSR32bvec,z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-3,3},{y,-3,3},Axes->True,PlotPoints->25]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]

Absolute stability regions for embedded pair Sofroniou-Spaletta 3(2)

SS32 coefficients (a-b-c-e-numbers):

Needs["DifferentialEquations`NDSolveUtilities`"];

abce=NDSolve`EmbeddedExplicitRungeKuttaCoefficients[3,Infinity]

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[abce[[1]],abce[[2]],z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-3,3},{y,-3,3},Axes->True,PlotPoints->25]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]

Absolute stability regions for embedded pair Dormand-Prince 5(4)

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

Needs["DifferentialEquations`NDSolveUtilities`"];

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

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[DP54amat,DP54bvec,z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-6,6},{y,-6,6},Axes->True,PlotPoints->50]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]

Absolute stability regions for embedded pair Bogacki-Shampine 5(4)

BS54 coefficients (a-b-c-e-numbers):

Needs["DifferentialEquations`NDSolveUtilities`"];

abce=NDSolve`EmbeddedExplicitRungeKuttaCoefficients[5,Infinity]

Linear stability function in z:

f1=RungeKuttaLinearStabilityFunction[abce[[1]],abce[[2]],z]

RegionPlot for detecting regions of absolute stability:

RegionPlot[{Abs[f1/.{z->x+I y}]<1,Abs[f2/.{z->x+I y}]<1},{x,-6,6},{y,-6,6},Axes->True,PlotPoints->55]

Intersections with negative x-axis:

NSolve[Abs[f1]==1,z]