Ugrás a tartalomhoz

Impulzív jelenségek modelljei

Karsai János

Typotex

10.2. A közönséges differenciálegyenletek esete

10.2. A közönséges differenciálegyenletek esete

Tekintsük az

x'=f(t,x), x?R^n
x ' = f ? ( t , x ) , x ? R n

differenciálegyenletet, ahol f(t,0)=0

f ? ( t , 0 ) = 0 . Az egyenlet x(t,t_0,x_0) x ? ( t , t 0 , x 0 ) megoldásait valamely, az origót belsejében tartalmazó G tartomány ?G ? G határpontjaiból indulva (síkon az origó körüli zárt görbe) határozzuk meg. Most csupán a módszer kísérleti támogatásával foglalkozunk. Rouche, Habets és Laloy [24] könyve mélyebb elméleti és gyakorlati ismereteket nyújt a Ljapunov módszerrel kapcsolatosan.

A módszer lépései: a lineáris fékezett oszcillátor

A Ljapunov módszer kísérleti lépéseit (mint korábban más módszereket is) a jól ismert lineáris fékezett oszcillátor példáján mutatjuk meg (lásd a 6.2. és 13.1. fejezeteket). Az alkalmazások során mindig a leghatékonyabb vizualizációs eszközöket alkalmazzuk. Az egyenlet 2D rendszer alakjában

x^'=y, y^'=–x–ay, ?///>///0. x ' = y , y ' = - x - a ? y , ? ///>/// 0. (10.2.1)

A rendszer teljes energiáját a

V(x,y)=(1/2)(x^2+y^2)
V ? ( x , y ) = 1 2 ? ( x 2 + y 2 )

függvény írja le. Ezért a megoldásokat az energia valamely szintvonalából indítjuk. A rendszer lineáris, tetszőleges körvonal, az egységsugarú kör is megfelelő.

  • A rendszer megadása

xyvar={x,y}; A=(0 1 / –1 –a); eqnparm={a›1};

xyvar = { x , y } ; A = ( 0 1 - 1 - a ) ; eqnparm = { a 1 } ;

xydot=Flatten[A.(x / y)];

xydot = Flatten [ A . ( x y ) ] ;

t0=0;t1=6;dt=1; x1=y1=–1.5; x2=y2=1.5;

t0 = 0 ; t1 = 6 ; dt = 1 ; x1 = y1 = - 1.5 ; x2 = y2 = 1.5 ;

  • A Ljapunov függvény megadása és tulajdonságai

V[x_,y_]:=(y^2/2)+(x^2/2);

V [ x_ , y_ ] := y 2 2 + x 2 2 ;

DV=Simplify[Grad[V[x,y],xyvar].xydot]/.eqnparm

DV = Simplify [ Grad [ V [ x , y ] , xyvar ] . xydot ] /. eqnparm

–y^2

- y 2

A Ljapunov függvény rendszer szerinti deriváltjára az ismert eredményt kaptuk. A derivált csak negatív szemidefinit, de akár a sajátértékek negativitása, akár a La Salle-féle invarianciaelv [24] biztosítja a zérómegoldás aszimptotikus stabilitását.

  • A V(x,y) függvény szintvonalai és a rendszer vektormezője E

    E

Egy általános kép kialakításához ábrázoljuk a V(x,y)=0.2,0.5,1

V ? ( x , y ) = 0.2 , 0.5 , 1 szintvonalakat és a vektormezőt.

level={0.2,0.5,1};

level = { 0.2 , 0.5 , 1 } ;

pV=ContourPlot[V[x,y],{x,x1,x2},{y,y1,y2}, Contours›level,ContourShading›False];

pV = ContourPlot [ V [ x , y ] , { x , x1 , x2 } , { y , y1 , y2 } , Contours level , ContourShading False ] ;

plfield=PlotVectorField[xydot/.eqnparm,{x,x1,x2},{y,y1,y2}];

plfield = PlotVectorField [ xydot /. ? eqnparm , { x , x1 , x2 } , { y , y1 , y2 } ] ;

Show[plfield,pV,Axes›True];

Show [ plfield , pV , Axes True ] ;

