Harjoitus 9 -- ratkaisut

1

Sijoitetaan yhtälön ratkaisufunktio funktioksi r.

eq=y''[x]+4x y'[x]+(102+4(x^2))y[x]==0;
r[x_]=y[x]/.Flatten[DSolve[{eq,y[0]==0,y'[0]==1},y[x],x]]

Simplify-funktio ei osaa käsitellä kompleksisia eksponentteja. ComplexExpand laskee ne auki Eulerin kaavalla (e^(i x))=cos x+i sin x:

r[x_]=ComplexExpand[r[x]]//Simplify

Funktio onkin reaalinen!

Plot[r[x],{x,-3,3},PlotRange->All]

Funktion kuvaaja on sinikäyrä, jota moduloi Gaussin funktio (e^(-(x^2))). Tämän kaltaisia funktioita saattaa tulla vastaan esim. aaltoliikeopissa, optiikassa tai kvanttimekaniikassa (aaltopaketti).

2

Differenssimenetelmässä lähdetään liikkeelle siitä, että määrätään ne pisteet (muuttujan arvot), joissa yhtälölle etsitään ratkaisua. Yhtälön ratkaisu halutaan selvittää alueessa -3<x<3, kun tiedetään funktion käyttäytyminen pisteessä x=0. Ratkaisuprosessin on alettava pisteestä x=0 ja edettävä positiiviseen ja negatiivisene suuntaan. Prosessin alkuun voidaan päästä monella tavalla.

Olkoot (Subscript[x, 0])=0 ja (Subscript[x, n])=3. Käytetään ensimmäisen indeksin arvona nollaa, sillä siitä on helppo lähteä sekä positiiviseen että negatiiviseen suuntaan. Alueessa 0<x≤3 ratkaisua edustaa (Subscript[y, i]), i=1,...,n ja alueessa -3≤x<0 ratkaisu on (Subscript[y, i]) kun i=-1,...,-n.  Muodostetaan differenssiyhtälö ja ratkaistaan (Subscript[y, i+1]):
	(Subscript[y, i+1])=(1/((1/(h^2))+((2(Subscript[x, i]))/h)))((Subscript[y, i])((2/(h^2))-102-4((Subscript[x, i])^2))+(Subscript[y, i-1])(((2(Subscript[x, i]))/h)-(1/(h^2))))	(1)
ja (Subscript[y, i-1]):
	(Subscript[y, i-1])=(1/((1/(h^2))-((2(Subscript[x, i]))/h)))((Subscript[y, i])((2/(h^2))-102-4((Subscript[x, i])^2))-(Subscript[y, i-1])(((2(Subscript[x, i]))/h)+(1/(h^2))))	(29
Alkuehtoina annettiin y(0)=0 ja y'(0)=1, eli (Subscript[y, 0])=0 ja (Subscript[y, ±1])=(Subscript[y, 0])±y'(0)h=(Subscript[y, 0])±h. Huom. nyt (Subscript[x, i])=(Subscript[x, 0])+i h, eikä (Subscript[x, 0])+(i-1)h, kuten aiemmin.

(huom. kannattaa kirjoittaa desimaalipisteet eksplisiittisesti näkyviin. Muuten Mathematica alkaa murskata murtolukuja ja laskussa kestää kauan.)

n=1000;x[0]=0.0;x[n]=3.0;x[-n]=-3;h=((x[n]-x[0])/(n-1));
Do[
x[i]=x[0]+i h;
x[-i]=x[0]-i h;,
{i,1,n-1}]

Alkuehdot:

y0=0;dy0=1;
y[0]=y0;
y[1]=y[0]+dy0 h;
y[-1]=y[0]-dy0 h;

Varsinainen ratkaisuprosessi: Positiiviseen suuntaan etenevä ratkaisu voidaan toteuttaa yhtälöllä (1), ja negatiivinenkin voitaisiin toteuttaa omassa Do-lausekkeessaan suoraan yhtälöllä (2). Toinen tapa on yhdistää Do-lausekkeet ja kirjoittaa negatiiviseen suuntaan mentäessä i->-i, i-1->-(i-1)=-i+1 jne.

Do[
y[i]=(1/((1/(h^2))+((2x[i-1])/h)))(y[i-1]((2/(h^2))-102-4((x[i-1])^2))+y[i-2](-(1/(h^2))+((2x[i-1])/h)));
y[-i]=(1/((1/(h^2))-((2x[-i+1])/h)))(y[-i+1]((2/(h^2))-102-4((x[-i+1])^2))-y[-i+2]((1/(h^2))+((2x[-i+1])/h)));,
{i,2,n}]

Negatiiviseen suuntaan etenevä Do-lauseke voisi myös olla muotoa (huom. kun indeksin arvot pienenevät, on annettava askelpituus)

Do[
y[i]=(1/((1/(h^2))-((2x[i+1])/h)))(y[i+1]((2/(h^2))-102-4((x[i+1])^2))-y[i+2]((1/(h^2))+((2x[i+1])/h)));,
{i,-2,-n}]

Piirretään likiratkaisun kuvaaja ja tarkka ratkaisu samaan kuvaajaan. Likiratkaisua ei erota tarkasta.

Show[ListPlot[Table[{x[i],y[i]},{i,-n,n}],PlotRange->All],Plot[r[t],{t,-3,3},PlotRange->All]]

Pienennetään n:n arvoa ja katsotaan, kuinka likiratkaisun käy. Muodostetaan ratkaisuprosessista funktio, joka ottaa argumentikseen tutkittavan n:n arvon.

examinesolution[m_]:=(
Clear[x,y];
n=m;
x[0]=0.0;x[n]=3.0;x[-n]=-3;h=((x[n]-x[0])/(n-1));
Do[x[i]=x[0]+i h;x[-i]=x[0]-i h,{i,1,n-1}];
y0=0;dy0=1;y[0]=y0;y[1]=y[0]+dy0 h;y[-1]=y[0]-dy0 h;
Do[
y[i]=(1/((1/(h^2))+((2x[i-1])/h)))(y[i-1]((2/(h^2))-102-4((x[i-1])^2))+y[i-2](-(1/(h^2))+((2x[i-1])/h)));
y[-i]=(1/((1/(h^2))-((2x[-i+1])/h)))(y[-i+1]((2/(h^2))-102-4((x[-i+1])^2))-y[-i+2]((1/(h^2))+((2x[-i+1])/h)));,
{i,2,n}];
Show[ListPlot[Table[{x[i],y[i]},{i,-n,n}],PlotRange->All],Plot[r[t],{t,-3,3},PlotRange->All]])

GraphicsGrid toimii samalla tavalla kuin GraphicsArray, jota on käytetty jo aiemmin.

GraphicsGrid[{Table[examinesolution[i],{i,{500,200,100}}],
Table[examinesolution[i],{i,{50,20,18}}]}]

Kun n≃18, alkaa ratkaisu leviämään reuna-alueella käsiin. Tätä vastaava h:n arvo on

h18=(3/(18-1))//N

Tarkka ratkaisu oli muotoa
	(1/10)(e^(-(x^2)))sin 10x,
ts. se oskilloi aaltoluvulla k=10. Tästä saadaan aallonpituus λ=((2π)/k):

λ=((2π)/10)//N

Nyt h:n arvo on n. 1/3...1/4 aallonpituutta.

λ/h18

Jos oskilloivasta funktiosta otetaan näytteitä askelpituudella h<λ/4, oskillaatiot näkyvät vielä hyvin. Jos askelpituutta tästä vielä kasvatetaan, värähtelyjen muoto vääristyy ja saattaa jopa kadota kokonaan.

g[i_]:=Show[Plot[Sin[z],{z,0,4π},PlotStyle->Red],ListPlot[Table[{z,Sin[z]},{z,0,4π,((2π)/i)}],Joined->True]];

GraphicsGrid[{Table[g[i],{i,{9,8.2,7.5}}],Table[g[i],{i,{6.3,5.02,4.1}}],Table[g[i],{i,{λ/h18,2.8,2.2}}]}]

Edellisen kuvakollaasin viimeistä kuvaa vastaava askelpituus näkyy vasemmassa alakulmassa.