Harjoitus 5 -- Ratkaisut

1

Ei kommenttia.

2

Tutkittava funktio oskilloi äärettömän tiheään nollan lähellä. PlotPoints-asetus määrää, kuinka tiheästi Plot-funktio ottaa piirrettävästä funktiosta "näytteitä" kuvaajaa varten.

f[x_]:=Sin[(1/x)];
Plot[f[x],{x,0.01,2},PlotRange->{-1,1},AxesOrigin->{0,0},PlotPoints->100]

Luodaan pisteet (Subscript[y, i]) (käytetään y:tä muuttujana, sillä x on varattu funktion f argumentiksi):

a=0.1;b=2.0;n=20;h=((b-a)/(n-1));y[1]=a;y[n]=b;
Do[y[i]=a+(i-1) h,{i,2,n-1}];

Show[Plot[f[x],{x,a,b},PlotRange->{-1,1},AxesOrigin->{0,0}],ListPlot[Table[{y[i],f[y[i]]},{i,1,n}]]]

Ongelmallinen tapaus, sillä tämä h:n arvo riittäisi mainiosti alueeseen x≥0.5, mutta ei lähempänä origoa. Tämän voisi korjata luomalla epätasavälisen pisteistön: (Subscript[x, i+1])-(Subscript[x, i])=(Subscript[h, i]) tms.

Luodaan taulukko derivaattojen arvoista:

df=Table[{y[i],((f[y[i+1]]-f[y[i-1]])/(2h))},{i,2,n-1}];
ddf=Table[{y[i],((f[y[i+1]]-2f[y[i]]+f[y[i-1]])/(h^2))},{i,2,n-1}];

