T4 mikroharjoitukset -97, osa II

teht. 1

Remove["Global`*"];

ratk=DSolve[ {-u1''[x]+u1[x]==40x,u1[0]==0,
-u2''[x]+u2[x]==4x,u2'[1]==-1,
u1[1/4]==u2[1/4],u1'[1/4]==u2'[1/4]}, {u1,u2},x]

u1=u1/.ratk[[1]]; u2=u2/.ratk[[1]];

u[x_]:=u1[x] /; x<=1/4
u[x_]:=u2[x] /; x>1/4

fii[e1_,e2_,x_]:=Module[{},
	Which[x<-e1,0,x<0,1+x/e1,x<e2,1-x/e2,True,0]]

f[x_]:=If[x<1/4,40x,4x];

Aloitetaan karkealla diskretoinnilla, vain 2 solmupistetta. Ratkaisu etenee seuraavasti:
valitaan n (tai h), luodaan kantafunktiot, kasataan matriisit ja yhtalon oikea puoli ja
lopuksi ratkaistaan kertoimet ja muodostetaan elementtiratkaisu paloittain lineaarisena
funktiona.

n=2; h=1/n; Do[v[k]=fii[h,h,x-k h],{k,n}]
R=Table[0,{n},{n}]; M=R;
Do[R[[k,k]]=2/h,{k,n-1}]; R[[n,n]]=1/h;
Do[R[[k,k+1]]=-1/h; R[[k+1,k]]=-1/h,{k,n-1}];
Do[M[[k,k]]=2h/3,{k,n-1}]; M[[n,n]]=h/3;
Do[M[[k,k+1]]=h/6; M[[k+1,k]]=h/6,{k,n-1}];
b=Table[0,{n}];
Do[b[[k]]=NIntegrate[v[k] f[x],{x,(k-1)h,(k+1)h}],{k,n-1}];
b[[n]]=-1+NIntegrate[v[n] f[x],{x,(n-1)h,1}];
c=LinearSolve[R+M,b]; uh=Sum[c[[k]] v[k],{k,n}];
Plot[{u[x],uh},{x,0,1},PlotRange->All]

Lasketaan virheet pisteissa 1/4 ja 4/5 seka L2-normissa. Laitetaan tulokset listoihin.

e14={Abs[N[u[1/4]]-uh/.x->1/4]};
e45={Abs[N[u[4/5]]-uh/.x->4/5]};
unorm=NIntegrate[u[x]^2,{x,0,1}];
eL2={Sqrt[
  1-2/unorm*(
    Sum[c[[k]] NIntegrate[u[x] v[k],{x,(k-1)h,(k+1)h}],{k,n-1}]
    +c[[n]] NIntegrate[u[x] v[n],{x,1-h,1}] )
  +c.M.c/unorm]};

Tuplataan n muutamia kertoja ja katsotaan miten virheet kehittyvat.

Do[n=2 n; h=1/n;
  Do[v[k]=fii[h,h,x-k h],{k,n}];
  R=Table[0,{n},{n}]; M=R;
  Do[R[[k,k]]=2/h,{k,n-1}]; R[[n,n]]=1/h;
  Do[R[[k,k+1]]=-1/h,{k,n-1}]; Do[R[[k+1,k]]=-1/h,{k,n-1}];
  Do[M[[k,k]]=2h/3,{k,n-1}]; M[[n,n]]=h/3;
  Do[M[[k,k+1]]=h/6,{k,n-1}]; Do[M[[k+1,k]]=h/6,{k,n-1}];
  b=Table[0,{n}];
  Do[b[[k]]=NIntegrate[v[k] f[x],{x,(k-1)h,(k+1)h}],{k,n-1}];
  b[[n]]=-1+NIntegrate[v[n] f[x],{x,(n-1)h,1}];
  c=LinearSolve[R+M,b]; uh=Sum[c[[k]] v[k],{k,n}];
  e14=AppendTo[e14,Abs[N[u[1/4]]-uh/.x->1/4]];
  e45=AppendTo[e45,Abs[N[u[4/5]]-uh/.x->4/5]];
  eL2=AppendTo[eL2,Sqrt[
    1-2/unorm*(
      Sum[c[[k]] NIntegrate[u[x] v[k],{x,(k-1)h,(k+1)h}],{k,n-1}]
      +c[[n]] NIntegrate[u[x] v[n],{x,1-h,1}] )
    +c.M.c/unorm]],
{4}]

