Harjoitus 4 -- Ratkaisut

1

Ei kommenttia.

2

a=Table[{x,Sin[x]}, {x,-2*Pi,2*Pi, 0.5}]

b=Interpolation[a, InterpolationOrder->0]
c=Interpolation[a, InterpolationOrder->1]
d=Interpolation[a, InterpolationOrder->2]

Plot[{a[x],b[x],c[x]},{x,-1.9*Pi,1.9*Pi}]

3

Tutkittava funktio:

f[x_ ]:=(x/(1+x ((Sin[x])^2)));
Plot[f[x],{x,0,3π },PlotRange->All]

Luodaan tasavälinen pisteistö välille 0<x<3π. Tehdään se ensin kiinnitetyllä n:n arvolla. Pisteistön luominen erikseen ei ole välttämätöntä tässä -- voisi käyttää suoraan Table-komentoa lukuparitaulukon luomiseen. Jos välillä [(Subscript[x, 1]),(Subscript[x, n])] on n pistettä, vierekkäisten pisteiden välimatka on (((Subscript[x, n])-(Subscript[x, 1]))/(n-1)) (osoita!).

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

Muodostetaan taulukko lukupareista ((Subscript[x, i]),f((Subscript[x, i]))):

t=Table[{x[i],f[x[i]]},{i,n}];
Show[ListPlot[t],Plot[f[x],{x,0,3π},PlotRange->All]]

Taulukon t olisi voinut muodostaa viittaamatta muuttujiin x[i].

t=Table[{x,f[x]},{x,0,3π,((3π)/(n-1))}];

Luodaan tähän taulukkoon interpoloiva polynomi harjoituksissa annetulla kaavalla. Käytetään tavallisia toisto- ja valintarakenteita (ks. edellisistä harjoituksista summan ja tulon laskemisen periaate). Nyt laskettavan summan jokainen termi on tulo, joten tarvitaan kaksi sisäkkäistä toistorakennetta. Huom. muuttuja tulo on alustettava ykköseksi vasta sisemmän For-silmukan alussa!

p[y_]:=Module[(*paikalliset muuttujat*){i,j,summa,tulo},
(*1. For-silmukka laskee summaa*)
For[i=1;summa=0,i≤n,i++,
(*2. For-silmukka laskee tuloa*)
For[j=1;tulo=1,j≤n,j++,
If[j≠i,tulo*=((y-x[j])/(x[i]-x[j]))]
];
summa+=f[x[i]]*tulo;];
(*viimeisen käskyn tulos on funktion arvo*)summa]

Polynomin p aste on n-1. Expand-funktio laskee kaikki kertolaskut auki.

p[y]//Expand//N

Plot[p[y],{y,0,3π},PlotRange->{-π,4π}]

Polynomin voi muodostaa myös Mathematican valmiilla funktioilla Sum ja Product, mutta toistorakenteilla tehtäessä periaate on sama kuin lähes kaikissa ohjelmointikielissä, joten se kannattaa opetella.

q[y_]:=Sum[f[x[i]]Product[If[j≠i,((y-x[j])/(x[i]-x[j])),1(*jos j==i, niin sitä j:n arvoa vastaava tulon tekijä on 1*)],{j,1,n}],{i,1,n}];

q[y]//Expand//N

Plot[q[y],{y,0,3π},PlotRange->{-π,4π}]

Tässä voitaisiin käyttää myös taulukkoa t. Merkintä t[[i,j]] antaa taulukon t i:nnen alkion (joka on kahden alkion kokoinen lista) j:nnen alkion (j=1 antaa (Subscript[x, i]):n ja j=2 antaa f((Subscript[x, i])):n).

q[y_]:=Sum[t[[i,2]]Product[If[j≠i,((y-t[[j,1]])/(t[[i,1]]-t[[j,1]])),1],{j,1,n}],{i,1,n}];

Muodostetaan interpoloiva polynomi eri n:n arvoilla. Huom. tässä Evaluate-funktion käyttäminen todella kannattaa, piirtäminen kestää muuten törkeän kauan, jos For-silmukat täytyy käydä läpi jokaisella y:n arvolla erikseen. PlotStyle->{Hue[x]}-komento piirtää kuvaajan x:n määräämällä värisävyllä (ks. Hue (paina F1)).

i=0;
Do[
Clear[x];x[1]=0;x[n]=3π; h=((x[n]-x[1])/(n-1));
Do[x[i]=x[1]+(i-1)h,{i,2,n-1}];
p[y_]:=Module[{i,j,summa,tulo},
For[i=1;summa=0,i≤n,i++,
For[j=1;tulo=1,j≤n,j++,
If[j≠i,tulo*=((y-x[j])/(x[i]-x[j]))]];
summa+=f[x[i]]*tulo;];
summa];
g[n]=Plot[Evaluate[p[y]],{y,0,3π},PlotRange->{-π,4π},PlotStyle->{Hue[((i++)/10)]}];,
{n,{2,5,10,20,50}} (*n saa annetussa listassa olevat arvot*)]

Nyt interpoloivien polynomien kuvaajat ovat indeksoiduissa muuttujissa g[n].

Show[Table[g[n],{n,{2,5,10,20,50}}]]

Huomataan, että suurilla n:n arvoilla interpoloiva polynomi oskilloi hyvin voimakkaasti datapisteistön reunoilla.

