Még mindig differenciálegyenletek

Előző házi feladatok

Ezekhez kapcsolódó új feladatok jönnek, illetve az órán is beszéltem a megoldásukról. Csak az utolsóhoz másolom be Holló Gábor megoldásait egy kis (előremutató:)) kiegészítéssel. Az új házi feladatok a szöveg között találhatók.

Tanulmányozzuk a logisztikus egyenletet azzal a feltétellel, hogy szüretelünk is, vagy másképp lehalásszuk a halak egy részét, mégpedig állandó sebességgel, azaz tekintsük ezt az egyenletet:

x'[t]==x[t](1-x[t])-h

ahol h pozitív állandó. Ennek értékétől hogyan függ a rendszer viselkedése?

Mi a helyzet periodikus lehalászás esetén, azaz itt:

x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t])

Rajzoljuk meg a fázisegyenest az alábbi egyenlethez:

x'[t]==x[t](1-x[t])(x[t]-4)

Rajzoljuk meg a fázisegyenest az alábbi egyenlethez:

α'[t]==A Exp[-(E/(R ((Subscript[T, 0])+β t)))](((1-α[t]))^n)

Itt minden együttható pozitív, kivéve n -et, ez lehet 0,1,2,3.

Rajzoljuk meg a fázissíkot kétváltozós lineáris rendszerek összes kvalitatíve különböző esetére. Tehát az egyenlet:

{x'[t]==a x[t]+b y[t],y'[t]==c x[t]+d y[t]}

és az A={{a,b},{c,d}} mátrix sajátértékei valósak, különbözők, egybeesők, komplexek stb.

A={{1,0},{0,1}};
MatrixForm[A]
Eigenvalues[A]
StreamPlot[{A[[1,1]] x+A[[1,2]]y,A[[2,1]]x+A[[2,2]] y},{x,-1,1},{y,-1,1}]

Először is: a mátrixvektor-szorzás használható.

A.{x,y}

Persze, hogy nem változik, hiszen most A az egységmátrix!

Másodszor: Mindig ugyanazt csináljuk; szedjük össze amink van, és írjuk egy programot, ami automatizálja a folyamatot.

phasePortrait[mat_]:=Module[{},Column[{MatrixForm[mat],StreamPlot[mat.{x,y},{x,-1,1},{y,-1,1},PlotLabel->"Eigenvaules: "<>ToString[Eigenvalues[mat]]<>"\n"]}]]

Próbáljuk ki az előző feladaton. (A felirat után beraktam egy új sort, hogy ne tapadjon az ábárhoz.)

phasePortrait[A]

phasePortrait[{{-1,0},{0,-1}}]

Megtáltosodunk. A többit egy csapásra intézzük el.

Map[phasePortrait,{{{1,0},{0,-1}},{{1,-1},{1,1}},{{1,1},{0,1}},{{1,-1},{1,-1}},{{0,-1},{1,0}},{{0,0},{0,-1}}}]

1. HF.: Egészítsük ki a programot azzal, hogy írja ki a stacionárius pont típusát is: nyereg, fókusz stb.

A Cell|DeleteAllOputput utasítást már menet közben is kiadom, hogy az időszakos mentés ne menjen oly lassan.

Fázisegyenes, bifurkáció, Poincaré-leképezés szüretre (lehalászás), fázissík

Fázisegyenes

Megcsináljuk lépésenként, majd összerakjuk.

pnds=ParametricNDSolve[{x'[t]== x[t](1-x[t]),x[0]==x0},x,{t,-2,2},{x0}]

pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]]

rhs=ParametricPlot[{-x(1-x),x},{x,-.2,1.2},PlotStyle->Directive[Red,Thick]]

Figyeljük meg, hogyan készült az elforgatott görbe.

plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}]

sp=StreamPlot[{1,x(1-x)},{t,-1,1},{x,-2,2}]

ph=Graphics[{Thick,Line[{{0,-2},{0,2}}],PointSize[Large],Point[{0,0}],Point[{0,1}]}]

2. HF.: Rajzoljunk nyilat az egyenesre, amely mutatja a trajektóriák haladásának irányát. Szorgalmi: a megoldásokra is!

GraphicsRow[{sp,plr,ph}]

