Differenciálegyenletek kezelése Mathematicával

Bemutató matematikus és fizikus hallgatóknak. Moson Péter és Nagy Katalin előadásához összeszedte: Tóth János. 2009. december

Megoldás és trajektória

Inga

nds=NDSolve[{x'[t]==y[t],y'[t]==-Sin[x[t]],x[0]==0,y[0]==1},{x,y},{t,0,15}][[1]];
ParametricPlot3D[Evaluate[{{0,x[t],y[t]},{t,x[t],-1.5},{t,1.5,y[t]},{t,x[t],y[t]}}/.nds],{t,0,15},PlotStyle->{Directive[Thick,Blue],Directive[Thick,Green],Directive[Thick,Brown],Directive[Thick,Red]},BoxRatios->{4,1,1},AxesLabel->{"t","x(t)","y(t)"},PlotLabel->"Az inga egyenletének megoldása,\na megoldás koordinátafüggvényei  és\n trajektóriája a fázissíkon"]

Első integrál

φ[p_,q_]:=((q^2)/2)-Cos[p];f[p_,q_]:={q,-Sin[p]};

D[φ[p,q],{{p,q}}].f[p,q]

ContourPlot[((y^2)/2) - Cos[x], {x, -5, 8}, 
       {y, -2, 2},AspectRatio->Automatic]

Brüsszelátor

brusszel={x'[t]==1+a ((x[t])^2)y[t]-(5+1)x[t],y'[t]==a(5 x[t]-a ((x[t])^2)y[t]),x[0]==8,y[0]==3}/.{a->1.96};

Framed[TraditionalForm[Column[brusszel]]]

ClearAll[nds];
nds=NDSolve[brusszel,{x,y},{t,0,100}][[1]];
ParametricPlot3D[Evaluate[{{90,x[t],y[t]},{t,x[t],0.5},{t,1.6,y[t]},{t,x[t],y[t]}}/.nds],{t,90,100},
BoxRatios->{1.5,1,1},
PlotRange->{{90,100},{0.4,1.6},{0.5,4.}},
AxesLabel->{"t","x[t]","y[t]"},  
PlotStyle->{Directive[Thick,Blue],Directive[Thick,Brown],Directive[Thick,Green],Directive[Thick,Red]},PlotLabel->Style["A Brüsszelátor megoldása,\na megoldás koordinátafüggvényei és\n trajektóriája a fázissíkon","Subsection"],LabelStyle->Directive[Orange,FontFamily->"Helvetica"]]

Szimbolikus megoldás (bemelegítésként)

Nagyjából azt tudja, amit Kamke.

E. Kamke : Differentialgleichungen : Lösungsmethoden und Lösungen, 1. Gewöhnliche Differentialgleichungen, 2. Partielle  Differentialgleichungen

Egyszerű példák

Inhomogén lineáris

ClearAll[a,b,c,egy,ds,kezd1,kezd2];
egy=y'[x]==2x-y[x];kezd1=y[3]==-3;kezd2=y[2]==-3;

First[DSolve[{egy,kezd1},y[x],x]]

ds=Flatten[DSolve[{egy,kezd2},y,x]]

y'[3]/.ds

ClearAll[ds];
ds[c_]:=Simplify[DSolve[{egy,y[1]==c},y[x],x][[1,1,2]]]

ds[c]

