Ugrás a tartalomhoz

Impulzív jelenségek modelljei

Karsai János

Typotex

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

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

Az R^n

R n tér minden x_0 x 0 pontjában megrajzolva az f(x_0) f ? ( x 0 ) (f: R^n›R^n) ( 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 G: R^n›R 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 x'=sin(t)x

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

field1D={1,xSin[t]};

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

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

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 x'=y, y'=–x–y

x ' = y , y ' = - x - y differenciálegyenletét. A trajektóriák a síkon az {y,–x–y} { y , - x - y } vektormezőhöz illeszkednek.

field2D={y,–x–y};

field2D = { y , - x - y } ;

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

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

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

Az x'=y,

x ' = y , y'=–x–y y ' = - x - y differenciálegyenlet megoldásainak grafikonjai a 3D térben az {1,y,–x–y} { 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).

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

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

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

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 {1,f(t,x)}

{ 1 , f ? ( t , x ) } iránymező 4D ábrázolást tenne szükségessé. Tekintsük az x'=y, y'=–x–y, z'=y–z x ' = y , y ' = - x - y , z ' = y - z egyenletrendszert. A megoldások trajektóriái az {y, –x–y, y–z} { y , - x - y , y - z } vektormezőhöz illeszkednek.

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

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

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

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

eqn=x^('')[t]+ax^'[t]+x[t]==0;

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

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

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]}}

{ { 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 a=1

a = 1 esetben is komplex alakot kapunk:

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

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

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

{ { 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};

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

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

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])}}

{ { 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};

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

Join[eqn,ic]

Join [ eqn , ic ]

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

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

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

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]))}

{ 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};

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

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

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

Join[eqn2d,ic2d]

Join [ eqn2d , ic2d ]

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

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

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

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)}}

{ { 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};

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

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

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

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

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

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

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

  • A 2D rendszer alak

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

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

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

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

Join[eqn2d,ic2d]

Join [ eqn2d , ic2d ]

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

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

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

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

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

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

  • A megoldás trajektóriája

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

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]

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

  • A megoldás koordinátái

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

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

  • A koordináták deriváltjai

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

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

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

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

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

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

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

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.