T4 mikrolaskarit 20-21.2 1997

teht. 1

Remove["Global`*"];

f=x (1-x) Exp[x]; u[n_]:=Sin[n Pi x]; l[n_]:=(n Pi)^2+1;
kerroin=2*Integrate[f*u[n], {x,0,1}]

Maaritellaan f:n Fourier-kertoimet n:n funktiona:

a[n_]:=kerroin/.{Cos[n Pi]->(-1)^n,Sin[n Pi]->0}; a[n]

Piirretaan nyt f ja viiden termin summa samaan kuvaan.

s5=Sum[a[n] u[n],{n,5}]; Plot[{f,s5},{x,0,1}]

Approksimaatio on varsin tarkka; katsotaan viela tehtya virhetta.

Plot[f-s5,{x,0,1}]

Ratkaistaan nyt tehtava tarkasti. Vastaus saadaan nk. puhtaana funktiona kun toiseksi argumentiksi laitetaan y[x]:n sijaan y.

DSolve[{-y''[x]+y[x]==f,y[0]==0,y[1]==0},y,x]

Annetaan ratkaisulle nimi tarkka1. Tata tarvitaan viela tehtavassa 3!

tarkka1=y/.%[[1]]

Piirretaan tarkka ratkaisu ja viiden termin sinisarja samaan kuvaan. Katsotaan myos erotusta ja huomataan, etta tarkkuus (myos suhteellinen!) on viela parempi kuin f:lle.

c[n_]:=a[n]/l[n]; s5=Sum[c[n] u[n],{n,5}];

Plot[{s5,tarkka1[x]},{x,0,1}]; Plot[tarkka1[x]-s5,{x,0,1}];

teht. 2

f=x^(-1/3); Integrate[f^2,{x,0,1}]

u[n_]:=Sin[n Pi x]; l[n_]:=(n Pi)^2;
a[n_]:=2*Integrate[f u[n], {x,0,1}]; a[n]

Nama Fourier-kertoimet ovat ikavia laskea:

N[a[30]]
N[a[30],60]
N[a[30],70]

Viiden termin summa ei nyt nayta aivan samalta kuin f.

DSolve[{-y''[x]==f,y[0]==0,y[1]==0},y,x]

tarkka2=y/.%[[1]]

Ratkaisun Fourier-osasumma on paljon tarkempi, koska u on silea toisin kuin f.

c[n_]:=a[n]/l[n]; s5=Sum[c[n] u[n],{n,5}];
Plot[{tarkka2[x],s5},{x,0,1}]; Plot[tarkka2[x]-s5,{x,0,1}]

teht. 3

Clear[f]; f[x_]:=x (1-x) Exp[x];

Ratkaistaan tehtavien 3 ja 4 yhtaloryhmat TridiagonalSolvella, joka soveltuu juuri tahan tarkoitukseen (ja vain tahan...)

<<LinearAlgebra/Tridiagonal.m;
?TridiagonalSolve

Sitten vain valitaan h, muodostetaan matriisin diagonaali ja diagonaalin viereinen nauha seka yhtalon oikea puoli ja ratkaistaan yhtaloryhma. Huomaa, etta diagonaalille tulee 2+h^2.

h=0.2; n=Round[1/h-1]; a=Table[-1,{n-1}]; b=Table[2+h^2,{n}];
r=h^2 Table[f[i h],{i,n}]; u=TridiagonalSolve[a,b,a,r]

Lisataan viela reunaehdot.

PrependTo[u,0]; AppendTo[u,0]

ListPlot tahtoo syotteen tassa muodossa:

x=Table[i h,{i,0,n+1}]; lista=Transpose[{x,u}]

kuva1=ListPlot[lista,Prolog->{PointSize[0.02]},
  DisplayFunction->Identity];
kuva2=Plot[tarkka1[x],{x,0,1},DisplayFunction->Identity]; 
Show[kuva1,kuva2,DisplayFunction->$DisplayFunction]

Lasketaan sitten arvio vakiolle c. Laitetaan tama ensimmaiseksi alkioksi listaan, johon lisataan myohemmat estimaatit.

c={ Max[Table[Abs[u[[i]]-N[tarkka1[i*h]]],{i,n}]]/h^2 }

Puolitetaan h ja  toistetaan samat vaiheet.

h=h/2; n=Round[1/h-1]; a=Table[-1,{n-1}]; b=Table[2+h^2,{n}];
r=h^2 Table[f[i h],{i,n}]; u=TridiagonalSolve[a,b,a,r];
PrependTo[u,0]; AppendTo[u,0];
x=Table[i h,{i,0,n+1}]; lista=Transpose[{x,u}];
kuva1=ListPlot[lista,Prolog->{PointSize[0.02]},
  DisplayFunction->Identity];
kuva2=Plot[tarkka1[x],{x,0,1},DisplayFunction->Identity]; 
Show[kuva1,kuva2,DisplayFunction->$DisplayFunction]

