Differentialgleichungen Demonstrationsübung IV

Systeme und Differentialgleichungen höherer Ordnung

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.

Van der Pol's equation:  z''(t)=μ(1-(z^2))z'-z

demovanderPol1;

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

Special Definition of differential equation:

{a,b}={0,100}; (* Interval *)
μ=2/10;
z0=1;v0=-1;
system=NDSolveProblem[{{z'[t]==v[t],v'[t]==μ(1-((z[t])^2))v[t]-z[t]},{z[0]==z0,v[0]==v0},{z,v},{t,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

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

ParametricPlot[Evaluate[{z[t],v[t]}/.sol],{t,a,b}]

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

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

Solution function with adaptive Heun-Euler 2(1) and low accuracy and precision:

solHE=NDSolve[system,StartingStepSize->Automatic,AccuracyGoal->1,PrecisionGoal->2,MaxSteps->Infinity,MaxStepSize->Infinity,MaxStepFraction->Infinity,Method->methodHE,NormFunction->Infinity]

ParametricPlot[Evaluate[{z[t],v[t]}/.solHE],{t,a,b}]

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

coordsHE = First[InterpolatingFunctionCoordinates[First[z/.solHE]]];
Length[coordsHE]

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

Plot the step-sizes as a function:

StepDataPlot[solHE]

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

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 Dormand-Prince 5(4):

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

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

ParametricPlot[Evaluate[{z[t],v[t]}/.solDP54],{t,a,b}]

coordsDP54 = First[InterpolatingFunctionCoordinates[First[z/.solDP54]]];
Length[coordsDP54]

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

StepDataPlot[solDP54]

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

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

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

ParametricPlot[Evaluate[{z[t],v[t]}/.solBS54],{t,a,b}]

coordsBS54 = First[InterpolatingFunctionCoordinates[First[z/.solBS54]]];
Length[coordsBS54]

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

StepDataPlot[solBS54]

Comparison of methods (utility) -> Version 7

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

{a,b}={0,100}; (* Interval *)
μ=2/10;
z0=1;v0=-1;
system1=NDSolveProblem[{{z'[t]==v[t],v'[t]==μ(1-((z[t])^2))v[t]-z[t]},{z[0]==z0,v[0]==v0},{z[t],v[t]},{t,a,b},{},{},{}}]

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

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

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

Lorenz equations

demoLorenz1;

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

Special Definition of differential equation:

{a,b}={0,100}; (* Interval *)
system=NDSolveProblem[{{x'[t]==-10 (x[t]-y[t]),y'[t]==-x[t] z[t]+28x[t]-y[t],z'[t]==x[t] y[t]-8/3z[t]},{x[0]==0,z[0]==0,y[0]==1},{x,y,z},{t,a,b},{},{},{}}]

Reference solution by built-in high precision adaptive method:

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

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol],{t,a,b},PlotPoints->1000,ColorFunction->(ColorData["Rainbow"][#4]&)]

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

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

Solution function with adaptive Heun-Euler 2(1) and low accuracy and precision:

solHE=NDSolve[system,StartingStepSize->Automatic,AccuracyGoal->1,PrecisionGoal->1,MaxSteps->Infinity,Method->methodHE,NormFunction->Infinity]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solHE],{t,a,b},PlotPoints->1000,ColorFunction->(ColorData["Rainbow"][#4]&)]

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

coordsHE = First[InterpolatingFunctionCoordinates[First[y/.solHE]]];
Length[coordsHE]

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

Plot the step-sizes as a function:

StepDataPlot[solHE]

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

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 Dormand-Prince 5(4):

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

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

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solDP54],{t,a,b},PlotPoints->1000,ColorFunction->(ColorData["Rainbow"][#4]&)]

coordsDP54 = First[InterpolatingFunctionCoordinates[First[y/.solDP54]]];
Length[coordsDP54]

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

StepDataPlot[solDP54]

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

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

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

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.solBS54],{t,a,b},PlotPoints->1000,ColorFunction->(ColorData["Rainbow"][#4]&)]

coordsBS54 = First[InterpolatingFunctionCoordinates[First[y/.solBS54]]];
Length[coordsBS54]

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

StepDataPlot[solBS54]

Comparison of methods (utility) -> Version 7

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

{a,b}={0,100}; (* Interval *)
system1=NDSolveProblem[{{x'[t]==-10 (x[t]-y[t]),y'[t]==-x[t] z[t]+28x[t]-y[t],z'[t]==x[t] y[t]-8/3z[t]},{x[0]==0,z[0]==0,y[0]==1},{x[t],y[t],z[t]},{t,a,b},{},{},{}}]

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

data=CompareMethods[system1,refsol,{methodDP54},AccuracyGoal->16,PrecisionGoal->16,MaxSteps->Infinity,NormFunction->Infinity];

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

data=CompareMethods[system1,refsol,{methodBS54},AccuracyGoal->16,PrecisionGoal->16,MaxSteps->Infinity,NormFunction->Infinity];

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

Planar 3-body system (only demo)

demoPlanar3Body1

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

demovanderPol1:=Manipulate[
ParametricPlot[
Evaluate[{z[t],v[t]}/.
Quiet@NDSolve[
{z'[t]==v[t],
v'[t]==μ(1-((z[t])^2))v[t]-z[t],
z[0]==zv0[[1]],
v[0]==zv0[[2]]},
{z[t],v[t]},
{t,0,tt}]],
{t,0,tt},
ImageSize->Automatic,PlotRange->4,AxesLabel->{TraditionalForm[z[t]],TraditionalForm[v[t]]}],
{{tt,10,"time t"},.1,100,Appearance->"Labeled"},
{{μ,0.75,"parameter μ"},0,2,0.1,Appearance->"Labeled"},
{{zv0,{0,0}},{-4,-4},{4,4},Locator}]

threeBodyTrajectories[{{x10_,y10_},{x20_,y20_},{x30_,y30_}},
{{vx10_,vy10_},{vx20_,vy20_},{vx30_,vy30_}},{m1_, m2_, m3_},T_, plotType:("x"|"v"),
plotOptions___] := 
Module[{nds,Tmax,prolog,funsToPlot},
nds=NDSolve[
{x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10,y1[0]==y10,x2[0]==x20,y2[0]==y20,x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,vx2[0]==vx20,vy2[0]==vy20,vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3,vx1,vx2,vx3,vy1,vy2,vy3},{t,0,T}];
If[Head[nds]=!=NDSolve,
Tmax = nds[[1, 1, 2, 1, 1, 2]];funsToPlot=If[plotType==="x",{{x1[t],y1[t]},{x2[t],y2[t]},{x3[t],y3[t]}},
{{vx1[t],vy1[t]},{vx2[t],vy2[t]},{vx3[t],vy3[t]}}]/. nds[[1]];
prolog = {PointSize[0.01],Transpose[{{RGBColor[1,.2,0],RGBColor[.5,.8,0],RGBColor[.2,.6,1]},Point /@ (funsToPlot /. t -> 0)}]};
ParametricPlot[Evaluate[funsToPlot],{t,0,Tmax},
PlotStyle->{RGBColor[1,.2,0],RGBColor[.5,.8,0],RGBColor[.2,.6,1]},Frame->False,PlotRange->All,AspectRatio -> 1,
Prolog -> prolog, Frame -> True, MaxRecursion->ControlActive[3,9],Axes -> False, FrameTicks -> False,plotOptions],
Text["Degenerate initial conditions."]]
]//Quiet;

demoPlanar3Body1:=Manipulate[
threeBodyTrajectories[{xy10,xy20,xy30},
{vxy10,vxy20,vxy30},{10^m1Exp, 10^m2Exp,10^m3Exp}, T,xv,ImageSize->{350,350}] ,
{{xv,"x","position/velocity"},{"x" -> "position","v" -> "velocity"}},
{{T,3,"time"},0.001, 10},
{{xy10,{ 0.97000436,-0.24308753},"(Subscript[(OverVector[x]), 1])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{xy20,{-0.97000436,0.24308753},"(Subscript[(OverVector[x]), 2])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{xy30,{0,0},"(Subscript[(OverVector[x]), 3])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{vxy10,{0.93240737/2,0.86473146/2},"(Subscript[(OverVector[v]), 1])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{vxy20,{0.93240737/2,0.86473146/2},"(Subscript[(OverVector[v]), 2])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{vxy30,{-0.93240737,-0.86473146},"(Subscript[(OverVector[v]), 3])[0]"},{-2,-2},{2,2},ImageSize->Small},
{{m1Exp,0,"(Subscript[m, 1])"},-3,3},{{m2Exp,0,"(Subscript[m, 2])"},-3,3},{{m3Exp,0,"(Subscript[m, 3])"},-3,3},
ControlPlacement->{Top,Top,Left,Left,Left,Right,Right,Right,Bottom,Bottom,Bottom},SaveDefinitions->True,AutorunSequencing->{3,5,7}]

demoLorenz1:=Manipulate[Clear[x,y,z,t];
Module[{soln},
With[{β=bb,σ=σσ,ρ=rr,x0=xx0,y0=yy0,z0=zz0},soln=Quiet@NDSolve[SetPrecision[{x'[t]==σ (-x[t]+y[t]),y'[t]==ρ x[t]-y[t]-x[t] z[t],z'[t]==x[t] y[t]-β z[t],x[0]==x0,y[0]==y0,z[0]==z0},Infinity],{x[t],y[t],z[t]},{t,0,tt},PrecisionGoal->ControlActive[4,8],WorkingPrecision->ControlActive[MachinePrecision,20]];ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln],{t,0,soln[[1,1,2,0,1,1,2]]},PlotPoints->ControlActive[100,200],MaxRecursion->ControlActive[4,6],Axes->False,Boxed->False,ImageSize->{400,400},PlotRange->All,ColorFunction->(ColorData["Rainbow"][#4]&)]]],{{tt,10,"time"},1,200,.01,ImageSize->Small,Appearance->"Labeled"},
Delimiter,
Style["parameters",Bold],
{{bb,8/3,"β"},0,10,.01,ImageSize->Small,Appearance->"Labeled"},{{σσ,10,"σ"},0,100,.01,ImageSize->Small,Appearance->"Labeled"},{{rr,28,"ρ"},0,100,.01,ImageSize->Small,Appearance->"Labeled"},Delimiter,Style["initial conditions",Bold],
{{xx0,1,"(Subscript[x, 0])"},0,200,.01,ImageSize->Small,Appearance->"Labeled"},{{yy0,5,"(Subscript[y, 0])"},0,200,.01,ImageSize->Small,Appearance->"Labeled"},{{zz0,10,"(Subscript[z, 0])"},0,200,.01,ImageSize->Small,Appearance->"Labeled"},ControlPlacement->Left,AutorunSequencing->{3,4,5},TrackedSymbols->Manipulate]