Sama lasku Mathematican valmiin funktion, InterpolatingPolynomial, avulla: Nyt tarvitaan vain taulukkoa datapisteistöstä, t, joka pitää muodostaa kullekin n:n arvolle erikseen. InterpolatingPolynomial-funktion 1. argumentti on datapisteistö ja 2. argumentti on muuttujan nimi, jota polynomissa tulee käyttää.

n=5;
t=Table[{x,f[x]},{x,0,3π,((3π)/(n-1))}];
r[y_]:=InterpolatingPolynomial[t,y];
r[y]//Expand//N
Plot[r[y],{y,0,3π},PlotRange->{-π,4π}]

Havaitaan, että valmis funktio muodostaa polynomin lähes täsmälleen samalla tavalla.

Muilla n:n arvoilla:

i=0;
Do[
t=Table[{x,f[x]},{x,0,3π,((3π)/(n-1))}];
r[y_]:=InterpolatingPolynomial[t//N,y];
w[n]=Plot[r[y],{y,0,3π},PlotRange->{-π,4π},PlotStyle->{Hue[((i++)/10)]}];,
{n,{2,5,10,20,50}}];

Show[Table[w[n],{n,{2,5,10,20,50}}]]

Vertailun vuoksi vielä vierekkäin. Eroa ei huomaa.

GraphicsArray[{Show[Table[w[n],{n,{2,5,10,20,50}}]],Show[Table[g[n],{n,{2,5,10,20,50}}]]}]

4

Paloittain määritelty interpoloiva funktio: (Subscript[q, i])(x)=(Subscript[a, i])(x-(Subscript[x, i]))+(Subscript[b, i]), kun (Subscript[x, i])<x<(Subscript[x, i+1]). Kertoimet (Subscript[a, i]) ja (Subscript[b, i]) määräytyvät ehdoista (Subscript[q, i])((Subscript[x, i]))=(Subscript[b, i])=(Subscript[f, i]) ja (Subscript[q, i])((Subscript[x, i+1]))=(Subscript[a, i])((Subscript[x, i+1])-(Subscript[x, i]))+(Subscript[b, i])=(Subscript[f, i+1]), josta (Subscript[a, i])=(((Subscript[f, i+1])-(Subscript[f, i]))/((Subscript[x, i+1])-(Subscript[x, i])))=(((Subscript[f, i+1])-(Subscript[f, i]))/h).

Muodostetaan erikseen indeksoidut muuttujat x[i] x:n arvoille. f(x):n arvot saadaan näistä laskemalla f[x[i]]. Tämän voisi tehdä myös käyttämällä taulukkoa t.

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

Muodostetaan paloittain määritelty funktio Q, joka kulkee pisteiden {x[i],f[x[i]]} kautta. Periaate on seuraava: Verrataan funktion argumentin arvoa järjestyksessä jokaiseen (Subscript[x, i]):hin, kun i=1,...,n. Jos x≥(Subscript[x, i]) ja x≤(Subscript[x, i+1]), niin asetetaan Q(x)=(Subscript[q, i])(x), ja (Subscript[q, i]) määritetään edellä laskettujen (Subscript[a, i]):n ja (Subscript[b, i]):n avulla. Merkitään argumenttia y:llä, sillä x on varattu taulukon nimeksi. Seuraava toteutus asettaa funktion Q arvoksi 0 datapisteistön ulkopuolella.

Q[y_]:=Module[{i,qx=0},
For[i=1,i≤n-1,i++,
If[y≥x[i]&&y≤x[i+1],
qx=((f[x[i+1]]-f[x[i]])/h)(y-x[i])+f[x[i]]//N;Break[]]];
qx];

Huom. nyt plottauksessa ei saa käyttää Evaluate-komentoa, sillä Q:n lauseke ei ole sama kaikilla y:n arvoilla.

Plot[Q[y],{y,0,3π}]

Tehdään sama Mathematican valmiilla funktiolla Interpolation. Se ottaa 1. argumentikseen datapisteistön (taulukon t), ja antaa tulokseksi ns. "puhtaan funktion". Puhtaista funktioista ei tarvitse tällä kurssilla tietää muuta kuin että niitä käytetään kirjoittamalla funktion nimen perään hakasulkeisiin argumentti. Interpolation-funktiolle voidaan antaa optiona interpolointiaste (ts. polynomien (Subscript[q, i]) aste).

t=Table[{x,f[x]},{x,0,3π,((3π)/(n-1))}];

s=Interpolation[t,InterpolationOrder->1]
Plot[s[y],{y,0,3π}]

Interpolaatio 3. asteen polynomeilla (se on oletusasteluku Interpolation-funktiossa):

s=Interpolation[t]
Plot[s[y],{y,0,3π}]

Havaitaan, että derivaatta ei ole jatkuva. Spline-interpolaatio tekee 1. ja 2. derivaatoista jatkuvia. Spline-funktio täytyy ladata paketista Splines:

<<Splines`

Nyt s[y] antaa lukuparin {sx[y],sy[y]}, joka voidaan piirtää ParametricPlot-komennolla.

s=SplineFit[t,Cubic];
ParametricPlot[s[y],{y,0,3π},AspectRatio->1/GoldenRatio]