Jól látható, hogy a vektorok a szintvonalakat kívülről befelé metszik, kivéve ha y=0

y = 0 , ahol a vektorok érintik őket. Sokkal jobban vizsgálható a szintvonalak és a vektorok viszonya, ha a szintvonalakból induló vektorokat ábrázoljuk. A szintértékek origóhoz való közelítésével az alábbi utasítások segítségével igen informatív animációt kaphatunk. E E

pV1=Table[ContourFieldPlot2D[V[x,y]–level[[i]],xydot/.eqnparm, {x,x1,x2},{y,y1,y2},{{PlotPoints›20},{ScaleFactor›0.4}, {AspectRatio›Automatic,PlotRange›{{x1,x2},{y1,y2}}}}], {i,1,Length[level]}];

pV1 = Table [ ContourFieldPlot2D [ V [ x , y ] - level [ i ] , xydot /. ? eqnparm , { x , x1 , x2 } , { y , y1 , y2 } , { { PlotPoints 20 } , { ScaleFactor 0.4 } , { AspectRatio Automatic , PlotRange { { x1 , x2 } , { y1 , y2 } } } } ] , { i , 1 , Length [ level ] } ] ;

Show[pV1];

Show [ pV1 ] ;

  • A rendszer megoldása

A rendszert a V(x, y)=1

V ? ( x , y ) = 1 körvonal pontjaiban oldjuk meg:

Icond=Table[?2N[{Cos[u],Sin[u]},2],{u,0,2?,(?/8)}];

Icond = Table [ 2 ? N [ { Cos [ u ] , Sin [ u ] } , 2 ] , { u , 0 , 2 ? ? , ? 8 } ] ;

Traj[t_]=ODESolve[xydot/.eqnparm,xyvar,Icond,{t,t0,t1}];

Traj [ t_ ] = ODESolve [ xydot /. eqnparm , xyvar , Icond , { t , t0 , t1 } ] ;

  • A Ljapunov függvény változása a megoldások mentén