DSolve[{y'[x]==a x+b y[x]},y[x],x]

Riccati

ric=DSolve[{y'[x]==(x^2)+((y[x])^2)},y[x],x][[1,1,2]];

Together[FunctionExpand[ric]]

Az elemi függvény fogalmának csak történeti jelentősége van. (HF: Program arra, h elemi-e egy fv?)

Elemi kvalitatív vizsgálat

Iránymező, irányvonalak

Ezt vizsgáljuk: y'(x)=2x-y(x), azaz (OverTilde[x])'(t)=1     (OverTilde[y])'(t)=2(OverTilde[x])(t)-(OverTilde[y])(t)

sp=StreamPlot[{1,2 x-y},{x,-2,4},{y,-2,4}]

l[k_]:=2 x-k;(*irányvonal, amit k meredekséggel metszenek a megoldások*)
iranyok=Plot[{l[-3],l[0],l[1],l[2]},{x,-3,4},PlotRange->{-3,4},PlotStyle->{{Thickness[0.008`],Black},{Thickness[0.008`],Red},{Thickness[0.008`],Black},{Thickness[0.008`],Black}}]

ds[c]

pl[c_]:=Plot[Evaluate[ds[c]],{x,-2,4},PlotStyle->{Thickness[0.008],((1+Sign[c])/2)Green+((1-Sign[c])/2)Blue},PlotRange->{-2,4}];

Show[{pl[2],pl[-3]}]

Show[sp,iranyok,pl[1],pl[2],pl[3],pl[-1],pl[-2],pl[-4],PlotRange->{-2,4},Frame->True,Axes->True]

HF. : Tanulmányozzuk a fenti egyenletet, vagy másokat az EquationTrekker függvény segítségével.

Görbesereg differenciálegyenlete

Eliminate[{y[x]==C (x^3),D[y[x]==C (x^3),x]},{C}]

Eliminate[{(((x-a))^2)+b ((y[x])^2)==1,D[(((x-a))^2)+b ((y[x])^2)==1,x],D[(((x-a))^2)+b ((y[x])^2)==1,{x,2}]},{a,b}]

curvefamilyODE[indep_,eq_,params_]:=Module[{},Eliminate[FoldList[D[#1,indep]&,eq,params],params]]

Philippov 17.

curvefamilyODE[x,y[x]==C (E^(k x)),{C,k}]

Philippov 18. OK, better than in the problem book.

curvefamilyODE[x,y[x]==(((x-C))^3),{C}]

Philippov 19. OK

curvefamilyODE[x,y[x]==C (x^3),{C}]

Philippov 20. OK

curvefamilyODE[x,y[x]==Sin[x+C],{C}]

Eliminate[{y[x]==Sin[x+C],D[y[x]==Sin[x+C],x]},C]

Reduce[{y[x]==Sin[x+C],D[y[x]==Sin[x+C],x]},{C}]

Eliminate[%,C]

Philippov 21. OK

curvefamilyODE[x,(x^2)+C ((y[x])^2)==2y[x],{C}]

Philippov 22. OK

curvefamilyODE[x,((y[x])^2)+C x==(x^3),{C}]

Philippov 23. OK

curvefamilyODE[x,y[x]==C(((x-C))^2),{C}]

Szebb így?

Map[Factor,%,1]

Philippov 24.

curvefamilyODE[x,C y[x]==Sin[C x],{C}]

Upon taking the derivative and the square of the original equation one obtains

(((C y'[x]))^2)==(C^2) ((Cos[C x])^2)

thus

(C^2)==(1-((y'[x])^2))/((y[x])^2.)

Substituting C into the origianl equation one gets

y'[x]==Cos[((x(Sqrt[1-((y'[x])^2)]))/(y[x]))]

How to obtain this?

Philippov 25. OK

curvefamilyODE[x,y[x]==a (x^2)+b (E^x),{a,b}]

Simplify[%]

Hogyan lehet egyenletet egy oldalra rendezni?

Simplify[Subtract@@%]==0

Philippov 26. OK

curvefamilyODE[x,(((x-a))^2)+b ((y[x])^2)==1,{a,b}]

Kijon osztas nelkul is kezzel ?? ??

Philippov 27.

curvefamilyODE[x,Log[y[x]]==a x+b y[x],{a,b}]

FullSimplify[%[[2]]]

y[x]>0 is missing.

Philippov 28. OK

curvefamilyODE[x,y[x]==a (x^3)+b (x^2)+c x,{a,b,c}]

Philippov 29. OK

curvefamilyODE[x,x==a ((y[x])^2)+b y[x]+c,{a,b,c}]

KKM, p. 60, Ex. 1 OK

curvefamilyODE[x, ((x^2)/(a^2))-(((y[x])^2)/1)==1,{a}]

KKM, p. 60, Ex. 2

curvefamilyODE[x, y[x]==a(1-(E^(-(x/a)))),{a}]

a≠0, although we do not construct the whole theory. y'[x]>0

KKM, p. 60, Ex. 3

curvefamilyODE[x, x Cos[α]+y[x] Sin[α]-1==0,{α}]

One may, if (s)he wishes, ask for simplifications with not too impressing results.

Simplify[%[[1]]]

%%[[2]]

How to obtain this:

Solve[Cos[α]+y'[x] Sin[α]==0,y'[x]]

Map[Minus,Reverse[%[[1,1]]]]

Reverse[%%[[1]]]

szab={Cos[α]->((Cot[α])/(Sqrt[1+((Cot[α])^2)])),Sin[α]->(1/(Sqrt[1+((Cot[α])^2)]))}

x Cos[α]+y[x] Sin[α]-1==0/.szab/.Map[Minus,Reverse[%44[[1,1]]]]

FullSimplify[%]

Ezt szeretném:

x y'[x]+y[x]==(Sqrt[1+((y'[x])^(2 ))])

KKM, 238. OK

curvefamilyODE[x, y[x]==a/x,{a}]

KKM, 239. OK

curvefamilyODE[x, (x^2)-((y[x])^2)==a x,{a}]

KKM, 240.

curvefamilyODE[x, y[x]==a (E^(x/a)),{a}]

a≠0, y'[x]>0

KKM, 241.

curvefamilyODE[x, y[x]==C x-C-C (x^2),{C}]

A konyvben:

((y'[x])^2)+y'[x]-x y'[x]+y[x]==0

Egyszerubb? Egyenerteku?

KKM, 242. OK

curvefamilyODE[x, y[x]==(E^x)(a x+b),{a,b}]

Irjuk igy:

d[y_]:=D[y,x]

d[y[x]]

(((d-1))^2)[y[x]]==0

Csak ez mukodjon, de hogyan megy?

KKM, 243.

curvefamilyODE[x, ((y[x])^2)==2C x+(C^2),{C}]

In problem books the equation is divided through by y[x] , excluding some solutions.

KKM, 244. OK

curvefamilyODE[x, y[x]==a (x^2)+b x+c,{a,b,c}]

KKM, 245. OK

curvefamilyODE[x, y[x]==(Subscript[C, 1]) x+((Subscript[C, 2])/x)+(Subscript[C, 3]),{(Subscript[C, 1]),(Subscript[C, 2]),(Subscript[C, 3])}]

KKM, 246.

curvefamilyODE[x, (((x-a))^2)+(((y[x]-b))^2)==1,{a,b}]

If you see a simplifying possibility, you may request it.

Map[Factor,%,1]

KKM, 247. OK

curvefamilyODE[x, y[x]==(Subscript[C, 1])(E^x)+(Subscript[C, 2])(E^(-x)),{(Subscript[C, 1]),(Subscript[C, 2])}]

KKM, 248.

curvefamilyODE[x, y[x]==a Sin[x+α],{a,α}]

We received more possibilities which have nothing to do with the solution; that is given as the first term.

Matveev, 62. OK

curvefamilyODE[x,y[x]==(1/3)(x^3)+C,{C}]

Matveev, 63.

curvefamilyODE[x,y[x]==(Sqrt[1-(x^2)])+C,{C}]

One would be happier to see a warning that x≠1, or that -1<x<1.

Matveev, 64.

curvefamilyODE[x,y[x]==Log[x]+C,{C}]

The warning x>1 would be more appropriate from the mathematical point of view—if one only considers real numbers.

Egyszerű típusok, pl.: egzakt egyenletek

Ezt akarjuk megoldani:

ex=(4(x^3)((y[x])^3)-2x y[x])+(3(x^4)((y[x])^2)-(x^2))y'[x]==0;

P[x_,y_]:=ex[[1,{1,2}]]/.{y[x]->y}
Q[x_,y_]:=((ex[[1,3]]/.{y[x]->y})/(y'[x]))

{P[x,y],Q[x,y]}

Keresztbevett parciálisok (+tartomány!!)

D[P[x,y],y]==D[Q[x,y],x]

cpP=ContourPlot[P[x,y]==0,{x,-4,4},{y,-1,1},ContourStyle->Directive[Red,Thick]]

Q[x,y]

cpQ=ContourPlot[-(x^2)+3 (x^4) (y^2)==0,{x,-4,4},{y,-1,1}]

Show[cpP,cpQ](*Egér!*)

cpQ=ContourPlot[Q[x,y]==0,{x,-4,4},{y,-1,1},ContourStyle->Directive[Blue,Thick]]

A primitív függvény meghatározása

A megoldás, F[x,y]=állandó, szintvonalak

F[x_,y_]=Integrate[Q[x,y],y]+Integrate[P[x,y]-Integrate[D[Q[x,y], x], y],x]

Az (1,1) ponton átmenő megoldások:

Union[Solve[F[x,y]==F[1,1],y]]

Plot3D[Evaluate[F[x,y]],{x,-4,4},{y,-1,1}]

(*néhány szintvonal*)
ContourPlot[F[x,y],{x,-5,5},{y,-5,5},Contours->20,ContourShading->False]

ContourPlot[F[x,y],{x,-1.5`,1.5`},{y,-1.5`,1.5`},Contours->20,ContourShading->False,PlotPoints->50]

ContourPlot[F[x,y],{x,-1.5`,0},{y,0,1.5`},Contours->20,ContourShading->False,PlotPoints->50]

Integráló tényező

Lehetne ezt is. HF

Közelítő módszerek

A pontos megoldás

ClearAll[pontos,y,x];
pontos=DSolve[{y'[x]==x y[x], y[0]==2},y[x],x][[1]]

Euler-módszer (vö. általános tételek)

NestList[ff,x0,2]

euler[f_,h_,x0_,y0_,n_:3]:=Module[{lep,x,y},lep[{x_,y_}]:={x+h,y+h f[x,y]};NestList[lep,{x0,y0},n]]

ClearAll[F];
F[x_,y_]:=x y;

lp=ListPlot[euler[F,0.01`,0,1,100],Joined->True,PlotStyle->{RGBColor[1,0,0]}]

pont=Plot[y[x]/.pontos,{x,0,1},PlotLabel->"pontos"]

Show[pont,lp]

abra[n_]:=Module[{},lp=ListPlot[euler[F,((1.`)/n),0,1,n],Joined->True,PlotStyle->{RGBColor[1,0,0]}];pl=Plot[2(E^((x^2)/2)),{x,0,1}];
Show[pl,lp,PlotLabel->"Lépésszám: "<>ToString[n]]]

GraphicsGrid[{{abra[1],abra[100]},{pont}},ImageSize->300]

Megoldás hatványsorokkal

Az előző példa

ClearAll[x,y,eq,ini];eq=y'[x]==x y[x]; ini= y[0]==2;

y[x_,n_:5]:=2+(∑_(k=1)^n)(Subscript[a, k])(x^k)+((O[x])^(n+1))

le=LogicalExpand[eq]

sol=Solve[le][[1]]

y[x]/.sol//Normal

Series[2(E^((x^2)/2)),{x,0,4}]

Bonyolultabb példa

HF: Feltételek? Mikor létezik, hol konvergens?

ClearAll[hatv,y];
hatv[eq_,inival_,n_]:=Module[{},y[x_]:=inival+(∑_(k=1)^n)(Subscript[a, k])(x^k)+((O[x])^(n+1));

y[x]/.Solve[LogicalExpand[eq]]]

hatv[y'[x]==x y[x],2,4]

h=hatv[((y'[x])^2)-y[x]==x,1,3]

koz=Plot[Evaluate[Normal[h]],{x,-0.5,2}]

ClearAll[y,x,nex];nex=NDSolve[{((y'[x])^2)-y[x]==x,y[0]==1},y,{x,-0.5,2},MaxSteps->(10^5)];

pl=Plot[y[x]/.nex,{x,-0.5,2},PlotStyle->Directive[Red,Thick]]

Show[koz,pl]

HF: Szimbolikus megoldás?

Composition, InverseSeries

Magasabbrendű Newton:

InverseSeries[Series[f[x],{x,a,3}]]//Normal/.{x->0}

ComposeSeries[Series[Exp[x],{x,0,10}],Series[Sin[x],{x,0,10}]]

Series[Exp[Sin[x]],{x,0,10}]

Taylorsor-módszer

eq=y'[x]==x y[x]; ini= y[0]==2;

HF

Szukcesszív approximáció (vö. általános tételek)

Az eszköz

NestList[φ,(Subscript[x, 0]),3]

nl=NestList[Cos,-6.7,12]

ListPlot[nl,PlotStyle->{PointSize[0.02`]}]

Alkalmazzuk a differenciálegyenletre

Az egyenértékű integrálegyenlet:

(∫_0^x)y'[ξ]\[DifferentialD]ξ==(∫_0^x)ξ y[ξ]\[DifferentialD]ξ/.{y[0]->1}

Az operátor, amit iterálunk:

A[φ_]:=1+Integrate[x φ, x]

nl=NestList[A , 1,2]

Összehasonlítás

ClearAll[pontos,y,x];
pontos=DSolve[{y'[x]==x y[x], y[0]==2},y[x],x][[1]]

Plot[Evaluate[{y[x]/.pontos,nl[[3]]}],{x,0,4},PlotStyle->{Red,Blue}]

Találjuk ki a hatványsor összegét

Becslések az értelmezési tartományra?

NDSolve: Határozatlan integrál numerikusan

sol = NDSolve[{y'[x] == Sin[Exp[x]],y[0] == 0},  y, {x, 0, 3}]

y[3]/.sol

First[y[3]/.sol]

Plot[Evaluate[y[x] /. sol], {x,0,3}, Epilog->Point[{0,0}]]

Megj. merev egyenletek

Alkalmazások

Csak hogy legyen kedvünk tovább nézni...

Volterra-Lotka

vl=NDSolve[{y'[t]==y[t](x[t]-1),x'[t]==x[t](2-y[t]),x[0]==1,y[0]==2.7},{x,y},{t,0,10}]

pp=ParametricPlot[Evaluate[{x[t],y[t]}/.First[vl]],{t,0,10},Epilog->{PointSize[Medium],Point[{1,2}]}]

sp=StreamPlot[{y(x-1),x(2-y)},{x,0,1.6},{y,1.4,3}]

stp=StreamPlot[{y(x-1),x(2-y)},{x,0,4},{y,0,4}]

Show[pp,sp]

Lorenz

s=NDSolve[{x'[t]==-3 (x[t]-y[t]),y'[t]==-x[t] z[t]+26.5 x[t]-y[t],z'[t]==x[t] y[t]-z[t],x[0]==z[0]==0,y[0]==1},{x,y,z},{t,0,200},MaxSteps->∞]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s],{t,0,200},PlotPoints->1000,ColorFunction->(Hue[#4]&)]

Gavrilov-Shilnikov

eqns = {
x'[t]==x[t](2+μ-10(((x[t])^2)+((y[t])^2))) + ((z[t])^2) + ((y[t])^2) +2 y[t],
y'[t]==-((z[t])^3)-(1+y[t])(((z[t])^2)+((y[t])^2)+2 y[t])-4 x[t]+μ y[t],
z'[t]==(1+y[t]) ((z[t])^2) +((x[t])^2)-ϵ};

Table[Block[{ϵ=0.0357,T=1000},
sol = First[NDSolve[{eqns, x[0]==y[0]==z[0]==1},{x,y,z},{t,0,T}, MaxSteps->+∞]];
ParametricPlot3D[Evaluate[{x[t],y[t],z[t]} /. sol],{t,T/4,T}, BoxRatios->1,PlotPoints->1000,PlotRange->All]],{μ,.455,.456,.001}]

Hővezetés, periodikusan fűtjük a talajt

nho=NDSolve[{D[u[t,x],t]==D[u[t,x],x,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Plot3D[Evaluate[u[t,x]/.nho],{t,0,10},{x,0,5}, PlotRange->All,AxesLabel->{t,x,u}]

Hullámegyenlet

s=NDSolve[{(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x],u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10]},u,{t,0,40},{x,-10,10}]

DensityPlot[Evaluate[First[u[t,x]/.s]],{t,0,40},{x,-10,10},PlotPoints->30]

n-testproblémák

Redukált 3 testprobléma

eqn=z''[t]==-z[t]/(((((z[t])^2)+((((1/2) (1+0.1 Sin[2 π t])))^2)))^(3/2));

nred=NDSolve[{eqn,z[0]==1,z'[0]==0},z,{t,0,30}];

Plot[Evaluate[z[t]/.nred],{t,0,30}]

Több kezdeti feltétellel

ClearAll[eqn,z,s,sol];eqn=D[z[t,s],{t,2}]==-z[t,s]/(((((z[t,s])^2)+((((1/2) (1+0.1 Sin[2 π t])))^2)))^(3/2));

sol = NDSolve[{eqn,z[0,s]==s,Derivative[1,0][z][0,s]==0},z,{t,0,10}, {s,1,2}]

Plot3D[Evaluate[z[tp, z0]/.sol],{tp,0,10},{z0,1,2},AxesLabel->{t,zkezd}]

Nemlineáris Schrődinger

A soliton profile perturbed by a periodic potential:

L=50;
Framed[TraditionalForm[Column[{I D[u[t,x],t]+D[u[t,x],x,x]+2 ((Abs[u[t,x]])^2) u[t,x]==0.1 (1-Cos[((π x)/L)]),u[0,x]==Sech[x] Exp[I x],u[t,-L]==u[t,L]}]]]

sol=NDSolve[{I D[u[t,x],t]+D[u[t,x],x,x]+2 ((Abs[u[t,x]])^2) u[t,x]==0.1 (1-Cos[((π x)/L)]),u[0,x]==Sech[x] Exp[I x],u[t,-L]==u[t,L]},u,{t,0,10},{x,-L,L}]

Plot3D[Evaluate[Abs[u[t,x]/.First[sol]]],{t,0,10},{x,-L,L},PlotPoints->60,MaxRecursion->3,Mesh->None]

Lineáris rendszerek és magasabbrendűek

Mátrix exponenciális függvénye

MatrixForm[A={{1,-1,1},{1,1,-1},{2,-1,0}}];
ini={1,-1,0};

ds=DSolve[Thread/@{{x'[t],y'[t],z'[t]}==A.{x[t],y[t],z[t]},{x[0],y[0],z[0]}==ini},{x[t],y[t],z[t]},t][[1]]

me=MatrixExp[A t].ini

({x[t],y[t],z[t]}/.ds)-me//Simplify

A fázissík az origó körül

Elmélet

ClearAll[sp,x,y];sp[λ_,μ_,opts___]:=StreamPlot[{λ x,μ y},{x,-3,3},{y,-3,3},opts]

i=sp[1,2,PlotLabel->"λ>0, μ>0\nInstabilis csomó"];
s=sp[-1,-2,PlotLabel->"λ<0, μ<0\nSabilis csomó"];
n=sp[1,-2,PlotLabel->"λ>0, μ<0\nNyereg"];

GraphicsGrid[{{i,s},{n}},ImageSize->400]

elf[λ_,opts___]:=StreamPlot[{λ x+y,λ y},{x,-3,3},{y,-3,3},opts]

ie=elf[1,PlotLabel->"λ>0\nInstabilis elfajult csomó"];
se=elf[-1,PlotLabel->"λ<0\nSabilis elfajult csomó"];

GraphicsRow[{ie,se},ImageSize->400]

kom[α_,β_,opts___]:=StreamPlot[{α x-β y,β x+α y},{x,-3,3},{y,-3,3},opts]

if1=kom[1,-2,PlotLabel->"α>0\nInstabilis fókusz"];
if2=kom[1,2,PlotLabel->"α>0\nInstabilis fókusz"];
if3=kom[1,0,PlotLabel->"α>0\nInstabilis fókusz"];
sf=kom[-1,-2,PlotLabel->"α<0\nSabilis fókusz"];
cf=kom[0,2,PlotLabel->"α=0\nCentrum"];
GraphicsGrid[{{if1,if2},{if3,sf},{cf}},ImageSize->400]

sp[λ_,μ_,opts___]:=StreamPlot[{λ x,μ y},{x,-3,3},{y,-3,3},opts]

i=sp[0,2,PlotLabel->"λ=0, μ>0"];
s=sp[0,-2,PlotLabel->"λ=0, μ<0"];
n=sp[0,0,PlotLabel->"λ=0, μ=0"];
masik=StreamPlot[{y,0},{x,-3,3},{y,-3,3}];

GraphicsGrid[{{i,s},{n,masik}},ImageSize->400]

Első példa, kézzel

eq={x'[t]==-2x[t]+y[t],y'[t]==x[t]-y[t]};

rhs = eq[[All, 2]] /. z_[t] -> z

var={x,y}

MatrixForm[A=D[rhs,{var}]/.{x->0,y->0}]

Eigenvalues[A]

N[%]

Stablis csomó

Program

Thread[{x,y}->0]

stabQ[rhs_,var_:{x,y}]:=Eigenvalues[D[rhs,{var}]]

stabQ[{-2 x+y,x-y}]

sokrhs={{-2x+2y,x-y},{-x+2y,x+y},{3x,2x+y},{x+3y,-6x-5y},{-4x+(1/2)y,2x-4y},{-2x-5y,2x+2y},{x+y,y},{x+2y,-2x+y}};

tipus={elfajult,nyereg,"instabilis\ncsomó","stabilis\nfókusz","stabilis\ncsomó","centrum","intabilis\nelfajult\ncsomó","intabilis\nfókusz"};

Grid[{sokrhs,stabQ/@sokrhs,tipus},Frame->All]

GraphicsGrid[Partition[StreamPlot[#,{x,-3,3},{y,-3,3},ImageSize->200]&/@sokrhs,2]]

Itt is kereshetünk Ljapunov-függvényt, kvadratikusat.

Inga: pendule simple

Wronski-determináns

DSolve[y''[x]-x y[x]==0,y,x]

Wronskian[{AiryAi[x],AiryBi[x]},x]

Peremértékfeladatok (nemlin)

Plot[Evaluate[{y[x], y'[x]} /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 1,y[10] == 0}, y,x]],{x,0,10}]

Table[Plot[Evaluate[{y[x], y'[x]} /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 0,y[10] == 0}, y,x, Method->{"Shooting", "StartingInitialConditions"->{y[0] == 0, y'[0] ==i}}]],{x,0,10}],{i,1,2}]

Nemlineáris rendszerek

Verifikálás: nemlin dersz megoldását

1. Vizsgáljuk meg, hogy az x'=y     y'=((y^2)/x)  egyenletrendszernek megoldása-e a φ(t):=(e^(2t))     ψ(t):=(e^t) függvénypár!

Először is az egyenletrendszer jobb oldalának értelmezési tartományaként Ω:=ℝ× (SuperPlus[ℝ])× ℝ vehető (miért?). Azután, ha

φ[t_]:=(E^(2t));ψ[t_]:=(E^t);

akkor

{φ'[t],ψ[t],ψ'[t],(((ψ[t])^2)/(φ[t]))}

tehát nem megoldás. Azt majd alább megadjuk.

Belátjuk, h első integrál, amit megadtunk

Bizonyítsuk be, hogy a   ψ(t,p,q):=arctg((p/q))-t képlettel megadott függvény első integrálja az x'=((x^2)/y)     y'=-((y^2)/x) egyenletnek!

Először is az egyenletrendszer jobb oldalának értelmezési tartományaként Ω:=ℝ× (SuperPlus[ℝ])× (SuperPlus[ℝ]) vehető (miért?). Azután, ha

ψ[t_,p_,q_]:=ArcTan[(p/q)]-t;f[p_,q_]:={((p^2)/q),-((q^2)/p)}

akkor, mivel:

D[ψ[t,p,q],{{t,p,q}}]//Simplify

ezért

{-1,(q/((p^2)+(q^2))),-(p/((p^2)+(q^2)))}.{1,((p^2)/q),-((q^2)/p)}//Simplify

Használjuk fel az első integrált a rendszer megoldásához!

arctg(((x(t))/(y(t))))-t=(Subscript[c, 1]), tehát x(t)=y(t)tg(t+(Subscript[c, 1])), amit a másodikba írva kapjuk y'(t)=-((y(t))/(tg(t+(Subscript[c, 1])))), ezt a lineáris egyenletet már meg tudjuk oldani.

Megoldunk néhányat, sajnos, trükkökkel, de első integrál néha ránézésre adódik

Oldjuk meg az x'=y     y'=((y^2)/x)  egyenletrendszert!

Mivel x y'=(y^2) és x'y=(y^2), ezért ((y/x))'=0, tehát (y/x)=(Subscript[c, 1]), amit az első egyenletbe írva kapjuk, hogy x'=(Subscript[c, 1]) x, tehát x(t)=(Subscript[c, 2])(e^((Subscript[c, 1])t)). Ezt az első egyenletbe téve kapjuk, hogy y(t)=(Subscript[c, 1])(Subscript[c, 2])(e^((Subscript[c, 1])t)).

TraditionalForm[Column[ds=DSolve[{x'[t]==y[t],y'[t]==(((y[t])^2)/(x[t]))},{x[t],y[t]},t][[1]]]]

Differenciál-algebrai rendszerek

Index stb. Szinguláris perturbáció

DAE = {(Subscript[x, 1])'[t]==(Subscript[x, 3])[t],(Subscript[x, 2])[t] (1-(Subscript[x, 2])[t])==0,(Subscript[x, 1])[t] (Subscript[x, 2])[t]+(Subscript[x, 3])[t] (1-(Subscript[x, 2])[t])==t};

Find the solution with (Subscript[x, 1])'(0)==1:

sol = NDSolve[{DAE,(Subscript[x, 1])'[0]== 1},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t,0,1}]

Stabilitás

Stabilitás explicit megoldásból a definíció alapján

1. Vizsgáljuk meg az x'(t)=((sin(x(t)))^2)egyenlet stacionárius pontjainak stabilitását!

A stacionárius pontok:  (Subscript[x, k]^( * ))=k π (k∈ℤ). Mivel a kezdetiérték-problémák megoldása (Subscript[x, k])(t)=arcctg(ctg((Subscript[x, 0]))-t)+k π, ezért minden megoldás korlátos, de látható, hogy a megoldások a (2k+1)π értékhez tartnak alulról, a 2k π értéktől meg távolodnak.

2.  Vizsgáljuk meg az x'=y     y'=x egyenlet stacionárius pontjainak stabilitását!

Az egyetlen stacionárius pont az origó. Mivel az x'=y     y'=x     x(0)=(Subscript[x, 0])    y(0)=(Subscript[y, 0]) kezdetiérték-probléma megoldása

TraditionalForm[Column[ds=DSolve[{x'[t]==-y[t], y'[t]== x[t],     x[0]== x0,    y[0]== y0},{x[t],y[t]},t][[1]]]]

ez az egész számegyenesen értelmezve van, és minden ϵ-hoz van olyan δ, hogy az origóhoz ϵ-nál közelebb induló megoldások nem távolodnak el jobban az origótól, mint δ:

|x0 cos(t)-y0 sin(t)|<|x0 cos(t)|+|y0 sin(t)|<|x0|+|y0|

|x0 sin(t)+y0 cos(t)|<|x0 sin(t)|+|y0 cos(t)|<|x0|+|y0|

Stab nemlin rsznél, linearizálással

Nem megy

eq={x'[t]==-((x[t])^3)-y[t],y'[t]==x[t]-((y[t])^3)};

rhs=eq[[All,2]]/.z_[t]->z

MatrixForm[A=D[rhs,{{x,y}}]/.{x->0,y->0}]

Eigenvalues[A]

Épp nem dönthető el.

StreamPlot[rhs,{x,-3,3},{y,-3,3}]

StreamPlot[rhs,{x,-.3,.3},{y,-.3,.3}]

Gyanús! Bendixson mutatja, h az első negyedben nem lehet periodikus megoldás. HF!

És most?

rhs={(x^3)-y+Sin[x+y],x+(y^3)+Sin[x]};

Solve[Thread[rhs==0],{x,y}]

FindRoot[Thread[rhs==0],{{x,1},{y,-1}}]

FindInstance[Thread[rhs==0],{x,y}]

Belátható, h nincs másik stac pont. HF!

MatrixForm[A=D[rhs,{{x,y}}]/.{x->0,y->0}]

Eigenvalues[A]

Instabilis, mert van pozitív valós részű sajátértéke.

StreamPlot[rhs,{x,-3,3},{y,-3,3}]

Megadott Ljapunov függvénnyel

eq={x'[t]==-((x[t])^3)-y[t],y'[t]==x[t]-((y[t])^3)};

rhs=eq[[All,2]]/.z_[t]->z

v[x_,y_]:=(x^2)+(y^2)

Expand[D[v[x,y],{{x,y}}].rhs]

StreamPlot[rhs,{x,-3,3},{y,-3,3}]

Keressünk Ljapunov-függvényt

1.

eq={x'[t]==2((y[t])^3)-((x[t])^5),y'[t]==-x[t]-((y[t])^3)+((y[t])^5)};

rhs=eq[[All,2]]/.z_[t]->z

Select[Union[NSolve[rhs=={0,0},{x,y}]],FreeQ[#[[-1]],Complex]&]

Csak az origót fogjuk vizsgálni.

Legyen

v[x_,y_]:=(x^α)+(y^β)

Expand[D[v[x,y],{{x,y}}].rhs]

Tűnjenek el a vegyes tagok:

Reduce[{2 α==β,α-1==1,β-1==3},{α,β},Integers]

joparam={α->2,β->4};

Expand[D[v[x,y],{{x,y}}].rhs]/.joparam

Azonos ezzel:

-2((x^6)+2(y^6)(1-(y^2)))

és ez negatív az origóhoz közel.

StreamPlot[rhs,{x,-3,3},{y,-3,3}]

És valóban csak lokális a stabilitás! Jól látszik a másik két egyensúlyi helyzet, nyeregpontnak tűnnek, nem?

2.

eq={x'[t]==x[t]-y[t]-((x[t])^3),y'[t]==x[t]+y[t]-((y[t])^3)};

rhs=eq[[All,2]]/.z_[t]->z

Legyen

v[x_,y_]:=(x^α)+(y^β)

Expand[D[v[x,y],{{x,y}}].rhs]

Tűnjenek el a vegyes tagok:

Reduce[{ α==β,α-1==1,β-1==1},{α,β},Integers]

joparam={α->2,β->2};

Expand[D[v[x,y],{{x,y}}].rhs]/.joparam

Azonos ezzel:

2((x^2)(1-(x^2))+(y^2)(1-(y^2)))

és ez pozitív az origóhoz közel:

{ContourPlot[-2((x^2)(1-(x^2))+(y^2)(1-(y^2))),{x,-3,3},{y,-3,3}],Plot3D[-2((x^2)(1-(x^2))+(y^2)(1-(y^2))),{x,-3,3},{y,-3,3}]}

StreamPlot[rhs,{x,-3,3},{y,-3,3}]

És valóban instabilis, de mintha lenne egy periodikus megoldás! A Bendixson-tétel nem is zárja ki:

Tr[D[rhs,{{x,y}}]]

Bifurkáció I.

A

eq={x'==-3x+α y,y'==2x+y};

Manipulate[StreamPlot[{-3x+α y,2x+y},{x,-3,3},{y,-3,3}],{{α,-3/2},-2,-1,0.1}]

Eigenvalues[D[{-3x+α y,2x+y},{{x,y}}]/.{x->0,y->0}]

CharacteristicPolynomial[{{-3,α},{2,1}},λ]

Plot[{-1-(Sqrt[2]) (Sqrt[2+α]),-1+(Sqrt[2]) (Sqrt[2+α])},{α,-2,3}]

Bifurkáció II. (Andronov-Hopf)

Manipulate[DynamicModule[{nds,x,y},nds=NDSolve[{x'[t]==1+a ((x[t])^2)y[t]-(5+1)x[t],y'[t]==a(5 x[t]-a ((x[t])^2)y[t]),x[0]==8,y[0]==3},{x,y},{t,0,100}][[1]];
Dynamic@ParametricPlot3D[Evaluate[{{90,x[t],y[t]},{t,x[t],0.5},{t,1.95,y[t]}}/.nds],{t,90,100},BoxRatios->{1,1,1},PlotRange->{{90,100},{0.4,1.95},{0.5,4.5}},AxesLabel->{"t","x[t]","y[t]"},  PlotStyle->Directive[Thick,Red]]],{a,2.03,1.9,0.001,Appearance->"Labeled"},LocalizeVariables->False]

Paraméterbecslés

sol=First[x/.NDSolve[{x''[t]+.33 x'[t]+.72x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]];
times=N[Range[0,100]/5];
data=Transpose[{times,sol[times]+RandomReal[.1,101]}];
lp=ListPlot[data,PlotRange->All]

model[γ_?NumberQ,a_?NumberQ,b_?NumberQ]:=(model[γ,a,b]=First[x/.NDSolve[{x''[t]+γ x'[t]+a x[t]+b x[t]^3==0,x[0]==2,x'[0]==0},x,{t,20}]]);

fit=FindFit[data,model[γ,a,b][x],{{γ,.1},{a,.1},{b,1}},x,PrecisionGoal->4,AccuracyGoal->4]

{Show[Plot[model[γ,a,b][x]/.fit,{x,0,20}],lp],ListPlot[Transpose[{times,data[[All,2]]-(model[γ,a,b][times]/.fit)}]]}

Fit

Késleltetett argumentumú egyenletek

Mackey-Glass

Manipulate[
Module[{sol,x,t},
sol = First[ NDSolve[{x'[t] == a ((x[t - τ])/(1 + ((x[t - τ])^10))) -b  x[t], x[t /; t ≤ 0]==1/2},x,{t,0,500}]];
ParametricPlot[Evaluate[{x[t], x[t - τ]} /. sol],{t,300,500}, PlotRange->{{0,2},{0,2}}]],{{a,.2},0,1},{{b,.1},0,1},
{{τ,17},1,20}]

Diszkrét dinamikai rendszerek

Innen kezdődik

Logistic[r_]:=Function[x,r x(1-x)];Exponential[r_]:=Function[x,x (E^(r(1-x)))];
Orbit[map_,x0_,n_]:=NestList[map,x0,n];

IterativeProcess[map_,x0_,{min_,max_},opts___]:=Module[{fr,orb},
orb=Orbit[map,x0,50];
fr=MapThread[Line[{{#1,#1},{#1,#2},{#2,#2}}]&,{Drop[orb,-1],Drop[orb,1]}];
Show[Plot[{map[x],x},{x,min,max},opts],Graphics[{fr}]]];

GraphicsArray[{ListPlot[Orbit[Logistic[2.6],0.1,100],AxesLabel->{r,"(SuperStar[N])"},PlotLabel->"r = 2.6"],IterativeProcess[Logistic[2.6],0.1,{0,1},AxesLabel->{N,"r N(1-N)"},PlotLabel->"r = 2.6"]},ImageSize->700]

Ez a szép

abrapar[r_]:=GraphicsArray[{ListPlot[Orbit[Exponential[r],0.1,100],AxesLabel->{"r","(SuperStar[N])"},PlotLabel->"r = "<>ToString[r]],IterativeProcess[Exponential[r],0.1,{0,2},AxesLabel->{"N","((NE)^(r(1-N)))"},PlotLabel->"r = "<>ToString[r]]},ImageSize->700]

abrapar[2.0]

abrapar[2.2]

abrapar[2.53]

abrapar[2.65]

abrapar[2.67]

abrapar[3.3]

Újra logisztikus

GraphicsArray[{ListPlot[Orbit[Logistic[3.2],0.1,100]],IterativeProcess[Logistic[3.2],0.1,{0,1}]}]

GraphicsArray[{ListPlot[Orbit[Logistic[3.9],0.1,100]],IterativeProcess[Logistic[3.9],0.1,{0,1}]}]

Bifurcation[map_,{p_,pi_,pf_,np_},n_]:=Module[{pp},ListPlot[Flatten[Table[({pp,#})&/@Drop[Orbit[map/.p->pp,0.5,n+50],50],{pp,pi,pf,((pf-pi)/np)}],1],Axes->False,Frame->True,PlotStyle->{PointSize[0.003]}]]

Bifurcation[Logistic[r],{r,0,4,100},2000]

Show[Graphics[Rectangle[{0,0},{1,1}]],Bifurcation[Logistic[r],{r,-2,4,250},200],Graphics[Rectangle[{0.3,0.5},{0.7,1}]],Bifurcation[Logistic[r],{r,0,4,100},200]]

Show[Graphics[{Rectangle[{0,0},{1,1}],Rectangle[{0.3,0.5},{0.7,1}]}],Bifurcation[Logistic[r],{r,-2,4,250},200],Bifurcation[Logistic[r],{r,0,4,100},200]]

Parciális differenciálegyenletek

Hővezetés

nds=NDSolve[{D[u[t,x],t]==D[u[t,x],{x,2}],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Plot3D[Evaluate[u[t,x]/.nds],{t,0,10},{x,0,5}, PlotRange->All,AxesLabel->{t,x,"u  "},LabelStyle->Medium]

nds=NDSolve[{D[u[t,x],t]==D[u[t,x],{x,2}],u[t,0]==0,u[0,x]==Sin[x]},u,{t,0,10},{x,0,2π}]

Plot3D[Evaluate[u[t,x]/.nds],{t,0,10},{x,0,2π}, PlotRange->All,AxesLabel->{t,x,"u  "},LabelStyle->Medium]

Még hő

hpde={D[u[t,x],t] ==  D[u[t,x],x,x], u[0,x] == UnitStep[1/2 - x], u[t,0] == 1, u[t,1] == 0};

Discontinuities in the initial value may result in too many spatial grid points:

Plot3D[Evaluate[u[t,x] /. NDSolve[hpde, u, {t,0,.01},{x,0,1}]],{x,0,1}, {t,0,.01}]

NDSolve

Variációszámítás

<<VariationalMethods`

Euler-Lagrange az ingára

EulerEquations[(1/2) m (l^2) ((θ'[t])^2)+m g l Cos[θ[t]],θ[t],t]

DSolve[%,θ[t],t]

ClearAll[r,φ,U];FirstIntegrals[(1/2) m (((r'[t])^2)+((r[t])^2) ((φ'[t])^2))-U[r],{r[t],φ[t]},t]

Most akkor tanuljuk az elméletet vagy ne?

Tiszta fogalmak (vö. mozi és könyv)

Előkészítjük a gépnek (transzformációk)

Finomodó lépésköz??

(τ/(h^2))<C

HF: a zárthelyi feladatok megoldása