Ugrás a tartalomhoz

Impulzív jelenségek modelljei

Karsai János

Typotex

5.5. Általános megoldóprogram: IDERKSolve

5.5. Általános megoldóprogram: IDERKSolve

Az általános impulzív rendszerek közelítő megoldására készült az IDERKSolve programcsomag. Emlékeztetőül, az ilyen rendszer alakja az alábbi:

x'=f(t,x), ha (S_i(t,x(t))?0) (i=1,...K), x(t+0)=I_i(t,x(t–0)); ha S_i(t,x(t))=0 (i=1,...K). x ' = f ? ( t , x ) , ha ( S i ? ( t , x ? ( t ) ) ? 0 ) ? ( i = 1 , ... K ) , x ? ( t + 0 ) = I i ? ( t , x ? ( t - 0 ) ) ; ha S i ? ( t , x ? ( t ) ) = 0 ? ( i = 1 , ... K ) . (5.5.1)

Elméletileg végtelen sok S_i

S i függvény adható meg, de a megoldóprogram számára csak véges sok definiálható.

Parancsösszefoglaló

A rendszer megoldásához meg kell adni a differenciálegyenletet, az impulzusok időpontját meghatározó véges sok S_i

S i függvényt, az I_i I i impulzusfüggvényeket és a kezdeti értékek listáját. Az impulzusok bekövetkezési idejének a meghatározásához az S_i(t,x(t))=0 S i ? ( t , x ? ( t ) ) = 0 egyenlet megoldása szükséges, ezért a differenciálegyenlet megoldására az NDSolve alkalmazása nem célszerű. Két impulzus között a differenciálegyenlet megoldását egy rögzített lépésnagyságú Runge–Kutta módszert alkalmazó rutin oldja meg (Maeder [19]), és az S_i(t,x(t)) S i ? ( t , x ? ( t ) ) függvények előjelváltását lépésenként figyeljük. Ha valamelyik S_i(t,x(t)) S i ? ( t , x ? ( t ) ) függvény előjelet vált, akkor a megfelelő I_i(t,x(t)) I i ? ( t , x ? ( t ) ) impulzust alkalmazzuk. A rendszerek viselkedését, illetve a közelítő megoldás valóságtartalmát nagyban befolyásolja, hogy az impulzust az S_i(t,x(t))=0 S i ? ( t , x ? ( t ) ) = 0 feltétel bekövetkezését megelőző vagy az azt követő lépésben alkalmazzuk. Ezért ehetőség van ennek megválasztására. Emiatt a független változó nem mindig változik az ugrások során (ugyanahhoz a pillanathoz tartozik az ugrás előtti és utáni értéke is) ezért a megoldásokat nem közelítjük InterpolatingFunction objektumokkal, hanem {{t_k,x_1(t_k),...,x_n(t_k)}....} { { t k , x 1 ( t k ) , ... , x n ( t k ) } ... . } alakú listákban adjuk meg.

A program indítása az alábbi módon történik:

IDERKSolve[xdot,impulse,var,x0list,{t,t0,t1,dt}]

Az eredmény {sol1,...,solM} alakú lista, ahol soli az i-edik kezdeti értékhez tartozó megoldás, aminek szerkezete: {{t_0,x_1(t_0),...,x_n(t_0)}, ..., {t_i,x_1(t_i),...,x_n(t_i)}...

{ t 0 , x 1 ( t 0 ) , ... , x n ( t 0 ) } , ... , { t i , x 1 ( t i ) , ... , x n ( t i ) } ... }. Az egyes paraméterek jelentése az alábbi:

var var {x1,x2,..., xn}, { x1 , x2 , ... , xn } , függő változók listája
xdot xdot {f1[t,x1,...,xn],...,fn[t,x1,...,xn]}, { f1 [ t , x1 , ... , xn ] , ... , fn [ t , x1 , ... , xn ] } , a differenciálegyenlet jobb oldala
x0list x0list {{x10,...,xn0},...{x1M,...,xnM}} { { x10 , ... , xn0 } , ... { x1M , ... , xnM } } ? kezdőértékek listája, M db n dimenziós vektor
t, dt t , dt független változó; közelítő algoritmus lépésköze
[t0,t1] [ t0 , t1 ] megoldás intervalluma
Impulse Impulse az impulzusok adatainak listája {{S1,I1,jumpstep1},...,{SK,IK,jumptepsK}} { { S1 , I1 , jumpstep1 } , ... , { SK , IK , jumptepsK } } alakban
Si Si az Si=0 Si = 0 feltételben szereplő tetszőleges skalárfüggvény
Ii Ii {Ii1,...,Iin}:az Si=0 { Ii1 , ... , Iin } : az ? Si = 0 felülethez tartozó impulzus
jumpstep–i jumpstep–i i-ik impulzus technikai paramétere, értéke 0 vagy 1 lehet
  • A jumpstep paraméter értelmezése:

