Impulzív jelenségek modelljei|Digitális Tankönyvtár

Főoldal > Tankönyvtár > Könyvek > Természettudományok > Matematika

Impulzív jelenségek modelljei

Karsai János

Typotex

Beágyazás

2.2. A Mathematica beépített függvényei

Ebben a pontban a differenciálegyenletek vizsgálatát segítő Mathematica eszközöket tekintjük át. Megjegyezzük, hogy több könyv és programcsomag ismert, amelyek kiterjesztik vagy elegánsabbá teszik az alapfunkciókat. Például a Gray, Mezzino, Pinsky [11] könyvéhez mellékelt programcsomag számos új formális megoldó-algoritmust tartalmaz, a VisualDSolve [26] csomagnak pedig igen elegáns megjelenítési lehetőségei vannak.

Parancsösszefoglaló

`PlotVectorField`		2D vektormező ábrázolása, Graphics`PlotField` csomag
`ListPlotVectorField`		Listával megadott 2D vektormező, Graphics`PlotField` csomag
`PlotVectorField3D`		3D vektormező ábrázolása, Graphics`PlotField3D` csomag
`ListPlotVectorField3D`		Listával megadott 3D vektormező, Graphics`PlotField3D` csomag
`DSolve`		ODE formális megoldása
`NDSolve`		ODE numerikus megoldása

Vektormezők

[D]

R n tér minden

[D]

x 0 pontjában megrajzolva az

[D]

f ? ( x 0 )

[D]

( f : R n › R n ) vektort, vektormezőt kapunk. Vektormezők valamely paraméteres görbecsalád érintővektorai (például differenciálegyenlet iránymezője) vagy valamely

[D]

G : R n › R függvény szintfelületeire ortogonális vektorok mezője. A Mathematica ezek mellett számos formálisan és adatokkal megadott vektormező ábrázolását teszi lehetővé a 2D és 3D terekben. Itt most csak a differenciálegyenletekkel kapcsolatos vektormezők (2.1. fejezetet) ábrázolásának legegyszerűbb módjait, a PlotVectorField és PlotVectorField3D utasítások használatát tekintjük át. Részletes információ a Graphics`PlotField` és a Graphics`PlotField3D` programcsomagok leírásában található.

1D differenciálegyenlet iránymezője

Tekintsük az

[D]

x ' = sin ? ( t ) ? x differenciálegyenletet. Ennek az egyenletnek az iránymezőjét az alábbi utasítással rajzolhatjuk meg.

[D]

field1D = { 1 , x ? Sin [ t ] } ;

[D]

PlotVectorField [ field1D , { t , 0 , 4 } , { x , 0 , 1.5 } , Axes › True , AxesLabel › { t , x } ] ;

2D autonóm differenciálegyenlet trajektóriáihoz tartozó vektormező

Tekintsük a csillapított harmonikus rezgés

[D]

x ' = y , y ' = - x - y differenciálegyenletét. A trajektóriák a síkon az

[D]

{ y , - x - y } vektormezőhöz illeszkednek.

[D]

field2D = { y , - x - y } ;

[D]