Plot[Evaluate[(V[#1[[1]],#1[[2]]]///&///)/@Traj[t]],{t,t0,t1}, PlotRange›{0,1}];

Plot [ Evaluate [ ( V [ #1 [ 1 ] , #1 [ 2 ] ] ///&/// ) /@ Traj [ t ] ] , { t , t0 , t1 } , PlotRange { 0 , 1 } ] ;

  • A megoldások és a Ljapunov függvény ábrázolása E

    E

Az alábbi utasítás eredményeként kapottak közül csak két lineárisan független megoldás grafikonját láthatjuk.

SolPlot[Traj[t]/.{u_,v_}›{u,v,V[u,v]},{t,t0,t1}, {x,y,V},PlotRange›All];

SolPlot [ Traj [ t ] /. { u_ , v_ } { u , v , V [ u , v ] } , { t , t0 , t1 } , { x , y , V } , PlotRange All ] ;

  • A megoldások trajektóriái

A 2D fázistérben alkalmazott korábbi ábrázolásokat kiegészítjük a V függvény szintvonalaival. Jól láthatók a trajektóriák és a szintvonalak metszéspontjai.

p00=PhasePlot[Traj[t],{t,t0,t1},AxesLabel›{x,y}];

p00 = PhasePlot [ Traj [ t ] , { t , t0 , t1 } , AxesLabel { x , y } ] ;

p0=Show[pV,p00];

p0 = Show [ pV , p00 ] ;

A fenti ábra alátámasztja az ismert tényt, hogy a vizsgált fékezett oszcillátor zérómegoldása aszimptotikusan stabilis.

  • Fáziskép animáció a fázissíkon E

    E

Már tapasztaltuk, hogy a fázisképek jól mutatják valamely adott halmazból induló megoldások változását. Az ismert animációt a Ljapunov függvény szintvonalaival kombinálva hatékony vizualizációs eszközhöz jutunk.

p2=PhaseMap[Traj[t],{t,t0,t1,dt},{Hue[0,0,0]}, PlotRange›{{x1,x2},{y1,y2}},Axes›True];

p2 = PhaseMap [ Traj [ t ] , { t , t0 , t1 , dt } , { Hue [ 0 , 0 , 0 ] } , PlotRange { { x1 , x2 } , { y1 , y2 } } , Axes True ] ;

p2a=Map[Show[p0,#]///&///,p2];

p2a = Map [ Show [ p0 , # ] ///&/// , p2 ] ;

A fázisképek és trajektóriák grafikus tömbben:

p2b=Show[GraphicsArray[Partition[p2a,3]]];

p2b = Show [ GraphicsArray [ Partition [ p2a , 3 ] ] ] ;

  • 3D megjelenítések

A 3D megjelenítések során a V(x,y)=konst.

V ? ( x , y ) = konst . felületek megjelenítése nem célszerű, mivel ezek a bezárt tartományt eltakarják. Ehelyett adott ? időpillanatokban a V szintvonalait (az általános esetben a V(?,x,y)=konst. V ? ( ? , x , y ) = konst . vonalakkal) ábrázoljuk.

pV3d=PlotContourLine3D[V[x,y],{x,x1,x2},{y,y1,y2}, {t,t0,t1,dt},Contours›level];

pV3d = PlotContourLine3D [ V [ x , y ] , { x , x1 , x2 } , { y , y1 , y2 } , { t , t0 , t1 , dt } , Contours level ] ;

field3d=PlotVectorField3D[Flatten[{1,xydot}]/.eqnparm,{t,t0,t1},{x,x1,x2},{y,y1,y2},PlotPoints›7];

field3d = PlotVectorField3D [ Flatten [ { 1 , xydot } ] /. ? eqnparm , { t , t0 , t1 } , { x , x1 , x2 } , { y , y1 , y2 } , PlotPoints 7 ] ;

showsol3d=SolPlot3DBW[Traj[t],{t,t0,t1}];

showsol3d = SolPlot3DBW [ Traj [ t ] , { t , t0 , t1 } ] ;

Show[field3d,showsol3d,pV3d,Axes›True];

Show [ field3d , showsol3d , pV3d , Axes True ] ;

  • Fáziskép animáció 3D-ben E

    E

phvol=PhaseVolBW[Traj[t],{t,t0,t1,dt},{Thickness[0.015],Hue[0,0,0]}];

phvol = PhaseVolBW [ Traj [ t ] , { t , t0 , t1 , dt } , { Thickness [ 0.015 ] , Hue [ 0 , 0 , 0 ] } ] ;

Show[pV3d,phvol];

Show [ pV3d , phvol ] ;

Még egy példa

10.2.1. Példa. Egy nemautonóm oszcillátor

Tekintsük az

x''+(a(t)+a_0)x'+sin(x)=0
x '' + ( a ? ( t ) + a 0 ) ? x ' + sin ? ( x ) = 0

egyenletet, amelyben a(t)?0

a ? ( t ) ? 0 , a_0///>///0 a 0 ///>/// 0 . Írjuk át rendszer alakba:

x^'=y, y^'=–(a(t)+a_0)y–sin(x). x ' = y , y ' = - ( a ? ( t ) + a 0 ) ? y - sin ? ( x ) . (10.2.2)

Legyen a(t)=sin^2(t)

a ? ( t ) = sin 2 ( t ) és a_0=0.1 a 0 = 0.1 .

xyvar={x,y};eqnparm={a[t]›Sin[t]^2,a0›0.1};

xyvar = { x , y } ; eqnparm = { a [ t ] Sin [ t ] 2 , a0 0.1 } ;

xydot={y,–(a[t]+a0)y–Sin[x]};

xydot = { y , - ( a [ t ] + a0 ) ? y - Sin [ x ] } ;

t0=0;t1=10;dt=2; x1=–1.6; x2=1.6; y1=–1.6; y2=1.6;

t0 = 0 ; t1 = 10 ; dt = 2 ; x1 = - 1.6 ; x2 = 1.6 ; y1 = - 1.6 ; y2 = 1.6 ;

A rendszer energiafüggvénye:

W(x,x')=(1–cos(x)+(y^2/2)).
W ? ( x , x ' ) = ( 1 - cos ? ( x ) + y 2 2 ) .

W[t_,x_,y_]:=(y^2/2)+1–Cos[x];

W [ t_ , x_ , y_ ] := y 2 2 + 1 - Cos [ x ] ;

Ennek deriváltja:

DW=?_t W[t,x,y]+Simplify[Grad[W[t,x,y],xyvar].xydot]

DW = ? t W [ t , x , y ] + Simplify [ Grad [ W [ t , x , y ] , xyvar ] . xydot ]

–y^2(a0+a[t])

- y 2 ? ( a0 + a [ t ] )

Annak ellenére, hogy deriváltja nem negatív definit, a stabilitási vizsgálatokban mégis az energia-függvényt használjuk leggyakrabban. Most tekintsünk egy másik, bár kissé bonyolulabb segédfüggvényt. A kérdés, hogy hogyan lehet ilyen alakú segédfüggvényt találni, teljesen jogos. A válasz csak részben megnyugtató: előzetes formális meggondolások alapján kísérletezéssel, illetve az irodalmi leírások ismeretében (lásd például [24]).

V[t_,x_,y_]:=(((y+a0Sin[x])^2/2)+(1+a0a[t]–a0^2)(1–Cos[x]))

V [ t_ , x_ , y_ ] := ( ( y + a0 ? Sin [ x ] ) 2 2 + ( 1 + a0 ? a [ t ] - a0 2 ) ? ( 1 - Cos [ x ] ) )

Ennek deriváltja:

DV=Simplify[?_t V[t,x,y]+Grad[V[t,x,y],xyvar].xydot]

DV = Simplify [ ? t V [ t , x , y ] + Grad [ V [ t , x , y ] , xyvar ] . xydot ]

–y^2a[t]+a0(–y^2–2a0ySin[x]–Sin[x]^2+yCos[x](y+a0Sin[x])–(–1+Cos[x])a^'[t])

- y 2 ? a [ t ] + a0 ? ( - y 2 - 2 ? a0 ? y ? Sin [ x ] - Sin [ x ] 2 + y ? Cos [ x ] ? ( y + a0 ? Sin [ x ] ) - ( - 1 + Cos [ x ] ) ? a ' [ t ] )

A konkrét esetben:

DV1=Simplify[?_t(V[t,x,y]/.eqnparm)+Grad[V[t,x,y],xyvar].xydot/.eqnparm]

DV1 = Simplify [ ? t ( V [ t , x , y ] /. eqnparm ) + Grad [ V [ t , x , y ] , xyvar ] . xydot /. eqnparm ]

–0.1y^2+0.2 Cos[t]Sin[t]–1. y^2Sin[t]^2+Cos[x](–0.2Cos[t]Sin[t]+y(0.1 y+0.01 Sin[x]))–0.02 ySin[x]–0.1 Sin[x]^2

- 0.1 ? y 2 + 0.2 Cos [ t ] ? Sin [ t ] - 1. y 2 ? Sin [ t ] 2 + Cos [ x ] ? ( - 0.2 ? Cos [ t ] ? Sin [ t ] + y ? ( 0.1 y + 0.01 Sin [ x ] ) ) - 0.02 y ? Sin [ x ] - 0.1 Sin [ x ] 2

Ellenőrizhető, hogy ha a'(t)?ß///<///2

a ' ? ( t ) ? ß ///</// 2 , akkor V deriváltja negatív definit. Az egyenletet az alábbi körvonal pontjaiban oldjuk meg.

Icond=Table[?2N[{Cos[u],Sin[u]},2],{u,0,2?,(?/8)}];

Icond = Table [ 2 ? N [ { Cos [ u ] , Sin [ u ] } , 2 ] , { u , 0 , 2 ? ? , ? 8 } ] ;

A Ljapunov függvény változása a megoldások mentén:

A rendszer nemautonóm, ezért a 2D ábrázolásokra nincs lehetőségünk.

Rajzoljuk meg a megoldások és a Ljapunov függvény megoldások menti viselkedését. Csak két jellemző megoldás grafikonját láthatjuk. E

E

Az integrálgörbék a Ljapunov függvény szintvonalaival:

Végül a fázisképek változása: E

E