A formálisan is pontos megoldáshoz az S_i(t,x(t))=0

S i ? ( t , x ? ( t ) ) = 0 feltétel ellenőrzése lenne szükséges, ami egyetlen közelítő megoldó eljárás alkalmazása esetén sem lehetséges. Ehelyett a differenciálegyenlet megoldásánál az S_i(t,x(t)) S i ? ( t , x ? ( t ) ) függvények előjelváltását figyeljük. Ekkor azonban felmerül a kérdés, hogy az impulzusfüggvényt az előjelváltás előtt vagy után alkalmazzuk. Az alábbi példák során látni fogjuk, hogy mindkét esetre szükség lehet a rendszer jellegétől függően. Ezért alkalmazzuk a jumpstep paramétert minden egyes felület esetén. Jelentése a következő:

jumpstep = 0: az Ii impulzusfüggvényt az Si előjelváltása előtti utolsó lépésre alkalmazzuk. jumpstep = 1: az Ii impulzusfüggvényt az Si előjelváltása utáni első lépésre alkalmazzuk.

  • Korlátozások, esetleges problémák

Ha valamely i-re az S_i(?,x(?))=0

S i ( ? , x ? ( ? ) ) = 0 feltétel teljesül, és a megoldás során t_k///<///?///<///t_(k+1) t k ///</// ? ///</// t k + 1 (t_k t k most a Runge–Kutta eljárás során vett k-ik lépés) és S_i(t_k,x(t_k))S_i(t_(k+1),x(t_(k+1)))?0 S i ( t k , x ? ( t k ) ) ? S i ( t k + 1 , x ? ( t k + 1 ) ) ? 0 , (S_i S i nem vált előjelet), vagy S_i(t_k,x(t_k))=0 S i ( t k , x ? ( t k ) ) = 0 , akkor a program nem feltétlenül működik megfelelően. Ez a korlátozás összhangban van azzal a jelenséggel, hogy ha a rendszer differenciálegyenlet részének egy megoldása valamely intervallumon kielégíti az S_i(t,x(t))=0 S i ( t , x ? ( t ) ) = 0 feltételt, akkor ezekben a (t,x(t)) ( t , x ? ( t ) ) pontokban a kezdetiérték-problémának nincs megoldása. Valódi korlátozást csak az az eset jelent, amikor a differenciálegyenlet x(t) x ? ( t ) megoldásának grafikonja érinti az S_i(t,x)=0 S i ( t , x ) = 0 felületet és az érintési pontban S_i(t,x(t)) S i ( t , x ? ( t ) ) nem vált előjelet.

Bizonyos esetekben előfordulhat, hogy a jumpstep paraméter beállítását kezdeti értékekhez kell igazítani.

A program használata

A program használatát példákon keresztül mutatjuk be. A rendszerek megoldásait most elsősorban technikai szempontból vizsgáljuk, matematikai problémákkal később foglalkozunk.

5.5.1. Példa. Falba ütköző rugó mozgása

Tekintsük a rugó mozgását leíró harmonikus oszcillátor egyenletét:

x'=y, y'=–x
x ' = y , y ' = - x

var={x,y}; xdot={y,–x};

var = { x , y } ; xdot = { y , - x } ;

A mozgás fizikailag korlátozva van. Ha a kitérés elér egy értéket, mondjuk x=0.5,

x = 0.5 , akkor rugalmatlan ütközés következik be. Tehát ha

x(t–0)–0.5=0, akkor y(t+0)=–0.9y(t–0).
x ? ( t - 0 ) - 0.5 = 0 , akkor y ? ( t + 0 ) = - 0.9 ? y ? ( t - 0 ) .

Először oldjuk meg a rendszert úgy, hogy az ugrás a feltétel teljesülése előtti lépésben következzen be (jumpstep=0). Az impulzus megadása ekkor az alábbi alakú:

Impulse={{x–0.5,{x,–0.9y},0}};

Impulse = { { x - 0.5 , { x , - 0.9 ? y } , 0 } } ;