alfa14=Table[Log[2,e14[[k]]/e14[[k+1]] ],{k,4}]

Selvasti O(h^2)

alfa45=Table[Log[2,e45[[k]]/e45[[k+1]] ],{k,4}]

Hmm... ei niin selva, mutta "keskimaarin O(h^2)"

alfaL2=Table[Log[2,eL2[[k]]/eL2[[k+1]] ],{k,4}]

Selvasti O(h^2).

teht. 2

Kokeillaan nyt epatasavalista diskretointia. Otetaan pisteita yhta paljon kuin edellisen tehtavan viimeisella iteraatiolla. Osoittautuu, ettei pisteita kannata uhrata paljon f:n epajatkuvuuskohtaan, koska tarkka ratkaisu on siellakin silea.
Seuraavassa n1 on pisteiden lkm 1/4:n vasemmalla puolella ja n2 oikealla. Listoihin pisteet ja h laitetaan diskretointipisteet ja osavalien pituudet. Pisteiden valinta tapahtuu heuristisella algoritmilla, jota ei tassa tarkemmin kommentoida. Todetaan vain, etta paa-assistentti tarjoaa pullaa ja kahvia sille, joka saa pienemman L2-virheen 32 pisteella.

Clear[h]; n=32;
n1=Round[1/3 n]; n2=n-n1;
h1=1/(2n1(n1+1)); h2=3/(2n2(n2+1));
pisteet=Table[0,{n}]; h=Table[0,{n}]; hh=h1; xx=1/4;
Do[pisteet[[n1-i]]=xx; h[[n1-i]]=hh; xx=xx-hh; hh=hh+h1,
  {i,0,n1-1}];
hh=h2; xx=1;
Do[pisteet[[n-i]]=xx; h[[n-i]]=hh; xx=xx-hh; hh=hh+h2,
  {i,0,n2-1}];

Do[v[i]=fii[h[[i]],h[[i+1]],x-pisteet[[i]]],{i,n-1}];
v[n]=fii[h[[n]],h[[n]],x-1];

R=Table[0,{n},{n}]; M=R;
Do[R[[k,k]]=1/h[[k]]+1/h[[k+1]],{k,n-1}];
   R[[n,n]]=1/h[[n]];
Do[R[[k,k+1]]=-1/h[[k+1]]; R[[k+1,k]]=-1/h[[k+1]],{k,n-1}];
Do[M[[k,k]]=(h[[k]]+h[[k+1]])/3,{k,n-1}];
   M[[n,n]]=h[[n]]/3;
Do[M[[k,k+1]]=h[[k+1]]/6; M[[k+1,k]]=h[[k+1]]/6,{k,n-1}];
b=Table[0,{n}];
Do[b[[k]]=NIntegrate[v[k] f[x],{x,pisteet[[k-1]],pisteet[[k+1]]}],
  {k,2,n-1}];
b[[1]]=NIntegrate[v[1] f[x],{x,0,pisteet[[2]]}];
b[[n]]=-1+NIntegrate[v[n] f[x],{x,1-h[[n]],1}];
c=LinearSolve[R+M,b]; uh=Sum[c[[k]] v[k],{k,n}];

Piirretaan tarkan ja elementtiratkaisun erotus, josta nahdaan onko virhe tasaisesti jakautunut.

Plot[{u[x]-uh},{x,0,1},PlotRange->All]

Virhe L2-normissa:

Sqrt[
 NIntegrate[(u[x]-uh)^2,{x,0,1}] ]

