Harjoitus 11 -- Ratkaisut

1

Ratkaistava yhtälö

eq=y''[x]+(k^2)y[x]==0;

Ratkaisuvälin pituus voi olla mitä tahansa (oltava numeerinen, jos käytetään NDSolvea)

L=1.0;

Aaltoluku k=((2π)/λ), missä λ=n(L/(2 )), n kokonaisluku (tämä on tarkka ratkaisu, oletetaan, että tätä ei tiedetä). "Arvataan"

k1=8.0;

Ratkaisuprosessi: ratkaistaan yhtälö alueessa [0,L/2] reunaehdoilla y(0)=0 ja y'(0)=1, ja alueessa [L/2,L] reunaehdoilla y(l)=0 ja y'(L)=±1 (ensin +1)

y0[x_]=y[x]/.Flatten[NDSolve[{eq/.k->k1,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y∞[x_]=y[x]/.Flatten[NDSolve[{eq/.k->k1,y[L]==0.0,y'[L]==1.0},y[x],{x,L/2,L}]];
Show[Plot[y0[x],{x,0,L/2}],Plot[y∞[x],{x,L/2,L}],PlotRange->All]

Ratkaisun hyvyyttä voi arvioida erotuksen (Subscript[y, 0])(L/2)-(Subscript[y, ∞])(L/2) avulla. Huom. nyt funktioiden derivaatoilla on aina sama arvo pisteessä L/2.

y0[L/2]-y∞[L/2]

Toinen alkuarvaus välinpuolitusmenetelmää varten:

k2=5.0;
y0[x_]=y[x]/.Flatten[NDSolve[{eq/.k->k2,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y∞[x_]=y[x]/.Flatten[NDSolve[{eq/.k->k2,y[L]==0.0,y'[L]==1.0},y[x],{x,L/2,L}]];
Show[Plot[y0[x],{x,0,L/2}],Plot[y∞[x],{x,L/2,L}],PlotRange->All]
y0[L/2]-y∞[L/2]

Kirjoitetaan erotuksen lasku omaksi funktiokseen:

Clear[Δ]

Δ[κ_]:=Module[{},
q0=Flatten[NDSolve[{eq/.k->κ,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y0[x_]=y[x]/.q0;
q∞=Flatten[NDSolve[{eq/.k->κ,y[L]==0.0,y'[L]==1.0},y[x],{x,L/2,L}]];
y∞[x_]=y[x]/.q∞;
y0[L/2]-y∞[L/2]]

Δ[8.0]
Δ[5.0]

Etsitään funktion Δ nollakohta (yhtälön ominaisarvo) tarkkuudella 1/1000000:

k1=8.0;k2=5.0;
While[Abs[k1-k2]>(10^(-6)),
k0=((k1+k2)/2);
If[Δ[k0]Δ[k2]<0,k1=k0,k2=k0]];
k0

Kirjoitetaan välinpuolitusmenetelmä omaksi funktiokseen (ks. harj. 7):

bisection[f_,a_,b_]:=Module[{z1,z2},
z1=a;
z2=b;
While[Abs[z1-z2]>(10^(-6)),
z=((z1+z2)/2);
If[f[z]f[z2]<0,z1=z,z2=z]];
z];

bisection[Δ,5.0,8.0]
bisection[Δ,8.0,14.0]

Vastaava ominaisfunktio. Normitetaan ratkaisu niin, että sen suurin arvo on 1. Tätä koodia voisi vähän optimoida.

Remove[plotsolution]

plotsolution[κ_,h_]:=Module[{},
q0=Flatten[NDSolve[{eq/.k->κ,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y0[x_]=y[x]/.q0;
q∞=Flatten[NDSolve[{eq/.k->κ,y[L]==0.0,y'[L]==1.0},y[x],{x,L/2,L}]];
y∞[x_]=y[x]/.q∞;
η[x_]:=y0[x]/;0≤x≤L/2;
η[x_]:=y∞[x]/;L/2<x≤L;
ηmax=Part[NMaximize[η[x],x],1];
ξ[x_]=η[x]/ηmax;
Show[Plot[ξ[x],{x,0,L},PlotStyle->Hue[h],PlotRange->All]]
]

plotsolution[k0,1]

Lähdetään etsimään lisää ratkaisuja. Edetään ykkösen pituisin askelin kohti suurempia k:n arvoja. Jos kahdella peräkkäisellä k:lla lasketut Δ:n arvot ovat erimerkkiset, käytetään välinpuolitusmenetelmää.

k1=1.0;
i=1;
While[i≤10,
If[Δ[k1]Δ[k1+1.0]<0,κ[i]=bisection[Δ,k1,k1+1.0];i++];
k1+=1.0;
]

Kymmenen ensimmäistä ominaisarvoa (ne ovat 2π:n välein):

?κ

Do[g[i]=plotsolution[κ[i],i/4];,{i,4}]
Show[Table[g[i],{i,4}]]

Etsitään seuraavaksi reunaehtoa y'(L)=-1 vastaavat ratkaisut. Nyt käy niin, että funktioiden (Subscript[y, 0]) ja (Subscript[y, ∞]) arvot pisteessä L/2 ovat aina samat, joten riittää vertailla derivaattojen arvoja.

Σ[κ_]:=Module[{},
q0=Flatten[NDSolve[{eq/.k->κ,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y0[x_]=y[x]/.q0;
q∞=Flatten[NDSolve[{eq/.k->κ,y[L]==0.0,y'[L]==-1.0},y[x],{x,L/2,L}]];
y∞[x_]=y[x]/.q∞;
y0'[L/2]-y∞'[L/2]]

k1=1.0;
i=1;
While[i≤10,
If[Σ[k1]Σ[k1+1.0]<0,λ[i]=bisection[Σ,k1,k1+1.0];i++];
k1+=1.0;
]

Aaltoluvun arvot ovat taas 2π:n välein, mutta alkaen tällä kertaa π:stä. (L:n matkalle mahtuu n+1/2 aallonpituutta)

?λ

plotsolution2[κ_,h_]:=Module[{},
q0=Flatten[NDSolve[{eq/.k->κ,y[0]==0.0,y'[0]==1.0},y[x],{x,0.0,L/2}]];
y0[x_]=y[x]/.q0;
q∞=Flatten[NDSolve[{eq/.k->κ,y[L]==0.0,y'[L]==-1.0},y[x],{x,L/2,L}]];
y∞[x_]=y[x]/.q∞;
η[x_]:=y0[x]/;0≤x≤L/2;
η[x_]:=y∞[x]/;L/2<x≤L;
ηmax=Part[NMaximize[η[x],x],1];
ξ[x_]=η[x]/ηmax;
Show[Plot[ξ[x],{x,0,L},PlotStyle->Hue[h],PlotRange->All]]
]

Do[g[i]=plotsolution2[λ[i],i/4];,{i,4}]
Show[Table[g[i],{i,4}]]

Sama differenssimenetelmällä

Taulukko x:n arvoista, joilla ratkaisua etsitään: (luodaan pariton määrä pisteitä, jotta niistä voidaan määrätä keskimmäinen, jossa ratkaisun hyvyyttä voidaan tutkia.

Clear[x,y]

n=101;x[1]=0.0;x[n]=L;h=((x[n]-x[1])/(n-1));
Do[x[i]=x[1]+(i-1)h,{i,2,n-1}]

Differenssiyhtälö:
	(((Subscript[y, i+1])-2(Subscript[y, i])+(Subscript[y, i-1]))/(h^2))+(k^2)(Subscript[y, i])=0,
josta
	(Subscript[y, i±1])=(2-(h^2)(k^2))(Subscript[y, i])-(Subscript[y, i∓1])

Ratkaisu yhdelle k:n arvolle, reunaehdolla y'(L)=1

y[1]=0.0;dy0=1.0;y[2]=y[1]+dy0 h;
y[n]=0.0;dyL=1.0;y[n-1]=y[n]-dyL h;

k=7.0;

Funktion asrvojen erotus välin keskipisteessä, l≡((n+1)/2), on
	(Subscript[y, 0])((Subscript[x, l]))-(Subscript[y, ∞])((Subscript[x, l]))=(2-(h^2)(k^2))(Subscript[y, l-1])-(Subscript[y, l-2])-(2-(h^2)(k^2))(Subscript[y, l+1])+(Subscript[y, l+2])
	=(2-(h^2)(k^2))((Subscript[y, l-1])-(Subscript[y, l+1]))-((Subscript[y, l-2])-(Subscript[y, l+2]))

l=((n+1)/2);

Do[
y[i]=(2-(h^2)(k^2))y[i-1]-y[i-2];
y[n-i+1]=(2-(h^2)(k^2))y[n-i+2]-y[n-i+3];,
{i,3,((n-1)/2)+1}]
(*funktion arvojen erotus keskipisteessä*)
(2-(h^2)(k^2))(y[l-1]-y[l+1])-(y[l-2]-y[l+2])

ListPlot[Table[{x[i],y[i]},{i,1,n}]]

Kerätään käskyt funktioiksi. Näissä käytetään koko ajan samaa resoluutiota ja samoja reunaehtoja.

Δd[k_]:=Module[{},
Do[
y[i]=(2-(h^2)(k^2))y[i-1]-y[i-2];
y[n-i+1]=(2-(h^2)(k^2))y[n-i+2]-y[n-i+3];,
{i,3,((n-1)/2)+1}];
(2-(h^2)(k^2))(y[l-1]-y[l+1])-(y[l-2]-y[l+2])
];

plotd[k_]:=Module[{},
Do[
y[i]=(2-(h^2)(k^2))y[i-1]-y[i-2];
y[n-i+1]=(2-(h^2)(k^2))y[n-i+2]-y[n-i+3];,
{i,3,((n-1)/2)+1}];
ListPlot[Table[{x[i],y[i]},{i,n}],Joined->True]
]

Δd[12.0]
Δd[15.0]

Välinpuolitusmenetelmä:

plotd[bisection[Δd,12.0,15.0]]