Harjoitus 9 -- Ratkaisut

1

DSolve ja NDSolve antavat vastaukset sääntölistoina. Sääntöjä sovelletaan /.-operaattorilla. Kannattaa nimetä ratkaisu, tai poimia suoraan ratkaisufunktio. Flatten-funktio poistaa ylimääräiset aaltosulkeet. Sitä tarvitaan ainakin ParametricPlot-funktion yhteydessä.

q=Flatten[DSolve[{x1'[t]+2t x1[t]==Sin[t],x1[0]==0},x1[t],t]]

Plot[x1[t]/.q,{t,0,10}]

NDSolve toimii hyvin epälineaaristen diff. yhtälöiden kanssa. Se antaa vastauksen InterpolatingFunction-oliona (vrt. harj. 4). Jos diff. yhtälössä esiintyy n:nnen kertaluvun derivaattoja, NDSolve vaatii n kpl reunaehtoja (jotta kaikki yleisen ratkaisun vakiot saisivat numeerisen arvon).

T=3.5;
y[t_]=x2[t]/.
Flatten[
NDSolve[{x2'''[t]-((x2'[t])/(2+Cos[10((x2[t])^2)]))+1==0,
x2[0]==1,x2'[0]==1,x2''[0]==1},
x2[t],{t,-T,T}]
]
Plot[y[t],{t,-T,T}]

2

Yleisimmässä tapauksessa ratkaistava yhtälö on muotoa (merk. (Subscript[ω, d])->η)

eq=x''[t]+γ x'[t]+(ω^2)x[t]==f0 Cos[η t];

jonka yleinen ratkaisu on

sol=Flatten[DSolve[eq,x[t],t]]

Vakiot C[1] ja C[2] määräytyvät myöhemmin alkuehdoista. Annetaan parametreille γ, ω, f0 ja η arvot sääntöjen avulla.

a

y[t_]=x[t]/.Flatten[DSolve[{eq,x[0]==1,x'[0]==0},x[t],t]]
val={f0->0,γ->0,η->0,ω->1};
g1=Plot[y[t]/.val,{t,0,4π}];
g2=ParametricPlot[{y[t],y'[t]}/.val,{t,0,4π}];
GraphicsArray[{g1,g2}]

ParametricPlot-funktion piirtämässä kuvaajassa x-akseli on poikkeama tasapainoasemasta ja y-akseli on nopeus. Ympyrän säteen neliö on verrannollinen värähtelijän energiaan (pysyy tässä vakiona).

b

Sijoitetaan parametrien arvot jo ratkaisuvaiheessa.

val={f0->0,γ->0.1,η->0,ω->1};
y[t_]=x[t]/.Flatten[DSolve[{eq/.val,x[0]==1,x'[0]==1},x[t],t]];
g1=Plot[y[t],{t,0,8π}];
g2=ParametricPlot[{y[t],y'[t]},{t,0,8π}];
GraphicsArray[{g1,g2}]

Nyt värähtely vaimenee ajan funktiona kitkavoiman takia, joten (x[t],x'[t])-tasossa kuvattuna systeemi lähestyy spiraalimaisesti origoa, jossa siis kappale on paikoillaan tasapainoasemassa. Huom. kappaleen rata (x[t],x'[t])-koordinaatistossa ei ole nyt "ympyräspiraali", vaan elliptinen, mikä johtuu alkuehdoista.

c

T=80;
val={f0->1,γ->0.1,η->1,ω->1};
y[t_]=x[t]/.Flatten[DSolve[{eq/.val,x[0]==1,x'[0]==0},x[t],t]];
g1=Plot[y[t],{t,0,T}];
g2=ParametricPlot[{y[t],y'[t]},{t,0,T}];
GraphicsArray[{g1,g2}]

Ajava voima saa kappaleen värähtelemään. Värähtely voimistuu kunnes kitkavoima kasvaa tarpeeksi suureksi, että kasvu pysähtyy. Jos kitkaa ei olisi, kasvu jatkuisi loputtomiin. Faasiavaruudessa (eli (x[t],x'[t])-koordinaatistossa) tarkasteltuna kappaleen liike lähestyy ympyrärataa.

Värähtelyjen voimistumista kuvaa teho f(t)x'(t) (voima kertaa nopeus).

Plot[y'[t]f0 Cos[η t]/.val,{t,0,T}]

Keskimääräinen teho saadaan yo. lausekkeen aikakeskiarvona. Seuraavan Plot-komennon suorittaminen saattaa kestää kauan. Huom. taajuuden η arvon sijoittaminen (/.val) täytyy tehdä joka kohdassa erikseen.

P[t0_]:=((η/(2π))/.val)NIntegrate[y'[t]f0 Cos[η t]/.val,{t,t0,t0+((2π)/η)/.val}]

Plot[P[t0],{t0,0,T}]

d

T=80;
val={f0->1,γ->0.1,η->0.9,ω->1};
y[t_]=x[t]/.Flatten[DSolve[{eq/.val,x[0]==1,x'[0]==0},x[t],t]];
g1=Plot[y[t],{t,0,T}];
g2=ParametricPlot[{y[t],y'[t]},{t,0,T}];
GraphicsArray[{g1,g2}]

ω on systeemin "luonnollinen taajuus", joka riippuu esim. jousivakiosta ja värähtelevän kappaleen massasta. Edellisessä kohdassa ajavan voiman taajuus oli sama kuin luonnollinen taajuus, jolloin värähtelyjen kasvu oli hyvin voimakasta (resonanssi). Nyt ajavan voiman taajuus poikkeaa hieman luonnollisesta taajuudesta, ja värähtelyjen saavuttama amplitudi on selvästi pienempi.

g1=Plot[y'[t]f0 Cos[η t]/.val,{t,0,T}];
P[t0_]:=((η/(2π))/.val)NIntegrate[y'[t]f0 Cos[η t]/.val,{t,t0,t0+((2π)/η)/.val}]
g2=Plot[P[t0],{t0,0,T}];
GraphicsArray[{g1,g2}]

Nyt teho f(t)x'(t) menee välillä negatiiviseksi, mikä tarkoittaa sitä, että ajava voima jarruttaa värähtelevää kappaletta.

d-ekstra

Tutkitaan tilannetta, jossa ajavan voiman taajuus ja systeemin luonnollinen taajuus poikkeavat vielä enemmän toisistaan.

T=50;
val={f0->2.0,γ->0.1,η->10.0,ω->1.0};
y[t_]=x[t]/.Flatten[DSolve[{eq/.val,x[0]==1,x'[0]==0},x[t],t]];
g1=Plot[y[t],{t,0,T}];
g2=ParametricPlot[{y[t],y'[t]},{t,0,T}];
GraphicsArray[{g1,g2}]

g1=Plot[y'[t]f0 Cos[η t]/.val,{t,0,T}];
P[t0_]:=((η/(2π))/.val)NIntegrate[y'[t]f0 Cos[η t]/.val,{t,t0,t0+((2π)/η)/.val}]
g2=Plot[P[t0],{t0,0,T}];
GraphicsArray[{g1,g2}]

Saavutettava värähtelyjen amplitudi on hyvin pieni.

e

T=50;
val={f0->1.0,γ->0.9,η->1.5,ω->1.0};
y[t_]=x[t]/.Flatten[DSolve[{eq/.val,x[0]==0,x'[0]==0},x[t],t]];
g1=Plot[y[t],{t,0,T}];
g2=ParametricPlot[{y[t],y'[t]},{t,0,T}];
GraphicsArray[{g1,g2}]

Ajava voima saa aikaa värähtelyt. Amplitudi ei ole kovin suuri, sillä ajava voima ei ole resonanssissa (η≠ω).

3

Ratkaistavat yhtälöt (epälineaarisia, joten käytetään funktiota NDSolve)

eq1=N1'[t]==ϵ1 N1[t]-β1 ((N1[t])^2)-γ1 N1[t]N2[t];
eq2=N2'[t]==ϵ2 N2[t]-β2 ((N2[t])^2)-γ2 N2[t]N1[t];

Analyyttista ratkaisuakin voi yrittää, mutta ratkaisu saattaa kestää todella kauan.

DSolve[{eq1,eq2},{N1[t],N2[t]},t]

a

T=20;
val={ϵ1->0.3,ϵ2->0.4,β1->0.0,β2->0.0,γ1->0.0,γ2->0.0};
sol=NDSolve[{eq1,eq2,N1[0]==1,N2[0]==1}/.val,{N1[t],N2[t]},{t,0,T}];
Plot[Evaluate[{N1[t],N2[t]}/.sol/.val],{t,0,T}]

Koska kaikki ympäristön ja muiden lajien vaikutukset on nollattu, jäljelle jää rajoittamaton eksponentiaalinen kasvu.

b

T=20;
val={ϵ1->0.3,ϵ2->0.4,β1->0.2,β2->0.2,γ1->0.0,γ2->0.0};
sol=NDSolve[{eq1,eq2,N1[0]==1,N2[0]==1}/.val,{N1[t],N2[t]},{t,0,T}];
Plot[Evaluate[{N1[t],N2[t]}/.sol/.val],{t,0,T}]

Nyt ympäristö asettaa kummankin populaation kasvulle rajat, joten kasvu taittuu jossain pisteessä. Populaatiot kehittyvät edelleen toisistaan riippumattomasti.

c

T=20;
val={ϵ1->0.3,ϵ2->0.4,β1->0.2,β2->0.2,γ1->0.1,γ2->0.1};
sol=NDSolve[{eq1,eq2,N1[0]==1,N2[0]==1}/.val,{N1[t],N2[t]},{t,0,T}];
Plot[Evaluate[{N1[t],N2[t]}/.sol/.val],{t,0,T}]

Nyt populaatiot varjostavat toisiaan yhtä voimakkaasti, ja ympäristön rajoitteet ovat samanlaiset molemmille. Tällaisessa tilanteessa se populaatio, jonka syntyvyys on suurempi, pääsee niskan päälle.

d

T=70;
val={ϵ1->0.3,ϵ2->-0.4,β1->0.0,β2->0.0,γ1->0.1,γ2->-0.15};
sol=NDSolve[{eq1,eq2,N1[0]==1,N2[0]==1}/.val,{N1[t],N2[t]},{t,0,T}];
Plot[Evaluate[{N1[t],N2[t]}/.sol],{t,0,T}]

Nyt toinen populaatio elää ensimmäisen "kustannuksella" (esim. saalistaja). Kun 1. populaation koko kasvaa, toinenkin lähtee kasvuun, joka taas aiheuttaa 1. populaation pienenemisen. Kun 1. populaatio romahtaa, toinen alkaa "kuolla nälkään", jne.

ParametricPlot[Evaluate[{N1[t],N2[t]}/.sol/.val],{t,0,T},AspectRatio->1/GoldenRatio]

"The Circle of Life"