Vertailuarvo tehtavasta 1 (tasavalinen 32 pisteen diskretointi):

eL2[[5]]

teht. 3

psi[h_,x_]:=Module[{},
	Which[x<0,0,x<h,4x (h-x)/h^2,True,0]]

n=2; h=1/n;
Do[v[2k-1]=fii[h,h,x-k h]; v[2k]=psi[h,x-(k-1)h],{k,n}]

Plot[{v[1],v[2],v[3],v[4]},{x,0,1}]

R=Table[0,{2n},{2n}]; M=R;
Do[R[[2k-1,2k-1]]=2/h; R[[2k,2k]]=16/(3h),{k,n}];
   R[[2n-1,2n-1]]=1/h;
Do[R[[2k-1,2k+1]]=-1/h; R[[2k+1,2k-1]]=-1/h,{k,n-1}];
Do[M[[2k-1,2k-1]]=2h/3; M[[2k,2k]]=8h/15, {k,n}];
   M[[2n-1,2n-1]]=h/3;
Do[M[[2k-1,2k]]=h/3; M[[2k,2k-1]]=h/3, {k,n}];
Do[M[[2k-1,2k+1]]=h/6; M[[2k+1,2k-1]]=h/6, {k,n-1}];
Do[M[[2k-1,2k+2]]=h/3; M[[2k+2,2k-1]]=h/3, {k,n-1}];
b=Table[0,{2n}];
Do[b[[2k-1]]=NIntegrate[v[2k-1] f[x],{x,(k-1)h,(k+1)h}],{k,n-1}];
b[[2n-1]]=-1+NIntegrate[v[2n-1] f[x],{x,(n-1)h,1}];
Do[b[[2k]]=NIntegrate[v[2k] f[x],{x,(k-1)h,k h}],{k,n}];
c=LinearSolve[R+M,b]; uh=Sum[c[[k]] v[k],{k,2n}];
Plot[{u[x],uh},{x,0,1},PlotRange->All]

Erinomainen ratkaisu neljalla kantafunktiolla! Katsotaan nyt virheen kehittymista kuten
tehtavassa 1.

eL2={Sqrt[
  1-2/unorm*(
    Sum[c[[2k-1]] NIntegrate[u[x] v[2k-1],
    	{x,(k-1)h,(k+1)h}],{k,n-1}]
    +c[[2n-1]] NIntegrate[u[x] v[2n-1],{x,1-h,1}]
    +Sum[c[[2k]] NIntegrate[u[x] v[2k],{x,(k-1)h,k h}],{k,n}] )
  +c.M.c/unorm]};

Do[n=2*n; h=1/n;
  Do[v[2k-1]=fii[h,h,x-k h]; v[2k]=psi[h,x-(k-1)h],{k,n}];
  R=Table[0,{2n},{2n}]; M=R;
  Do[R[[2k-1,2k-1]]=2/h; R[[2k,2k]]=16/(3h),{k,n}];
     R[[2n-1,2n-1]]=1/h;
  Do[R[[2k-1,2k+1]]=-1/h; R[[2k+1,2k-1]]=-1/h,{k,n-1}];
  Do[M[[2k-1,2k-1]]=2h/3; M[[2k,2k]]=8h/15, {k,n}];
     M[[2n-1,2n-1]]=h/3;
  Do[M[[2k-1,2k]]=h/3; M[[2k,2k-1]]=h/3, {k,n}];
  Do[M[[2k-1,2k+1]]=h/6; M[[2k+1,2k-1]]=h/6, {k,n-1}];
  Do[M[[2k-1,2k+2]]=h/3; M[[2k+2,2k-1]]=h/3, {k,n-1}];
  b=Table[0,{2n}];
  Do[b[[2k-1]]=NIntegrate[v[2k-1] f[x],{x,(k-1)h,(k+1)h}],{k,n-1}];
  b[[2n-1]]=-1+NIntegrate[v[2n-1] f[x],{x,(n-1)h,1}];
  Do[b[[2k]]=NIntegrate[v[2k] f[x],{x,(k-1)h,k h}],{k,n}];
  c=LinearSolve[R+M,b]; uh=Sum[c[[k]] v[k],{k,2n}];
  eL2=AppendTo[eL2,Sqrt[
    1-2/unorm*(
      Sum[c[[2k-1]] NIntegrate[u[x] v[2k-1],
      	{x,(k-1)h,(k+1)h}],{k,n-1}]
      +c[[2n-1]] NIntegrate[u[x] v[2n-1],{x,1-h,1}]
      +Sum[c[[2k]] NIntegrate[u[x] v[2k],{x,(k-1)h,k h}],{k,n}] )
   +c.M.c/unorm]],
{4}]