A rendszer független változója t, és a [t0,t1]

[ t0 , t1 ] intervallumon oldjuk meg az {0,1} és {1,0} kezdeti feltételekkel, amelyeket az x0list változóban tárolunk.

t0=0; t1=6; dt=0.05; x0list={{0,1},{1,0}};

t0 = 0 ; t1 = 6 ; dt = 0.05 ; x0list = { { 0 , 1 } , { 1 , 0 } } ;

Az IDERKSolve program most már indítható:

sol=IDERKSolve[xdot,Impulse,var,x0list,{t,t0,t1,dt}];

sol = IDERKSolve [ xdot , Impulse , var , x0list , { t , t0 , t1 , dt } ] ;

Az első kezdeti értékhez tartozó megoldás sol[[1]], amelynek a kezdeti pillanathoz tartozó eleme sol[[1,1]] (={0.,1.,0.}).

A második kezdeti értékhez tartozó megoldás pedig sol[[2]]. Az egyes koordináták (x(t)

x ? ( t ) és y(t) y ? ( t ) ) mint pontsorozattal megadott görbék előállításához az egyszerű Coordinate függvényt használjuk.

Coordinate[data_,n_]:=({#1[[1]],#1[[n]]}///&///)/@data

Coordinate [ data_ , n_ ] := ( { #1 [ 1 ] , #1 [ n ] } ///&/// ) /@ data

Coordinate[data_,{n_,m_}]:=({#1[[n]],#1[[m]]}///&///)/@data

Coordinate [ data_ , { n_ , m_ } ] := ( { #1 [ n ] , #1 [ m ] } ///&/// ) /@ data

Ezután az a megoldásokat a ListPlot, a MultipleListPlot vagy az egyszerű

ListPlot1[data_List,opt___]:=Show[(ListPlot[#1,DisplayFunction›Identity,opt]///&///)/@data,DisplayFunction›$DisplayFunction]

ListPlot1 [ data_List , opt___ ] := Show [ ( ListPlot [ #1 , DisplayFunction Identity , opt ] ///&/// ) /@ data , DisplayFunction $DisplayFunction ]

utasításokkal ábrázolhatjuk. Az első kezdeti értékhez tartozó megoldás komponensei:

ListPlot1[{Coordinate[sol[[1]],2],Coordinate[sol[[1]],3]},PlotJoined›True,GridLines›Automatic,PlotRange›All];

ListPlot1 [ { Coordinate [ sol [ 1 ] , 2 ] , Coordinate [ sol [ 1 ] , 3 ] } , PlotJoined True , GridLines Automatic , PlotRange All ] ;

Nyilvánvaló, hogy a folytonos vonal az x(t)

x ? ( t ) , szakadásos görbe (a számítógépes ábrázolásnál a valódi szakadás helyett függőleges vonal jelenik meg) pedig az y(t) y ? ( t ) koordináta. A második kezdeti értékhez tartozó megoldás komponensei:

ListPlot1[{Coordinate[sol[[2]],2],Coordinate[sol[[2]],3]},PlotJoined›True,GridLines›Automatic,PlotRange›All];

ListPlot1 [ { Coordinate [ sol [ 2 ] , 2 ] , Coordinate [ sol [ 2 ] , 3 ] } , PlotJoined True , GridLines Automatic , PlotRange All ] ;

Megjegyezzük, hogy az első esetben a kitérés nem lehet 0.5-nél nagyobb, a másodikban pedig nem lehet 0.5-nél kisebb.

5.5.2. Példa. Módosított példa: az ugrás a feltételt követően történik

Minden, az előző példában megadott definíció változatlan marad, kivéve az Impulzus változó harmadik elemét:

var={x,y};xdot={y,–x}; Impulse={{x–0.5,{x,–0.9y},1}};

var = { x , y } ; xdot = { y , - x } ; Impulse = { { x - 0.5 , { x , - 0.9 ? y } , 1 } } ;

t0=0; t1=6; dt=0.05; x0list={{0,1},{1,0}};

t0 = 0 ; t1 = 6 ; dt = 0.05 ; x0list = { { 0 , 1 } , { 1 , 0 } } ;

Oldjuk meg a rendszert és ábrázoljuk az első megoldást.

sol=IDERKSolve[xdot,Impulse,var,x0list,{t,t0,t1,dt}];

sol = IDERKSolve [ xdot , Impulse , var , x0list , { t , t0 , t1 , dt } ] ;

ListPlot1[{Coordinate[sol[[1]],2],Coordinate[sol[[1]],3]},PlotJoined›True,GridLines›Automatic,PlotRange›All];

ListPlot1 [ { Coordinate [ sol [ 1 ] , 2 ] , Coordinate [ sol [ 1 ] , 3 ] } , PlotJoined True , GridLines Automatic , PlotRange All ] ;

Látható, hogy éppen a korlátozó felületről való (általunk elvárt) visszapattanás nem következik be. Ehelyett, az x=0.5

x = 0.5 érték elérésekor a megoldó program bizonytalanná válik, és a megoldás „pattogni” kezd.

5.5.3. Példa. Újabb módosítás: két ütközési felület

Most tegyük fel, hogy a rugó az x(t)=0.5

x ? ( t ) = 0.5 felületről visszapattan és az x(t)=0 x ? ( t ) = 0 kitérésnél egy pillanatra fékező felülettel találkozik:

x(t–0)–0.5=0, akkor y(t+0)=–ay(t–0), x(t–0)=0, akkor y(t+0)=by(t–0).
x ? ( t - 0 ) - 0.5 = 0 , akkor y ? ( t + 0 ) = - a ? y ? ( t - 0 ) , x ? ( t - 0 ) = 0 , akkor y ? ( t + 0 ) = b ? y ? ( t - 0 ) .

var={x,y}; xdot={y,–x};

var = { x , y } ; xdot = { y , - x } ;

Impulse={{x–0.5,{x,–0.9y},0},{x,{x,0.8y},1}};

Impulse = { { x - 0.5 , { x , - 0.9 ? y } , 0 } , { x , { x , 0.8 y } , 1 } } ;

A fékező, de visszapattanást nem okozó impulzus esetén a jumpstep paraméter értékét 1-re állítottuk.

t0=0; t1=6; dt=0.01; x0list={{0.4,1}};

t0 = 0 ; t1 = 6 ; dt = 0.01 ; x0list = { { 0.4 , 1 } } ;

Oldjuk meg a rendszert és ábrázoljuk az első megoldást.

sol=IDERKSolve[xdot,Impulse,var,x0list,{t,t0,t1,dt}];

sol = IDERKSolve [ xdot , Impulse , var , x0list , { t , t0 , t1 , dt } ] ;

ListPlot1[{Coordinate[sol[[1]],2],Coordinate[sol[[1]],3]},PlotJoined›True,GridLines›Automatic,PlotRange›All];

ListPlot1 [ { Coordinate [ sol [ 1 ] , 2 ] , Coordinate [ sol [ 1 ] , 3 ] } , PlotJoined True , GridLines Automatic , PlotRange All ] ;

5.5.4. Példa. Impulzusok rögzített pillanatokban

Az IDERKSolve program alkalmas az IDESolve program által hatékonyan kezelhető rögzített pillanatokban ható, általában pedig tetszőleges vegyes típusú impulzusokkal perturbált rendszerek megoldására is. Tekintsük a már az IDESolve program bemutatásánál megismert lineáris rendszert:

x'=–x, t?i, x(n+0)=x(n–0)+1.
x ' = - x , t ? i , x ? ( n + 0 ) = x ? ( n - 0 ) + 1 .

Az IDERKSolve program számára a rendszer a következőképpen írható le:

var={x}; xdot={–x}; Impulse={{Sin[?t],{x+1},1}};

var = { x } ; xdot = { - x } ; Impulse = { { Sin [ ? ? t ] , { x + 1 } , 1 } } ;

A jumpstep paraméter értékét 1-re állítottuk.

t0=0; t1=6; dt=0.01; x0list={{2.}};

t0 = 0 ; t1 = 6 ; dt = 0.01 ; x0list = { { 2. } } ;

Oldjuk meg a rendszert és ábrázoljuk az első megoldást.

sol=IDERKSolve[xdot,Impulse,var,x0list,{t,t0,t1,dt}];

sol = IDERKSolve [ xdot , Impulse , var , x0list , { t , t0 , t1 , dt } ] ;

ListPlot1[{Coordinate[sol[[1]],2]},PlotJoined›True,GridLines›Automatic,PlotRange›{0,3}];

ListPlot1 [ { Coordinate [ sol [ 1 ] , 2 ] } , PlotJoined True , GridLines Automatic , PlotRange { 0 , 3 } ] ;

Javasoljuk az olvasónak, próbálja ki, hogy mi történik, ha a jumpstep változó értékét 0-ra állítjuk.