Show[Plot[f'[x],{x,a,b},PlotRange->{-45,20},AxesOrigin->{0,0}],ListPlot[df]]

Show[Plot[f''[x],{x,a,b},PlotRange->{-100,1000},AxesOrigin->{0,0}],ListPlot[ddf]]

Tämä h:n arvo ei riitä mihinkään. Oikeastaan mikään h:n arvo (tai mikään n:n arvo) ei riitä, jos tarpeeksi lähelle origoa mennään.

Clear[a,b,n,h]

a=0.05;b=0.2;n=40;h=((b-a)/(n-1));y[1]=a;y[n]=b;
Do[y[i]=a+(i-1) h,{i,2,n-1}];

df=Table[{y[i],((f[y[i+1]]-f[y[i-1]])/(2h))},{i,2,n-1}];
ddf=Table[{y[i],((f[y[i+1]]-2f[y[i]]+f[y[i-1]])/(h^2))},{i,2,n-1}];

g1=Show[Plot[f'[x],{x,a,b},PlotRange->All],ListPlot[df]];
g2=Show[Plot[f''[x],{x,a,b},PlotRange->All],ListPlot[ddf]];
GraphicsArray[{g1,g2}]

3

Differentiaaliyhtälöiden ratkaisuun perehdytään myöhemmin. Tässä vähän maistiaisia:

Clear[n,t,λ]

ratk=DSolve[{n'[t]==-λ n[t],n[0]==n0}(*ratkaistava yhtälö(ryhmä) aaltosulkeisiin, myös alkuehdot tänne*),
n[t](* ratkaistava funktio*),
t (*riippumaton muuttuja*)]

Ratkaisu saadaan sääntönä (transformation rule). Myös parametrit (Subscript[n, 0]) ja λ kannattaa antaa sääntöjen kautta (näin niihin ei tule pysyvästi mitään arvoa). Seuraavassa sovelletaan ensin sääntöä ratk, jonka jälkeen sijoitetaan (Subscript[n, 0]):n ja λ:n arvot.

Plot[n[t]/.ratk/.{n0->(10^25),λ->0.1},{t,0,50}]

Kun differentiaaliyhtälöön sijoitetaan derivaatan approksimaatio, saadaan (Subscript[N, i+1])=(1-λ h)(Subscript[N, i]). Olkoon h=0.1, ja (Subscript[t, i])=(i-1)h, kun i=1,...,m.

h=0.1;m=200;
Do[t[i]=(i-1)h,{i,1,m}];

Taulukoidaan arvot (Subscript[n, i]): (merkitään n:ää μ:llä).

μ[1]=(10^25);λ=0.1;
Do[μ[i]=(1-λ h)μ[i-1],{i,2,m}];

Show[ListPlot[Table[{t[i],μ[i]},{i,m}]],Plot[n[t]/.ratk/.{n0->(10^25),λ->0.1},{t,0,m h},PlotStyle->Hue[1]]]

Likimääräinen ja tarkka ratkaisu käyvät hyvin yksiin. Suuremmilla h:n arvoilla (tai pienemmillä m:n arvoilla) havaitaan poikkeamaa:

Clear[t,μ];

m=40;t[1]=0;t[m]=40;h=((t[m]-t[1])/(m-1));
Do[t[i]=t[1]+(i-1)h,{i,1,m}];

μ[1]=(10^25);λ=0.1;
Do[μ[i]=(1-λ h)μ[i-1],{i,2,m}];

Show[ListPlot[Table[{t[i],μ[i]},{i,m}]],Plot[n[t]/.ratk/.{n0->(10^25),λ->0.1},{t,0,m h},PlotStyle->Hue[1]]]

Tässä ratkaisumenetelmässä virhe siis kasaantuu joka aika-askeleella.

4

Otetaan yksiköt mukaan. Jälleen aikamuuttujaan t on liitettävä yksikkö Second.

<<Units`

Clear[g,a,v,x,t,c]

c=0.1Meter (Second^(-3));
g[t_]:=c t Second;

Nopeus on kiihtyvyyden integraali. Merkitään integroimismuuttujaa s:llä. Alaviiva funktion määrittelyssä (g[t_]) tarkoittaa, että g hyväksyy argumentikseen minkä tahansa lausekkeen. Kun integroidaan ajan suhteen, on yksiköitä kerrottava sekunnilla.

v0=1Meter (Second^(-1));
v[t_]:=v0+Integrate[g[s],{s,0,t}]Second;

Paikka on nopeuden integraali.

x0=0;
x[t_]:=x0+Integrate[v[s],{s,0,t}]Second//Simplify;

Plot[Evaluate[x[t]/Meter],{t,0,10},ImageSize->Small]

Toinen lähestymistapa: (huom. sääntöjä (a->b) voidaan soveltaa myös funktioihin, ja niillä voidaan sijoittaa funktioihin arvoja)

Clear[g,v,x,t,s,x0,v0]

x[t_]:=x0+Integrate[v0+Integrate[g[s],{s,0,r}],{r,0,t}];

Plot[Evaluate[x[t]/.g[s]->0.1s/.{v0->1,x0->0}],{t,0,10},ImageSize->Small]

Määrätyn integraalin voi kirjoittaa esc-dintt-esc (määräämätön esc-int-esc):

g[t_]:=0.1t;x0=0;v0=1;
x[t_]:=x0+(∫_0^t)(v0+(∫_0^s)g[r]\[DifferentialD]r)\[DifferentialD]s

Plot[Evaluate[x[t]],{t,0,10},ImageSize->Small]

b) Puolisuunnikassääntö

Katso luennoista puolisuunnikassäännön määritelmä. Lasketaan paikka ajanhetkillä (Subscript[t, i]), i=1,...,n; (Subscript[t, 1])=0 ja (Subscript[t, n])=10s.

n=20;t[1]=0;t[n]=10;h=((t[n]-t[1])/(n-1));
Do[t[i]=t[1]+(i-1)h,{i,2,n-1}]

Oletetaan, että kiihtyvyys tunnetaan vain hetkillä (Subscript[t, i]).

Do[g[i]=0.1t[i],{i,1,n}]

Paikan laskemiseksi hetkellä (Subscript[t, i]) tarvitaan nopeuden arvot hetkillä (Subscript[t, j]), missä j=1,...,i. Taulukoidaan siis ensin nopeuden arvot kaikilla ajanhetkillä. Nopeus hetkellä (Subscript[t, 1]) saadaan alkuehdosta. Loput lasketaan puolisuunnikaskaavalla kiihtyvyydestä:
	v((Subscript[t, i]))=v((Subscript[t, 1]))+(Overscript[(∫_(Subscript[t, 1])), (Subscript[t, i])])a(t)dt≈(Subscript[v, 1])+h(∑_(j=1)^(i-1))(((Subscript[a, i])+(Subscript[a, i+1]))/2).
Koska nyt halutaan kaikki nopeuden arvot, yo. kaava on kätevintä kirjoittaa muotoon
	v((Subscript[t, i]))=(Subscript[v, 1])+h(∑_(j=1)^(i-2))(((Subscript[a, j])+(Subscript[a, j+1]))/2)+h(((Subscript[a, i-1])+(Subscript[a, i]))/2)=(Subscript[v, i-1])+h(((Subscript[a, i-1])+(Subscript[a, i]))/2).

v[1]=1; (*alkuehto v0=1m/s*)
Do[v[i]=v[i-1]+h ((g[i-1]+g[i])/2),{i,2,n}]

Summan voi tietenkin laskea jokaiselle i:n arvolle erikseen:

Do[Δv=0;Do[Δv+=((g[j]+g[j+1])/2),{j,1,i-1}];v[i]=v[1]+h Δv,{i,2,n}]

ListPlot[Table[{t[i],v[i]},{i,1,n}],ImageSize->Small]

Paikka voidaan nyt laskea nopeudesta, samalla tavalla kuin äsken nopeus kiihtyvyydestä.

x[1]=0;
Do[x[i]=x[i-1]+h((v[i-1]+v[i])/2),{i,2,n}]

ListPlot[Table[{t[i],x[i]},{i,1,n}],ImageSize->Small]