alfaL2=Table[Log[2,eL2[[k]]/eL2[[k+1]] ],{k,4}]

Virhe on O(h^3) joten toisen asteen polynomit tuovat merkittavan parannuksen.
Samaan tarkkuuteen kuin tehtavassa 1 paastaan huomattavasti pienemmalla joukolla
kantafunktioita.

teht. 4

Tutkitaan 3 pienimman ominaisarvon kehittymista h:n funktiona seka lineaarisilla etta
toisen asteen kantafunktioilla.

n=2; virheet1={}; virheet2={};
Do[n=2n; h=N[1/n];
  R1=Table[0,{n},{n}]; M1=R1;
  Do[R1[[k,k]]=2/h,{k,n-1}]; R1[[n,n]]=1/h;
  Do[R1[[k,k+1]]=-1/h; R1[[k+1,k]]=-1/h,{k,n-1}];
  Do[M1[[k,k]]=2h/3,{k,n-1}]; M1[[n,n]]=h/3;
  Do[M1[[k,k+1]]=h/6; M1[[k+1,k]]=h/6,{k,n-1}];

  R2=Table[0,{2n},{2n}]; M2=R2;
  Do[R2[[2k-1,2k-1]]=2/h; R2[[2k,2k]]=16/(3h),{k,n}];
     R2[[2n-1,2n-1]]=1/h;
  Do[R2[[2k-1,2k+1]]=-1/h; R2[[2k+1,2k-1]]=-1/h,{k,n-1}];
  Do[M2[[2k-1,2k-1]]=2h/3; M2[[2k,2k]]=8h/15, {k,n}];
     M2[[2n-1,2n-1]]=h/3;
  Do[M2[[2k-1,2k]]=h/3; M2[[2k,2k-1]]=h/3, {k,n}];
  Do[M2[[2k-1,2k+1]]=h/6; M2[[2k+1,2k-1]]=h/6, {k,n-1}];
  Do[M2[[2k-1,2k+2]]=h/3; M2[[2k+2,2k-1]]=h/3, {k,n-1}];
  
  lista1=Sort[Eigenvalues[Inverse[M1].R1]];
  lista2=Sort[Eigenvalues[Inverse[M2].R2]];
  AppendTo[virheet1,
    Table[lista1[[i]]-N[Pi^2(i-1/2)^2],{i,3}]];
  AppendTo[virheet2,
    Table[lista2[[i]]-N[Pi^2(i-1/2)^2],{i,3}]],
{4}]

MatrixForm[virheet1]

MatrixForm[virheet2]

Yllaolevissa taulukoissa on ominaisarvojen virheet ensin h:n arvolla 1/4 (1. rivit) jonka
jalkeen h puolittuu aina seuraavalle riville siirryttaessa. Nayttaisi, etta lineaarisilla virhe
puolittuu ja kvadraattisilla putoaa neljasosaan. Allaolevat taulukot vahvistavat asian.

suhteet1=Table[virheet1[[i,j]]/virheet1[[i+1,j]],{i,3},{j,3}];
suhteet2=Table[virheet2[[i,j]]/virheet2[[i+1,j]],{i,3},{j,3}];

MatrixForm[suhteet1]

MatrixForm[suhteet2]