VisualDSolve

First of all: many many thanks to Antonin for reviving this project: working from a floppy disk in the library in Prague!

http://store.wolfram.com/view/book/D0706.str

Equilibrium Points

This is done by looking at contours....

Clear[f,g];
f[x_,y_]:= x - (y^2) Cos[y];
g[x_,y_]:=-y +x Sin[x];
fcnVec[{x_,y_}]:= {f[x,y], g[x,y]};

Short[fZero=Cases[Normal[ContourPlot[f[x,y]==0,{x,-10,10},{y,-10,10}]],Line[z_]:>z,∞],15]

Show[Graphics[Point/@fZero],ContourPlot[g[x,y]==0,{x,-10,10},{y,-10,10},ContourStyle->{Thick,Red}],Frame->True, FrameTicksStyle->16]

(signs=Sign[Apply[g,fZero[[3]],{1}]])//Short

(changes=Rest[signs RotateRight[signs]])//Short

Next find the positions of the sign changes.

posns=Flatten[Position[changes,-1]]

Find the corresponding points.

seeds=fZero[[3,Flatten[posns]]]

Refine the points using Newton's method via FindRoot.

roots=Table[v/.FindRoot[fcnVec[v],{v,s}],{s,seeds}]

Check visually that all is correct.

Show[ContourPlot[{f[x,y]==0,g[x,y]==0}, {x,-10,10}, {y,-10,10}, ContourStyle->{{Thick,Red},{Thick,Black}}],Graphics[{Thickness[0.005],Blue,Table[Circle[r, 0.3], {r, roots}]}], Frame->True, FrameTicksStyle->16]

Clear[f,g,x,y];
f[x_,y_]:=-Cos[y]+2y Cos[(y^2)]Cos[2x];
g[x_,y_]:=-Sin[x]+2Sin[(y^2)]Sin[2x];
roots=FindRoots2D[{f[x,y], g[x,y]},  {x, -3.5, 4}, {y, -1.8, 4.2}];
Length[roots]

