Zajęcia IV : Równanie Van der Pola

Zad. 1 (Uproszczenie rozważanego zagadnienia)

Niech x (t) będzie rozwiązaniem równania
 (Overscript[x, ¨])+f'(x)(OverDot[x])+ax=0.
Zauważmy, że  x((t/(Sqrt[a]))) jest rozwiązaniem  równania
 (Overscript[x, ¨])+(1/(Sqrt[a]))f'(x)(OverDot[x])+x=0.
W drugą stronę, jeżeli x (t) jest rozwiązaniem równania
(Overscript[x, ¨])+(1/(Sqrt[a]))f'(x)(OverDot[x])+x=0
to x((Sqrt[a])t) jest rozwiązaniem równania
(Overscript[x, ¨])+f'(x)(OverDot[x])+ax=0.
Wobec tego zamiast badać równanie  (Overscript[x, ¨])+f'(x)(OverDot[x])+ax=0   możemy zająć się równaniem   (Overscript[x, ¨])+(OverTilde[f])'(x)(OverDot[x])+x=0.

Zad. 1I (Dalsze upraszczanie rozważanego zagadnienia)

Niech x (t) będzie rozwiązaniem rówania
 (Overscript[x, ¨])+f'(x)(OverDot[x])+x=0.
Zauważmy, że  x(-t) jest rozwiązaniem  równania
 (Overscript[x, ¨])-f'(x)(OverDot[x])+x=0.
W drugą stronę, jeżeli x (t) jest rozwiązaniem równania
 (Overscript[x, ¨])-f'(x)(OverDot[x])+x=0
to x(-t) jest rozwiązaniem równania
(Overscript[x, ¨])+f'(x)(OverDot[x])+x=0.
Wobec tego zamiast badać równanie  (Overscript[x, ¨])+μ f'(x)(OverDot[x])+x=0
dla μ∈ℝ   możemy zająć się równaniem   (Overscript[x, ¨])+μ f'(x)(OverDot[x])+x=0  
dla μ∈(SubPlus[ℝ])  (μ=0 pomijamy jako trywialne).
Dynamika  zadana przez rozwiązania równania z  parametrem μ  będzie taka sama z dokładnością do odwrócenia czasu jak dynamika zadana przez równanie z parametrem -μ.

Zad. III (Zapisujemy równanie drugiego rzędu jako układ równań rzędu pierwszego)

Zauważamy, że równanie  (Overscript[x, ¨])+f'(x)(OverDot[x])+x=0  jest równoważne układowi równań:
({{{, ({{(OverDot[x])=-f(x)+y,,  }, {(OverDot[y])=-x.,  }})}})

W naszym konkretnym przypadku możemy więc zajmować się następującym układem:
({{{, ({{(OverDot[x])=μ (x-((x^3)/3))+y,,   μ>0,}, {(OverDot[y])=-x. ,  }})}})

Zad. IV (Wstępne badanie punktów stacjonarnych)

f[x_,y_,μ_]:=μ(x-x^3/3)+y
g[x_]:=-x

Szukamy punktów stacjonarnych

Solve[{f[x,y,μ]==0,g[x]==0},{x,y}]

Moża zauważyć, że dla f (x) nieparzystych (0, 0) jest punktem stacjonarnym i jest to jedyny punkt stacjonarny.

Linearyzujemy zadany układ

{D[f[x,y,μ],{{x,y}}],D[g[x],{{x,y}}]}//MatrixForm

Zad. V (Badanie stabilności punktu stacjonarnego (0,0))

wlasne[x_,y_,μ_]:=Evaluate[{D[f[x,y,μ],{{x,y}}],D[g[x],{{x,y}}]}]

Manipulate[ww=N[Eigenvalues[wlasne[0,0,μ]]];
"Wartości własne linearyzacji układu w punkcie stacjonarnym (0,0) dla μ= "<>ToString[μ]<>"\n"<>
"1. "<>ToString[ww[[1]]]<>",\n"<>
"2. "<>ToString[ww[[2]]]<>".\n"<>If[Positive[Re[ww[[1]]]]&&Positive[Re[ww[[2]]]],"Część rzeczywista obu wartości własnych jest dodatnia.",If[Positive[Re[ww[[1]]]*Re[ww[[2]]]],"Część rzeczywista obu wartości własnych jest ujemna.",If[Negative[Re[ww[[1]]]*Re[ww[[2]]]],"Jedna z wartości własnych ma dodatnią częśc rzeczywistą a druga ujemną.",If[Re[ww[[1]]]==Re[ww[[2]]],"Obie wartości mają częśc rzeczywistą równą zero.","Tylko jedna wartość własna ma zerową część rzeczywistą."]]]]
,{μ,0,10}]

Wniosek: na mocy twierdzenia Grobmana - Hartmana dla μ > 0 punkt stacjonarny (0, 0) jest niestabilny.

Zad. VI (Portret fazowy)

a)

VectorPlot[{f[x,y,2],g[x]},{x,-3,3},{y,-3,3},
VectorScale->{Automatic,Automatic,Log[#5+1]&},
StreamPoints->Coarse,
StreamStyle->Orange,
FrameLabel->{"x","y"},
RotateLabel->False,
PlotLabel->"Portret fazowy badanego układu",
LabelStyle->Directive[Bold,Black,24],
Epilog->{Red,PointSize[0.02],Point[{0,0}]},
ImageSize->600]

b)

s1=NDSolve[{x'[t]==f[x[t],y[t],2],y'[t]==g[x[t]],x[0]==2,y[0]==3},{x,y},{t,0,40}];
Plot[Evaluate[{x[t],y[t]}/.s1],{t,0,40},
PlotLegends->LineLegend[{Style["x(t)",18,Blue],Style["y(t)",18,Red]},LegendFunction->Frame],
AxesStyle->Arrowheads[0.03],
AxesLabel->{Style["t",16],ToString[Style["x(t)",16,Blue],StandardForm]<>", "<>ToString[Style["y(t)",16,Red],StandardForm]},
PlotStyle->{Blue,Red},
GridLines->{None,Automatic},
GridLinesStyle->Directive[Gray,Dashed],
PlotLabel->Style["Wykres szeregów czasowych",26],
LabelStyle->Directive[Bold,Black,12],
ImageSize->600]

c 1)

s2=NDSolve[{x'[t]==f[x[t],y[t],2],y'[t]==g[x[t]],x[0]==2,y[0]==3},{x,y},{t,0,40}];
ParametricPlot[Evaluate[{x[t],y[t]}/.s2],{t,0,40},
PlotLegends->Placed[{Style["Badana trajektoria",Italic,18]},Below],
PlotRange->{{-2.5,3.3},{-3,3.5}},
AxesStyle->Arrowheads[0.03],
AxesLabel->{Style["x(t)",16],Style["y(t)",16]},
PlotStyle->{Red,Thick},
GridLines->{Automatic,Automatic},
GridLinesStyle->Directive[Gray,Dashed],
PlotLabel->Style["Wykres badanej trajektorii",26],
LabelStyle->Directive[Bold,Black,12],
Epilog->{{Red,PointSize[0.02],Point[{0,0}]},{Black,PointSize[0.02],Point[{2,3}]},Inset[Style["(x(0),y(0))=(2,3)",Black,Bold,16],{2.5,3.2}]},
ImageSize->600]

c 2)