Show[sp,plr,ph]

Szedjük össze. Három alapvető konstrukció van: Module, Block, With

Idemásolom a lépéseket.

pnds=ParametricNDSolve[{x'[t]== x[t](1-x[t]),x[0]==x0},x,{t,-2,2},{x0}]

pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]]

rhs=ParametricPlot[{-x(1-x),x},{x,-.2,1.2},PlotStyle->Directive[Red,Thick]]

plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}]

sp=StreamPlot[{1,x(1-x)},{t,-1,1},{x,-2,2}]

ph=Graphics[{Thick,Line[{{0,-2},{0,2}}],PointSize[Large],Point[{0,0}],Point[{0,1}]}]

GraphicsRow[{sp,plr,ph}]

Show[sp,plr,ph]

Értelmetlenül bepakolom egy Module-ba, csak azért, hogy lássuk annak a szintaxisát.

phaseLine1[rhsfunval_]:=Module[{},pnds=ParametricNDSolve[{x'[t]== x[t](1-x[t]),x[0]==x0},x,{t,-2,2},{x0}];
pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]];
rhs=ParametricPlot[{-x(1-x),x},{x,-.2,1.2},PlotStyle->Directive[Red,Thick]];
plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}];
sp=StreamPlot[{1,x(1-x)},{t,-1,1},{x,-2,2}];
ph=Graphics[{Thick,Line[{{0,-2},{0,2}}],PointSize[Large],Point[{0,0}],Point[{0,1}]}];
GraphicsColumn[{GraphicsRow[{sp,plr,ph}],
Show[sp,plr,ph,PlotRange->{{-0.5,1.0},{-0.5,1.5}}]}]]

phaseLine1[x(1-x)]

Mindegy, mit adunk be neki, ugyanazt rajzolja.

phaseLine1[x(1-(x^2))]

Most három dolgot módosítunk: bevezetjüük a valódi változót rhsfunval néven, lokális változókat írunk a Module megfelelő listájába, és végül az x szám helyébe az x(t) függvényértéket írjuk.

x(1-x)/.{x->x[t]}

phaseLine2[rhsfunval_]:=Module[{pnds,pl2,rhs,plr,sp,ph},pnds=ParametricNDSolve[{x'[t]== (rhsfunval/.{x->x[t]}),x[0]==x0},x,{t,-2,2},{x0}];
pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]];
rhs=ParametricPlot[{-rhsfunval,x},{x,-.2,1.2},PlotStyle->Directive[Red,Thick]];
plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}];
sp=StreamPlot[{1,rhsfunval},{t,-1,1},{x,-2,2}];
ph=Graphics[{Thick,Line[{{0,-2},{0,2}}],PointSize[Large],Point[{0,0}],Point[{0,1}]}];
GraphicsColumn[{GraphicsRow[{sp,plr,ph}],
Show[sp,plr,ph,PlotRange->{{-0.5,1.0},{-0.5,1.5}}]}]]

phaseLine2[x(1-x)]

phaseLine2[x(1-(x^2))]

Jó lenne a stacionárius pontokat automatikusan megkapni.

sol=Solve[x(1-x)==0,Reals]

x/.sol

... és lerajzolni.