PlotVectorField [ field2D , { x , - 2 , 2 } , { y , - 2 , 2 } , Axes -///>/// True , AxesLabel -///>/// { x , y } ] ;

2D differenciálegyenlet iránymezője 3D-ben

[D]

x ' = y ,

[D]

y ' = - x - y differenciálegyenlet megoldásainak grafikonjai a 3D térben az

[D]

{ 1 , y , - x - y } iránymezőhöz illeszkednek.

Vigyázat! A perspektivikus ábrázolás becsaphatja szemünket. Érdemes több nézőpontból megnézni az ábrát (lásd a ViewPoint opciót).

[D]

field2D3D = { 1 , y , - x - y } ;

[D]

PlotVectorField3D [ field2D3D , { t , 0 , 2 ? ? } , { x , - 1 , 1 } , { y , - 1 , 1 } , Axes › True , AxesLabel › { t , x , y } ] ;

3D differenciálegyenletekhez tartozó vektormezők

A 3D autonóm egyenletek esetén a trajektóriák vektormezője (hasonlóan a 2D rendszerekhez) ábrázolható. Az

[D]

{ 1 , f ? ( t , x ) } iránymező 4D ábrázolást tenne szükségessé. Tekintsük az

[D]

x ' = y , y ' = - x - y , z ' = y - z egyenletrendszert. A megoldások trajektóriái az

[D]

{ y , - x - y , y - z } vektormezőhöz illeszkednek.

[D]

field3D = { y , - x - y , y - z } ;

[D]

PlotVectorField3D [ field3D , { x , - 2 , 2 } , { y , - 2 , 2 } , { z , - 2 , 2 } , Axes › True , AxesLabel › { x , y , z } ] ;

Differenciálegyenletek szimbolikus megoldása

A Mathematica DSolve utasítása bizonyos integrálható típusú közönséges differenciálegyenletek és parciális differenciálegyenletek kezdeti- és peremérték problémáinak formális megoldásait határozza meg. Sok egyszerűsítési és transzformációs lehetőséget a DSolve nem ismer. Az érdeklődő olvasónak javasoljuk Gray, Mezzino, Pinsky [11] könyvét. Az ott használt ODE programcsomag számos egyenlettípus integrálását el tudja végezni.

Általános megoldás

Példaként tekintsük a fékezett harmonikus oszcillátor egyenletét (amelynek 2D vektormezőjét és 3D iránymezőjét már ábrázoltuk) mint másodfokú egyenletet és mint síkbeli rendszert is.

A másodfokú egyenlet alak

[D]

eqn = x '' [ t ] + a ? x ' [ t ] + x [ t ] == 0 ;

[D]

sol = DSolve [ eqn , x [ t ] , t ]

${{x[t]›e^((1/2)(–a–?(–4+a^2))t)C[1]+e^((1/2)(–a+?(–4+a^2))t)C[2]}}$

[D]

{ { x [ t ] › ? 1 2 ? ( - a - - 4 + a 2 ) ? t ? C [ 1 ] + ? 1 2 ? ( - a + - 4 + a 2 ) ? t ? C [ 2 ] } }

A megoldást helyettesítési szabály alakjában kapjuk. Az általános megoldásban szereplő konstansok C[1] és C[2] (n-edfokú vagy n-dimenziós esetben, C[1], C[2],..., C[n]). Vegyük észre, hogy mivel nem tudjuk az a együttható nagyságát, az exponenciális kifejezés lehet valós és komplex is. Az

[D]

a = 1 esetben is komplex alakot kapunk:

[D]

sol1 = DSolve [ eqn /. a › 1 , x [ t ] , t ]

${{x[t]›e^(–(–1)^(1/3)t)C[1]+e^((–1)^(2/3)t)C[2]}}$

[D]

{ { x [ t ] › ? - ( - 1 ) 1 / 3 ? t ? C [ 1 ] + ? ( - 1 ) 2 / 3 ? t ? C [ 2 ] } }

A 2D rendszer alak

$eqn2d={x^'[t]==y[t],y^'[t]+y[t]+x[t]==0};$

[D]

eqn2d = { x ' [ t ] == y [ t ] , y ' [ t ] + y [ t ] + x [ t ] == 0 } ;

[D]

sol2d = Simplify [ DSolve [ eqn2d , { x [ t ] , y [ t ] } , t ] ]

${{x[t]›–(1/6)ie^(–(–1)^(1/3)t)((3i–?3+(3i+?3)e^(i?3t))C[1]+2?3(–1+e^(i?3t))C[2]), y[t]›(1/6)ie^(–(–1)^(1/3)t)(2?3(–1+e^(i?3t))C[1]+(–3i–?3+(–3i+?3)e^(i?3t))C[2])}}$

[D]

{ { x [ t ] › - 1 6 ? ? ? ? - ( - 1 ) 1 / 3 ? t ? ( ( 3 ? ? - 3 + ( 3 ? ? + 3 ) ? ? ? ? 3 ? t ) ? C [ 1 ] + 2 ? 3 ? ( - 1 + ? ? ? 3 ? t ) ? C [ 2 ] ) , y [ t ] › 1 6 ? ? ? ? - ( - 1 ) 1 / 3 ? t ? ( 2 ? 3 ? ( - 1 + ? ? ? 3 ? t ) ? C [ 1 ] + ( - 3 ? ? - 3 + ( - 3 ? ? + 3 ) ? ? ? ? 3 ? t ) ? C [ 2 ] ) } }

A megoldás most mindkét koordinátát tartalmazza. Az első komponens formálisan nem azonos a másodrendű egyenlet megoldásával. Ellenőrizhető, hogy ugyanazt a megoldáshalmazt állítják elő.

Kezdetiérték-probléma megoldása

A kezdeti értékeket az egyenlet részeként kell a DSolve utasításnak átadni. Tekintsük ismét az oszcillátor példáját mindkét esetben. Konkrét értékeket adunk meg, de a megoldást szimbolikus értékek esetén is megkapjuk.

A másodfokú egyenlet alak

$eqn={x^('')[t]+x^'[t]+x[t]==0}; ic={x[0]==1,x^'[0]==0};$

[D]

eqn = { x '' [ t ] + x ' [ t ] + x [ t ] == 0 } ; ic = { x [ 0 ] == 1 , x ' [ 0 ] == 0 } ;

[D]

Join [ eqn , ic ]

${x[t]+x^'[t]+x^('')[t]==0,x[0]==1,x^'[0]==0}$

[D]

{ x [ t ] + x ' [ t ] + x '' [ t ] == 0 , x [ 0 ] == 1 , x ' [ 0 ] == 0 }

[D]

sol = Simplify [ ExpToTrig [ x [ t ] /. ? DSolve [ Join [ eqn , ic ] , x [ t ] , t ] ] , Trig › True ]

${(1/–3i+?3)((–i+?3–2iCos[?3t]+2Sin[?3t])(Cosh[(–1)^(1/3)t]–Sinh[(–1)^(1/3)t]))}$

[D]

{ 1 - 3 ? ? + 3 ? ( ( - ? + 3 - 2 ? ? ? Cos [ 3 ? t ] + 2 ? Sin [ 3 ? t ] ) ? ( Cosh [ ( - 1 ) 1 / 3 ? t ] - Sinh [ ( - 1 ) 1 / 3 ? t ] ) ) }

A 2D rendszer alak

$eqn2d={x^'[t]==y[t],y^'[t]+y[t]+x[t]==0};$

[D]

eqn2d = { x ' [ t ] == y [ t ] , y ' [ t ] + y [ t ] + x [ t ] == 0 } ;

[D]

ic2d = { x [ 0 ] == 1 , y [ 0 ] == 0 } ;

[D]

Join [ eqn2d , ic2d ]

${x^'[t]==y[t],x[t]+y[t]+y^'[t]==0,x[0]==1,y[0]==0}$

[D]

{ x ' [ t ] == y [ t ] , x [ t ] + y [ t ] + y ' [ t ] == 0 , x [ 0 ] == 1 , y [ 0 ] == 0 }

[D]

sol2d = Simplify [ ExpToTrig [ DSolve [ Join [ eqn2d , ic2d ] , { x [ t ] , y [ t ] } , t ] ] ]

${{x[t]›(ie^(–(–1)^(1/3)t)(3i+?3+2?3e^(i?3t))/6(–1+(–1)^(1/3))),y[t]›(ie^(–(1/2)(1+i?3)t)(–1+e^(i?3t))/?3)}}$

[D]

{ { x [ t ] › ? ? ? - ( - 1 ) 1 / 3 ? t ? ( 3 ? ? + 3 + 2 ? 3 ? ? ? ? 3 ? t ) 6 ? ( - 1 + ( - 1 ) 1 / 3 ) , y [ t ] › ? ? ? - 1 2 ? ( 1 + ? ? 3 ) ? t ? ( - 1 + ? ? ? 3 ? t ) 3 } }

Kezdetiérték-probléma numerikus megoldása

Közelítő megoldásokat az NDSolve program segítségével kapunk. Használata hasonló a DSolve programhoz. A megoldáshoz meg kell adni az intervallumot, amelyen a megoldást elő kívánjuk állítani. Több, a pontosság javítására szolgáló opciója van. Részletes leírás és az alkalmazott numerikus eljárás implementációjával kapcsolatos megjegyzések a Help-ben találhatók. A megoldásokat helyettesítési szabályként kapjuk, amelynek jobb oldala InterpolatingFunction objektum, amellyel kapcsolatban néhány példát később bemutatunk. Tekintsük ismét a fékezett oszcillátort mint másodfokú egyenletet és mint 2D rendszert is.

A másodfokú egyenlet alak

$eqn={x^('')[t]+x^'[t]+x[t]==0};$

[D]

eqn = { x '' [ t ] + x ' [ t ] + x [ t ] == 0 } ;

$ic={x[0]==1.,x^'[0]==0.}; t0=0; T=6;$

[D]

ic = { x [ 0 ] == 1. , x ' [ 0 ] == 0. } ; t0 = 0 ; T = 6 ;

[D]

NDSolve [ Join [ eqn , ic ] , x [ t ] , { t , t0 , T } ]

[D]

{ { x [ t ] › InterpolatingFunction [ { { 0. , 6. } } , ///<//////>/// ] [ t ] } }

A 2D rendszer alak

$eqn2d={x^'[t]==y[t],y^'[t]+y[t]+x[t]==0};$

[D]

eqn2d = { x ' [ t ] == y [ t ] , y ' [ t ] + y [ t ] + x [ t ] == 0 } ;

[D]

ic2d = { x [ 0 ] == 1. , y [ 0 ] == 0. } ; t0 = 0 ; T = 6 ;

[D]

Join [ eqn2d , ic2d ]

${x^'[t]==y[t],x[t]+y[t]+y^'[t]==0,x[0]==1.,y[0]==0.}$

[D]

{ x ' [ t ] == y [ t ] , x [ t ] + y [ t ] + y ' [ t ] == 0 , x [ 0 ] == 1. , y [ 0 ] == 0. }

[D]

nsol2d = NDSolve [ Join [ eqn2d , ic2d ] , { x [ t ] , y [ t ] } , { t , t0 , T } ]

[D]

{ { x [ t ] › InterpolatingFunction [ { { 0. , 6. } } , ///<//////>/// ] [ t ] , y [ t ] › InterpolatingFunction [ { { 0. , 6. } } , ///<//////>/// ] [ t ] } }

A megoldás trajektóriája

[D]

ParametricPlot [ Evaluate [ oscsol [ t ] ] , { t , 0 , 6 } ] ;

Az „InterpolatingFunction” objektum

Most a misztikus InterpolatingFunction objektummal foglalkozunk. Az adatok függvényekkel való közelítésének egyik alapvető módja az interpoláció, amikor bizonyos helyeken előírjuk a függvény és esetleg valamely deriváltjainak az értékét. A Mathematica több interpolációs utasítással rendekezik (lásd az Interpolation és a hozzá kapcsolódó utasításokat). Mindegyikük InterpolatingFunction objektumot eredményez. Mivel a differenciálegyenletek megoldásainak közelítő algoritmusai adathalmazokat állítanak elő, kézenfekvő az interpoláció azonnali elvégzése. Az InterpolatingFunction objektum bonyolult szerkezetű, de szerencsére ismeretére nincs gyakran szükség (a kivételek közé tartozik az impulzív rendszerek megoldására később bemutatandó IDESolve programcsomag). Használata nagyon kényelmes, ugyanúgy kezelhető, ábrázolható, deriváltható(!), mint a formulákkal megadott függvények. Példaként tekintsük most a fékezett oszcillátor imént kapott megoldásait.

$oscsol[t_]=First[{x[t],y[t]}/.nsol2d]$

[D]

oscsol [ t_ ] = First [ { x [ t ] , y [ t ] } /. nsol2d ]

A megoldás koordinátái

[D]

Plot [ Evaluate [ oscsol [ t ] ] , { t , 0 , 6 } ] ;

A koordináták deriváltjai

$Plot[Evaluate[oscsol^'[t]],{t,0,6}];$

[D]

Plot [ Evaluate [ oscsol ' [ t ] ] , { t , 0 , 6 } ] ;

$Plot[Evaluate[oscsol^('')[t]],{t,0,6}];$

[D]

Plot [ Evaluate [ oscsol '' [ t ] ] , { t , 0 , 6 } ] ;

$Plot[Evaluate[oscsol^((3))[t]],{t,0,6}];$

[D]

Plot [ Evaluate [ oscsol ( 3 ) [ t ] ] , { t , 0 , 6 } ] ;

$Plot[Evaluate[oscsol^((4))[t]],{t,0,6}];$

[D]

Plot [ Evaluate [ oscsol ( 4 ) [ t ] ] , { t , 0 , 6 } ] ;

Vegyük észre, hogy a második deriváltak még teljesen pontosak, a harmadik derivált kicsit már tüskés. A negyedik derivált viszont már használhatatlan.