AppendTo[c,Max[Table[Abs[u[[i]]-N[tarkka1[i*h]]],{i,n}]]/h^2]

Kuva nayttaa jo hyvalta, keskitytaan nyt c:n arvioon. Puolitetaan h muutamia kertoja ja katsotaan miten c konvergoi.

Do[
   h=h/2; n=Round[1/h-1]; a=Table[-1,{n-1}]; b=Table[2+h^2,{n}];
   r=h^2 Table[f[i h],{i,n}]; u=TridiagonalSolve[a,b,a,r];
   AppendTo[c,
     Max[Table[Abs[u[[i]]-N[tarkka1[i*h]]],{i,n}]]/h^2],
{5}]; 
c

teht. 4

Clear[f]; f[x_]:=x^(-1/3);

Aivan kuten edella, mutta matriisin diagonaalilla on nyt pelkka 2.

h=0.2; n=Round[1/h-1]; a=Table[-1,{n-1}]; b=Table[2,{n}];
r=h^2 Table[f[i h],{i,n}]; u=TridiagonalSolve[a,b,a,r];
PrependTo[u,0]; AppendTo[u,0];
x=Table[i h,{i,0,n+1}]; lista=Transpose[{x,u}];
kuva1=ListPlot[lista,Prolog->{PointSize[0.02]},
  DisplayFunction->Identity];
kuva2=Plot[tarkka2[x],{x,0,1},DisplayFunction->Identity]; 
Show[kuva1,kuva2,DisplayFunction->$DisplayFunction]

Estimaatti eksponentille alpha saadaan kun tunnetaan virheet arvoilla h ja h/2. Lasketaan naita useita puolittamalla aina h:ta.

alpha={}; vanha=Max[Table[Abs[u[[i]]-tarkka2[i*h]],{i,n}]]

Do[
   h=h/2; n=Round[1/h-1]; a=Table[-1,{n-1}]; b=Table[2,{n}];
   r=h^2 Table[f[i h],{i,n}]; u=TridiagonalSolve[a,b,a,r];
   uusi=Max[Table[Abs[u[[i]]-tarkka2[i*h]],{i,n}]];
   AppendTo[alpha,Log[2,vanha/uusi]]; vanha=uusi,
{6}];
alpha

Jono konvergoi kohti 5/3:a, mutta melko hitaasti. Huomaa tarkan ratkaisun termi x^(5/3)!

teht. 5

Tassa tehtavassa ratkaistaan tridiagonaalimatriisin ominaisarvot. Matriisi on nyt oikeasti luotava (pelkat nauhat ei riita). Kaytetaan tehtavapaperissa annettua keinoa.

<<LinearAlgebra/MatrixManipulation.m;

Lasketaan matriisiin virhe kolmen pienimman ominaisarvon virheet h:n arvoilla 1/10, 1/20, 1/40 ja 1/80.

h=0.1; virhe={};
Do[
   n=Round[1/h-1]; ve=Table[0,{n-2}]; PrependTo[ve,-1];
   PrependTo[ve,2]; A=Reverse[HankelMatrix[Reverse[ve],ve]];
   A[[n,n]]=1; A=A/h^2; spektri=Sort[Eigenvalues[A]];
   AppendTo[virhe,
            Table[N[spektri[[j]]-Pi^2 (j-0.5)^2],{j,3}]];
   h=h/2,
{i,4}];
MatrixForm[virhe]

Virhe puolittuu h:n puolittuessa (sita tarkemmin mita pienempi ominaisarvo ja h). Siis virhe on O(h). Kokeillaan sitten reunaehdon tarkempaa diskretointia.

h=0.1; virhe={};
Do[
   n=Round[1/h-1]; ve=Table[0,{n-2}]; PrependTo[ve,-1];
   PrependTo[ve,2]; A=Reverse[HankelMatrix[Reverse[ve],ve]];
   A[[n,n-1]]=-2/3; A[[n,n]]=2/3; A=A/h^2;
   spektri=Sort[Eigenvalues[A]];
   AppendTo[virhe,
            Table[N[spektri[[j]]-Pi^2 (j-0.5)^2],{j,3}]];
   h=h/2,
{i,4}];
MatrixForm[virhe]

Nyt virhe putoaa pikemminkin neljasosaan eli virhe on O(h^2). Kokeillaan viela diskretointi
u'(1)=(u(1+h)-u(1-h))/2h.

h=0.1; virhe={};
Do[
   n=Round[1/h-1]; ve=Table[0,{n-2}]; PrependTo[ve,-1];
   PrependTo[ve,2]; A=Reverse[HankelMatrix[Reverse[ve],ve]];
   A[[n,n-1]]=-2; A[[n,n]]=2; A=A/h^2;
   spektri=Sort[Eigenvalues[A]];
   AppendTo[virhe,
            Table[N[spektri[[j]]-Pi^2 (j-0.5)^2],{j,3}]];
   h=h/2,
{i,4}];
MatrixForm[virhe]

Selvasti O(h).