Show[ContourPlot[{f[x,y]==0,g[x,y]==0},{x,-3.5,4},{y,-1.8,4.2},ContourStyle->{{Thick,Red},{Thick,Blue}}],Graphics[{EdgeForm[Thickness[0.003]],Yellow,(Disk[#1,0.07]&)/@roots}],AspectRatio->Automatic, FrameTicksStyle->16]

A Comprehensive FindRoots2D Function

FindRoots2D::usage=
"FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to {f=0, g=0} in the given rectangle.";

Options[FindRoots2D]={PlotPoints->Automatic, MaxRecursion->Automatic};

FindRoots2D[funcs_,{x_,a_,b_},{y_,c_,d_},opts___]:=Module[{fZero,seeds,signs,fy},
fy=Compile[{x,y},Evaluate[funcs[[2]]]];
fZero=Cases[Normal[ContourPlot[funcs[[1]]==0,{x,a-((b-a)/97),b+((b-a)/103)},{y,c-((d-c)/98),d+((d-c)/102)},Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]], Line[z_] :> z,∞];
seeds=Flatten[(signs=
Sign[Apply[fy,#,{1}]];#1[[1+Flatten[Position[Rest[signs * RotateRight[signs]],-1]]]])&/@fZero,1];
If[seeds=={},{},Select[Union[
({x,y}/.FindRoot[{funcs[[1]],funcs[[2]]},{x,#[[1]]},{y,#[[2]]},Evaluate[FilterRules[{opts}, Options[FindRoot]]]]&)/@seeds, SameTest -> (Norm[#1-#2]<10^-6&)],a≤#[[1]]≤b&&c≤#[[2]]≤d&]]]

It would be nice if Mathematica had a FindAllRoots function that was a little more general than they have now. The contour method is certainly robust in 2-dim, and 1-dim is not really a problem.

Applications

Flipping a book

Application: Saving the MIR space station, by Michael Foale, astronaut and physicist.

Throwing a discus

Using the optimal throwing angles (release, attack) with a headwind and a tailwind. Headwind is better.

Duffing equation: Strange attractor

Book discusses how one can checking correctness of numerical solutions.

DuffingEquation[{1,1,0.15,0.3,1},x[t]]

Following illustrates the Poincaré section function in the package, but uses it in 3-D, breaking up the solutions mod 2π.

Clear[t,x,y];Manipulate[sol=x[t]/.First[NDSolve[{DuffingEquation[{1,1,0.15,0.3,1},x[t]],x[0]==a,x'[0]==b},x[t],{t,0,tmaxForDE},MaxSteps->∞]];
Show[
ParametricPlot3D[Evaluate[{Mod[t,2π],sol,(∂_t)sol}],{t,0,tmaxForDE},PlotStyle->{Purple,Thickness[0.004]}, Exclusions-> Table[2 π i, {i, 0, 100}]],
Graphics3D[{{PointSize[0.002],Red,Point[ps]},
{Blue,PointSize[0.02],Point[{0, a,b}],Point[Table[Evaluate[{ 2π,sol,(∂_t)sol}/.t->2 π i],{i,(tmaxForDE/(2 π))}]]},
{Gray,Line[Table[{{0,j,-1.2}, {0,j,1.2}}, {j,-1.5,1.5,.5}]],
Line[Table[{{0,-1.4,j}, {0,1.7,j}}, {j,-1,1,.25}]]}}],ImageSize->400,PlotRange->{{0.001, 6.28},{-1.4,1.7},{-0.7,1.2}},BoxRatios->{6.4, 3,2.6},
TicksStyle->14,Ticks->{Range[0,2π,π], Range[-1,1], Range[-.5,1,.5]},LabelStyle->16,AxesLabel->{"t (mod 2 π)", "x", "x'"},ViewPoint->3{-1.6, -2, .4}, AxesEdge->{{-1,-1},{-1,-1},{-1,1}}],
{{a,0, "x(0)"}, -1,1},{{b,0, "x'(0)"},-1,1},{{tmaxForDE,15,"max time for solution"},0.001,150},TrackedSymbols->{a,b,tmaxForDE},TrackedSymbols->Manipulate, Initialization :> (Needs["VisualDSolve`VisualDSolve`"];ps =Prepend[#, 2π]& /@  PoincareSection[DuffingEquation[{1,1,0.15,0.3,1},x[t]],{t,0,2 π 3000},{x,-1.5,1.5},{x',-1.3,1.3},{-1,1},TimeInterval->2 π,SavePoints->True,MaxSteps->∞][[2]])]

Some new demos

MakePhasePlotDemo[]

Chapter 18: Miscellaneous New Topics

Square wheels

Square wheel demo: general: LOTS of computation here!

Getting the wheel for a road is tricky, in that the height of the road determines whether the wheel will close up. Further, the closing up occurs if a certain integral is an integer (2, 3, 4).

Change height in following demo to see how the "closed-wheel condition" can fail.
Then get closed wheels automatically by using FindRoot on the integral.
Interpolation is used for the road.
So this demo combines several ideas in numerical computation.

Module[{c, tmaxForWheelPlot, xval, p,t,nper}, 
Manipulate[Module[{road, θSol,tmin, tt,p2,f,func,θS,transform,a,cc, newc, bound},
a=ArcSinh[1.];xval = Min[xval, 20 tmaxForWheelPlot];
road[yfunc_] :=road[yfunc]=Plot[yfunc,{t,-5,20},PlotStyle->{Red,Thick},PlotPoints->50, Filling->-2, FillingStyle->LightOrange];
θSol[fcn_, {tmin_, tmax_}]:=θSol[fcn]=Module[{θ},θ[t]/.First[NDSolve[{θ[0]==-(π/2),θ'[t]==-(1/fcn)},θ[t],{t,tmin,Min[8,20 tmax]}]]];
tmin=-1;
bound[num_, l_, h_]:= Max[Min[num, h], l];
p ={{-a, bound[p[[1,2]], -1.75, -.9]},{-.5, bound[p[[2,2]], -1.4,-.9]},{0, bound[p[[3,2]], -1., -.8]},{0.5, bound[p[[4,2]], -1.4, -.8]}};
p2= Append[p,{a, p[[1,2]]}];
f=Min[-0.2,Interpolation[p2,PeriodicInterpolation->True][t]];
If[IntegerQ[nper],
pmax = Abs[SortBy[N[p2], Last][[-1,-1]]];
tmaxForWheelPlot=2. a nper-1;
objective[cc_?NumericQ] := -1/Quiet@NIntegrate[(1/(2 π (f+cc))),{t,-a,a}];
mPrint[{c,pmax, {cc,c-0.25,Min[pmax-c,c+0.25]}}];
c=cc/.FindRoot[objective[cc]==nper, {cc,c-0.25,Min[Min[pmax-c,c+0.25], -0.05]}]; 
mPrint[c];
objective[cc_?NumericQ] :=(f=Min[-0.2,Interpolation[{0,cc}+# & /@ p2,PeriodicInterpolation->True][t]];
func[tval_] := f /. t-> tval; 
θS=θSol[func[t], {tmin,tmaxForWheelPlot }];
((Norm[(f {Cos[θS],Sin[θS]}/. t->-1)-(f {Cos[θS],Sin[θS]}/. t->tmaxForWheelPlot)])^2));
c=cc/.FindMinimum[objective[cc],{cc,c-0.1,Min[-0.05, c+0.1]}][[2]]; nper=False];
mPrint[c];p1={0,c}+# & /@ p2;
f=Min[-0.2,Interpolation[p1,PeriodicInterpolation->True][t]];
func[tval_] := f /. t-> tval; 
θS=θSol[func[t], {tmin, tmaxForWheelPlot }];
transform[pt_, tval_]:= TranslationTransform[{tval,0}][RotationTransform[-(π/2)-θS/.t->tval][pt]];
Column[{
Style[StringForm["Integral for wheel closure = (Subscript[∫, -a]^a)((-1)/(2 π road(x))) \[DifferentialD]x = ``, where a = arcsinh 1.",NumberForm[intval=-1/Quiet[NIntegrate[(1/(2 π func[tt])),{tt,-a,a}]],3]], FontFamily->"Helvetica", 12],
Show[road[func[t]],
wheel=ParametricPlot[Evaluate[transform[-f {Cos[θS],Sin[θS]}/. t->tt,xval]],{tt,tmin,Min[8,tmaxForWheelPlot]},PlotStyle->Thickness[0.006], PlotPoints->50];
Graphics[If[(Abs[intval-2]<0.02|| (Abs[intval-3]<0.01)||(Abs[intval-4]<0.01))&&(1.99< (tmaxForWheelPlot+1)/(2 a)<2.01 ||2.99<(tmaxForWheelPlot+1)/(2 a)<3.01 ||3.99<(tmaxForWheelPlot+1)/(2 a)<4.01 ), Graphics[{EdgeForm[], FaceForm[Lighter[Yellow]],Cases[Normal[wheel], Line[z_] :> Polygon[Drop[z,-10]], ∞]}],{}]],
wheel,
Graphics[{PointSize[0.018], Point[cen={xval,0}],PointSize[0.015],
tangentPt=-func[xval] {Cos[θS],Sin[θS]}/. t-> xval;tangentPtOrig={0,func[0] }; 
Point[transform[tangentPt,xval]],
Point[{xval,0}+RotationTransform[-(π/2)-θS/.t->xval][tangentPtOrig]],
Line[{cen,transform[ tangentPtOrig,xval]}],If[locQ, Locator/@ p, {}],
{Gray, Thickness[0.003],Line[{{{-a,-2},{-a,2}},{{a,-2},{a,2}},{{-2,0},{5.5,0}}}]}}],
PlotRange->{{-2, 5.5}, {-2,2}},PlotLabel->None,
Frame->True,FrameTicks->{{ Range[-2,2],{}},{{-1,0,1,2,3,4,5},{{-a,"-a"},{a,"a"}}}},ImageSize->400,ImagePadding->20,PlotRangeClipping->True,PlotRangePadding->0,AspectRatio->Automatic, Axes->False]}]],
{{c,0.0,"change road height"},-1,0.71},
{{tmaxForWheelPlot,5.3,"upper limit on wheel integration"},0.1, 8}, 
{{xval,0,"roll the wheel"},-1, 5},
{{locQ, True, "show locators"}, {True, False}},
{{nper, False, "try to find closed wheel over periods"}, {2,3,4}}, 
{{p, {{-a,-(Sqrt[2.])},{-0.5,-1.07},{0,-1.01},{0.5,-1.14}}}, Locator, Appearance->None},
 TrackedSymbols->{p, nper,xval, tmaxForWheelPlot,c, locQ}, SaveDefinitions->True, ControlPlacement->Top]]

More general wheels: LOTS of computation here!

Bike Tracks

Module[{xmin, xrange, ymin, yrange, locQ,p, tanQ, tantt,met},Manipulate[  Module[{ranPts,t,backTrack, frontTrack},

backTrack=Which[met==="Bezier",BezierFunction[p, SplineDegree->Length[p]-1][t], met==="Spline",Interpolation[Transpose[{Range[0,1,1/(Length[p]-1)],p}], Method->"Spline"][t]];
frontTrack=backTrack+randir (((∂_t)backTrack)/(Norm[(∂_t)backTrack]));
Show[
ParametricPlot[Evaluate[{frontTrack, backTrack}],{t,0.01,1},PlotStyle->{{Thick,Red},{Thick,Blue}},PlotRange->All],
Graphics[{
If[locQ,{EdgeForm[Thickness[0.0025]],FaceForm[LightYellow], Disk[#,0.13]& /@ p, MapIndexed[Text[Style[#2[[1]], 9],#]&,p]}, {}],
If[tanQ, {Thickness[0.008],PointSize[0.02],Point[backTrack], Arrow[{backTrack, frontTrack}]}/. t-> If[randir==1,tantt,1-tantt],{}]}], 
Frame->True, FrameTicks->False,ImagePadding->1,Axes->False,
PlotRange->{{xmin, xmin+xrange},{ymin, ymin+yrange}},ImageSize->300]],
{{p, {{0, 1},{0.6, 1.2},{1.2,1.8}, {3, .5}, {4,.4}, {4.5, .42}, {5.5, 2}} }, Locator, LocatorAutoCreate->True, Appearance->None}, 
{{randir,1}, {-1,1}, None},
{{locQ, False, "show locators"}, {True, False}},
{{met,"Bezier", "curve drawing method"}, {"Bezier", "Spline"}},
{{xmin,1.06, "xmin"}, -3,3, ImageSize->Small}, {{xrange,3.9, "range of x"}, 1, 10, ImageSize->Small}, {{ymin,0.28, "ymin"}, -3,3, ImageSize->Small}, {{yrange,1.64, "range of y"}, 1, 10, ImageSize->Small}, 
Row[{Spacer[58],Button["use random points",p=Sort[Table[{i,RandomReal[met/.{"Bezier"->{-0,3.}, "Spline"->{+.6, 1.5}}]}, {i,0, 6, (6/9)}]], ImageSize->Medium, BaseStyle->10],Spacer[10],
Button["use random direction",randir=Sign[RandomReal[{-1,1}]], ImageSize->Medium, BaseStyle->10]}],
 {{tanQ,False,"show the bicycle"}, {True, False}},
{{tantt,0.5,"bicycle movement parameter"},0.01,0.99, ImageSize->Small},TrackedSymbols->Manipulate, SaveDefinitions->True]]

Module[{front, front1, frontDer,backSol,locQ,xmin,xrange, ymin, yrange,t, b,p,met}, Manipulate[ L=1;
front[t_] := Evaluate[Which[met==="Bezier",BezierFunction[p, SplineDegree->Length[p]-1][t], met==="Spline", Interpolation[Transpose[{Range[0,1,1/(Length[p]-1)],p}], Method->"Spline"][t]]];
front1[t_?NumericQ] := front[t];
frontDer[t_?NumericQ] := Evaluate[D[front[t],t]];

backSol=b[t]/.First[NDSolve[
{b'[t]==((frontDer[t].(b[t]-front1[t])(b[t]-front1[t]))/(L^2)),b[0]=={-1,0}+front1[0]},b[t],{t,0,1}]];
Show[
ParametricPlot[front[t], {t,0,1}, PlotStyle->{Red, Thick}],
ParametricPlot[backSol, {t,0,1}, PlotStyle->{Blue, Thick}],
Graphics[{If[locQ,{EdgeForm[Thickness[0.0025]],FaceForm[LightYellow],Disk[#, 0.12]& /@ (p), MapIndexed[Text[Style[#2[[1]], 9],#]&,p]}, {}]}], 
Frame->True,FrameTicks->None,PlotRange->{{xmin, xmin+xrange}, {ymin, ymin+yrange}},ImagePadding->0+4,  
PlotRangePadding->{0,0.3},Axes->False,ImageSize->400],
{{p, {{1,0},{2.05,0},{2.95,0.12},{3.93,0.17},{5.065,0.58},{5.4,1.32},{5.29,2.18},{3.24,2.68},{0.93,3.03}} }, Locator, LocatorAutoCreate->True, Appearance->None}, 
{{locQ, True, "show locators"}, {True, False}},
{{met,"Bezier", "curve drawing method"}, {"Bezier", "Spline"}},
{{xmin,-0.2, "xmin"}, -3,3, ImageSize->Small}, {{xrange,6,"xrange"}, 1, 10, ImageSize->Small},
{{ymin,0,"ymin"}, -3,3, ImageSize->Small}, {{yrange,3,"yrange"}, 1, 10, ImageSize->Small}, 
 TrackedSymbols->Manipulate, SaveDefinitions->True, ControlPlacement->Left]]

Sharks-fish: a predator-prey model with satiation:

Module[{x,y, a,b,c,d,e,f, pr, axes},Manipulate[
pr={{x0,x1},{y0,y1}}={{-0.9, 0.3},{-0.5,0.5}};
{a,b,c,d,f}={1/2, 1,1,1,1};
axes={Red,Thickness[0.015],Line[{{{x0,0},{x1,0}},{{0,y0},{0,y1}}}]};
ept[z_List] :={PointSize[0.12],If[Re[#]<0,Red,Yellow],Point[{Re[#], Im[#]}]}&/@z;
GraphicsRow[{ phasePlaneImage[e],
GraphicsColumn[{Graphics[{Thickness[0.06], Circle[{1/6,0},(Sqrt[7])/6 ],Line[{{x0,0},{x1,0}}],axes,ept[(-2 a+b+2 c-d+d e-4 e f+{1,-1}(Sqrt[(((2 a-b-2 c+d-d e+4 e f))^2)-8 (-2 a c+b c+a d-a d e+4 a e f-2 b e f)]))/4.]},Frame->True, FrameTicks->None, Background->Lighter[Purple], PlotRange->pr,ImageSize->100,Axes->False],
Graphics[{Thickness[0.06],Line[{{x0,0},{x1,0}}],axes,
ept[{-c,((b c-a (c+e f))/(c+e f))}]}, PlotRange->pr,ImageSize->100, Frame->True,FrameTicks->None,Axes->False,Background->Lighter[Orange]]},
   Spacings-> 20]}, Spacings->-100],
{{e,1.2,"satiation parameter e"},1.2,0.125, -0.025, Appearance->"Labeled"}, TrackedSymbols->{e},ControlPlacement->Top, SaveDefinitions->True]]

Lakes of Wada for Forced, Damped Pendulum

An explanation is given for the basins of attraction of a forced damped pendulum: a remarkable collection of shapes:

Init

Get["VisualDSolve`"]

SetVisualDSolveFont["Times"]