Point[{0,#}]&/@(x/.sol)

Prepend[Point[{0,#}]&/@(x/.sol),PointSize[Large]]

Graphics[{Thick,Line[{{0,-2},{0,2}}],PointSize[Large],Point[{0,0}],Point[{0,1}]}]

Join[{Thick,Line[{{0,-2},{0,2}}]},Prepend[Point[{0,#}]&/@(x/.sol),PointSize[Large]]]

Graphics[Join[{Thick,Line[{{0,-2},{0,2}}]},Prepend[Point[{0,#}]&/@(x/.sol),PointSize[Large]]]]

Automatikusan határozzuk meg az egyenesszakasz két végpontját is.

phaseLine3[rhsfunval_]:=Module[{pnds,pl2,rhs,plr,sp,ph},pnds=ParametricNDSolve[{x'[t]== (rhsfunval/.{x->x[t]}),x[0]==x0},x,{t,-2,2},{x0}];
pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]];
rhs=ParametricPlot[{-rhsfunval,x},{x,-2,2},PlotStyle->Directive[Red,Thick]];
plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}];
sp=StreamPlot[{1,rhsfunval},{t,-2,2},{x,-2,2}];
sol=Solve[rhsfunval==0,x,Reals];from=Min[x/.sol];to=Max[x/.sol];
ph=Graphics[Join[{Thick,Line[{{0,from-1},{0,to+1}}]},Prepend[Point[{0,#}]&/@(x/.sol),PointSize[Large]]]];
GraphicsColumn[{GraphicsRow[{sp,plr,ph}],
Show[sp,plr,ph,PlotRange->{{-2,2},{-2,1.5}}]}]]

phaseLine3[x(1-x)]

phaseLine3[x(1-(x^2))]

phaseLine3[x]

phaseLine3[-x]

phaseLine3[0 x]

Itt még többször szerepel from és to, tehát több dolog illeszkedik automatikusan a megoldandó feladathoz.

phaseLine4[rhsfunval_]:=Module[{pnds,pl2,rhs,plr,sp,ph,from,to,sol},pnds=ParametricNDSolve[{x'[t]== (rhsfunval/.{x->x[t]}),x[0]==x0},x,{t,-2,2},{x0}];
pl2=Plot[Evaluate[Table[x[x0][t] /. pnds, {x0, -.1, 1.1, .2}]], {t, 0,2}, PlotRange -> {-2,2},AspectRatio->GoldenRatio,PlotStyle->Directive[Thickness[0.01]]];
rhs=ParametricPlot[{-rhsfunval,x},{x,-2,2},PlotStyle->Directive[Red,Thick]];
plr=Show[rhs,pl2,PlotRange->{-0.5,1.2}];
sol=Solve[rhsfunval==0,x,Reals];
from=Min[x/.sol];to=Max[x/.sol];
sp=StreamPlot[{1,rhsfunval},{t,-2,2},{x,from-1,to+1}];
ph=Graphics[Join[{Thick,Line[{{0,from-1},{0,to+1}}]},Prepend[Point[{0,#}]&/@(x/.sol),PointSize[Large]]]];
Column[{GraphicsRow[{sp,plr,ph}],Show[sp,plr,ph,PlotRange->{{from-1,to+1},{from-1,to+1}}]}]]

phaseLine4[x(1-x)]

phaseLine4[x(1-(x^2))]

phaseLine4[x(1-x)(x+2)]

3. HF.: Alkalmazzuk minél többre, hogy kiderüljenek a hibái. Fejlesszük, ami nem jó benne. (Kétváltozós rendszerek?)

Példul az alábbiakra az a paraméter rögzített értéke mellett.

(Subscript[f, 1])[a_,x_]:=(x^2)-a x

(Subscript[f, 2])[a_,x_]:=(x^3)-a x

f[a_,x_]:=(x^3)- x+a

f[a_,x_]:=a x+Sin[x] (a∈[-1,1])

Bifurkáció (egyszerű példák)

Hogyan függ a hosszútávú viselkedés a paraméterektől?

x'[t]==k x[t]

Ha k=0, minőségi változás következik be: bifurkáció.

rhsfun=(x^2)(a-x);

s=Solve[rhsfun==0,x]

Evaluate[x/.s]

Plot[Evaluate[x/.s],{a,-2,2},PlotStyle->{Thickness[0.01],Red,Thickness[0.01],Blue,Thickness[0.01],Brown},AxesLabel->{a,(SubStar[x])},LabelStyle->Directive[{Large,Bold}] ]

bif[rhsfun_]:=Module[{s},s=Solve[rhsfun==0,x];Plot[Evaluate[x/.s],{a,-2,2},PlotStyle->{Thickness[0.01],Red,Thickness[0.01],Blue,Thickness[0.01],Brown},AxesLabel->{a,(SubStar[x])},LabelStyle->Directive[{Large,Bold}] ]]

bif[(x^2)(a-x)]

bif[ x(1-x)-a]

bif[(x^2)-a x]

(Subscript[f, 2])[a_,x_]:=(x^3)-a x

f[a_,x_]:=(x^3)- x+a

f[a_,x_]:=a x+Sin[x] (a∈[-1,1])

4. HF.: Tekintsük ezt: f(a,b,x):=a x-(x^3)-b. Bifurkációs görbék: f(1,b,x),f(0,b,x),f(-1,b,x). Többi a esetén?

A többi szorgalmi:

Rögzítsük b értékét. Hogyan függ a-tól?

Derítsük fel az egész (a,b) síkot. Milyen bifurkációt tapasztalunk a határoló görbéknél?

M a helyzet a (0,0) pontban, ha a  (és esetleg b is) változik?

Külön a vegyészeknek:

Ugyanez a Schlögl-reakcióra: f(a,b,x):=1-a x+b (x^2)-(x^3) (Hogyan jött ki a dimenziómentes alak? Mi az új paraméterek jelentése?)

Hyperlink["Brusselator","http://demonstrations.wolfram.com/HopfBifurcationInTheBrusselator/"]

Az alábbihoz kell a ReactionKinetics csomag:

Quit[];

SetDirectory[NotebookDirectory[]];

Needs["ReactionKinetics1`"];//AbsoluteTiming

sch=GetProblem["Schlögl"]

DetailedBalanced[sch,{(Subscript[k, 1]),(Subscript[k, -1]),(Subscript[k, 2]),(Subscript[k, -2])}]

Poincaré-leképezés

Numerikusan; szimbolikusan szinte sohasem megy.

NDSolve

nd[p_]:=NDSolve[{x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t]),x[0]==p},x,{t,0,1}][[1,1,2]];

nd[1]

nd[1][1]

poincare[p_]:=nd[p][1]

poincare[0.44]

Plot[{p,poincare[p]},{p,0.32,1.32},AxesOrigin->{0,0}]

FindRoot[p==poincare[p],{p,0.8}]

ParametricNDSolve

poi=ParametricNDSolve[{x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t]),x[0]==p},x,{t,0,1},{p}]

(x[1]/.poi)[1]

Plot[Evaluate[{p,(x[p]/.poi)[1]}],{p,0.32,1.32},GridLines->Automatic]

WhenEvent

ClearAll[nd2];nd2=ParametricNDSolveValue[{x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t]),x[0]==p,WhenEvent[t==1,{Sow[x[t]],"StopIntegrate"}]},x,{t,0,1},{p}];

Reap[nd2[0.45]]

ClearAll[nd3];nd3=ParametricNDSolveValue[{x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t]),x[0]==p,WhenEvent[t==1,{Sow[x[t]],"StopIntegrate"}]},x,{t,0,1},{p},Method->{"ParametricCaching"->None}];

Reap[nd3[0.45]]

poinc[p_]:=First@First@Last@Reap[nd3[p]]

poinc[0.45]

A WhenEvent - be bonyolultabbat is írhatok.

Plot[{p,poinc[p]},{p,0.32,1.32}]

Még másik

poi2[p_]=x[p]/.poi

poi3[p_]=poi2[p][1]

Plot[Evaluate[poi3[p]],{p,0.32,1.32}]

Fixpontok

Plot[{Evaluate[poi3[p]],p},{p,0.32,1.32}]

NestList[poi3,1,4]

NestList[poi3,.41,4]

NestList[poi3,.31,4]

FixedPointList[poi3,1]

Milyen (SuperStar[h]) esetén létezik egyetlen fixpontja a Poincaré leképezésnek?

pnds=ParametricNDSolve[{x'[t]==5x[t](1-x[t])-h(1+Sin[2π t]),x[0]==p},x,{t,0,1},{p,h}]

y17=x[1,0.7]/.pnds

y17[1]

poincare[p_,h_]:=(x[p,h]/.pnds)[1]

poincare[1,0.7]

Manipulate[Plot[{p,poincare[p,h]},{p,0.32,1.32},AxesOrigin->{0,0}],{h,1.14,1.20,0.01}]

Manipulate[Plot[{p,poincare[p,h]},{p,0.32,1.32},AxesOrigin->{0,0}],{h,1.,1.30,0.01}]

5. HF.: Keressenek ilyesmit a demonstrations.wolfram.com honlapon, és röviden ismertessenek egy-egyet közülük:

demo : poincaré map -disk

demo : bifurcation