s3=NDSolve[{x'[t]==f[x[t],y[t],2],y'[t]==g[x[t]],x[0]==2,y[0]==3},{x,y},{t,0,40}];
Show[
ParametricPlot[Evaluate[{x[t],y[t]}/.s3],{t,0,40},
PlotLegends->Placed[{Style["Badana trajektoria",Italic,18]},Below],
PlotRange->{{-2.5,3.3},{-3,3.5}},
AxesStyle->Arrowheads[0.03],
AxesLabel->{Style["x(t)",16],Style["y(t)",16]},
PlotStyle->{Red,Thick},
GridLines->{Automatic,Automatic},
GridLinesStyle->Directive[Gray,Dashed],
PlotLabel->Style["Wykres badanej trajektorii",26],
LabelStyle->Directive[Bold,Black,12],
Epilog->{{Red,PointSize[0.03],Point[{0,0}]},{Black,PointSize[0.02],Point[{2,3}]},Inset[Style["(x(0),y(0))=(2,3)",Black,Bold,16],{2.5,3.2}]},
ImageSize->600],
VectorPlot[{f[x,y,2],g[x]},{x,-2.5,3.3},{y,-3,3.5},
VectorScale->{Automatic,Automatic,Log[#5+1]&},
StreamPoints->Coarse,
StreamStyle->Yellow,
FrameLabel->{"x","y"}
]
]

Zad. VII+VIII (Zbiorcze zestawienie wyników)

Manipulate[pom1=NDSolve[{x'[t]==f[x[t],y[t],μ],y'[t]==g[x[t]],x[0]==start1[[1]],y[0]==start1[[2]]},{x,y},{t,0,t1}];
pom2=NDSolve[{x'[t]==f[x[t],y[t],μ],y'[t]==g[x[t]],x[0]==start2[[1]],y[0]==start2[[2]]},{x,y},{t,0,t2}];
Show[
ParametricPlot[Evaluate[{x[t],y[t]}/.pom1],{t,0,Min[tt1,t1]},
PlotRange->{{-2.5,3.3},{-4.5,4.5}},
AxesStyle->Arrowheads[0.03],
AxesLabel->{Style["x(t)",16],Style["y(t)",16]},
PlotStyle->{Black,Thickness[0.01]},
GridLines->{Automatic,Automatic},
GridLinesStyle->Directive[Gray,Dashed],
PlotLabel->Style["Wykres badanych trajektorii",26],
Epilog->{
{Red,PointSize[0.02],Point[{0,0}]},
{Black,PointSize[0.02],Point[{start1[[1]],start1[[2]]}]},Inset[Style["("<>ToString[If[start1[[1]]==0,0,If[Round[start1[[1]],0.01]==0,"~0",Round[start1[[1]],0.01]]]]<>", "<>ToString[If[start1[[2]]==0,0,If[Round[start1[[2]],0.01]==0,"~0",Round[start1[[2]],0.01]]]]<>")",Bold,16],{start1[[1]],start1[[2]]+0.2}],
{Purple,PointSize[0.02],Point[{start2[[1]],start2[[2]]}]},
Inset[Style["("<>ToString[If[start2[[1]]==0,0,If[Round[start2[[1]],0.01]==0,"~0",Round[start2[[1]],0.01]]]]<>", "<>ToString[If[start2[[2]]==0,0,If[Round[start2[[2]],0.01]==0,"~0",Round[start2[[2]],0.01]]]]<>")",Purple,Bold,16],{start2[[1]],start2[[2]]-0.2}]
},
PlotLegends->Placed[{Style["Badana trajektoria 1",Italic,18]},Below],
ImageSize->500],
ParametricPlot[Evaluate[{x[t],y[t]}/.pom2],{t,0,Min[tt2,t2]},
PlotStyle->{Purple,Thickness[0.01]},
PlotLegends->Placed[{Style["Badana trajektoria 2",Italic,18]},Below]
],
VectorPlot[{f[x,y,μ],g[x]},{x,-2.5,3.3},{y,-4.5,4.5},
VectorStyle->Darker[LightBlue,0.1],
VectorScale->{Tiny,Scaled[.9],None},
StreamPoints->Coarse,
StreamStyle->Darker[LightBlue,0.1],
FrameLabel->{"x","y"}
],
LabelStyle->Directive[Black]
],
{{μ,2,Style["μ",16]},0,5,0.1},
{{start1,{0.2,0.2},Style["((Subscript[x, 1])(0),(Subscript[x, 2])(0))=",16]},},
{{start2,{-1,-2},Style["((Subscript[x, 1])(0),(Subscript[x, 2])(0))=",16]},Locator},
{{t1,2,Style["(Subscript[t, 1])-czas max. trajektorii 1",16]},{1,2,4,8,16,32,64}},
{{tt1,1,Style["Czas istnienia trajektorii 1",16]},1,t1,t1/64},
{{t2,2,Style["(Subscript[t, 2])-czas max. trajektorii 2",16]},{1,2,4,8,16,32,64}},
{{tt2,1,Style["Czas istnienia trajektorii 2",16]},1,t2,t2/64},
ControlPlacement->Right]

Zad. IX (Inny przykład)

Manipulate[pom1=NDSolve[{x'[t]==-(μ x[t]-Sin[x[t]])+y[t],y'[t]==-x[t],x[0]==start1[[1]],y[0]==start1[[2]]},{x,y},{t,0,64}];
Show[
ParametricPlot[Evaluate[{x[t],y[t]}/.pom1],{t,0,t1},
PlotRange->{{-30,30},{-30,30}},
AxesStyle->Arrowheads[0.03],
AxesLabel->{Style["x(t)",16],Style["y(t)",16]},
PlotStyle->Black,
GridLines->{Automatic,Automatic},
GridLinesStyle->Directive[Gray,Dashed],
PlotLabel->Style["Wykres badanych trajektorii",26],
Epilog->{
{Red,PointSize[0.02],Point[{0,0}]},
{Black,PointSize[0.02],Point[{start1[[1]],start1[[2]]}]},Inset[Style["("<>ToString[If[start1[[1]]==0,0,If[Round[start1[[1]],0.01]==0,"~0",Round[start1[[1]],0.01]]]]<>", "<>ToString[If[start1[[2]]==0,0,If[Round[start1[[2]],0.01]==0,"~0",Round[start1[[2]],0.01]]]]<>")",Bold,16],{start1[[1]],start1[[2]]+0.2}],
},
PlotLegends->Placed[{Style["Badana trajektoria",Italic,18]},Below],
ImageSize->500],
VectorPlot[{f[x,y,μ],g[x]},{x,-30,30},{y,-30,30},
VectorStyle->LightBlue,
VectorScale->{Tiny,Scaled[.9],None},
StreamPoints->Coarse,
StreamStyle->LightBlue,
FrameLabel->{"x","y"}
],
LabelStyle->Directive[Black]
],
{{μ,0.1,Style["μ",16]},0.01,0.1},
{{start1,{0.2,0.2},Style["((Subscript[x, 1])(0),(Subscript[x, 2])(0))=",16]},Locator},
{{t1,64,Style["Czas życia trajektorii",16]},{4,8,16,32,64}},
ControlPlacement->Right]

Zuważamy, że dla małego parametru μ > 0 układ zdaje się mieć "dużo" stabilnych coraz większych cykli granicznych. Wraz ze wzrostem parametru μ cykle te zanikają aż pozostaje tylko jeden taki cykl stabilny w otoczeniu (0, 0) (np dla μ = 0.1), jeżeli zwiększymy μ do wartości μ = 1 lub więcej wówczas nie obserwujemy cykli granicznych (punkt stacjonarny przyciąga wszystkie trajektorie).