[ModSim] H9 (LMM en shooting)

1994 days ago by Toon.Baeyens

Hoofdstuk 9

Met behulp van onderstaande techniek kunnen we lineaire meerstapsmethoden bepalen van een gewenste orde.

Stel dat we het probleem $y'(t) = f(t)$ wensen op te lossen. Hiervoor kunnen we een LMM ontwikkelen:

$$ y_{n+1} + a_1 y_{n} + a_2 y_{n-1} = h\left(b_0 f(t_n+h) + b_1 f(t_n) + b_2 f(t_n -h)\right) $$

We zijn op zoek naar coëfficienten $a_1$, $a_2$, $b_0$, $b_1$ en $b_2$ waarvoor de taylorontwikkeling naar $h$ van 

$$ y(t_n+h) + a_1 y(t_n) + a_2 y(t_n-h) - h\left(b_0 y'(t_n+h) + b_1 y'(t_n) + b_2 y'(t_n -h)\right) $$

Deze symbolische uitwerking kunnen we in Sage uitvoeren. We gaan op zoek naar de expliciete 2 staps Adams LMM methode. Dit is de methode van hoogste orde waarvoor $a_2 = 0$.

# We veronderstellen tn = 0 h = var('h') # een variabele y = function('y')(h) # een symbolische functie die van h afhangt f = diff(y, h) a = [1, var('a1'), 0] b = [0, var('b1'), var('b2')] # bij expliciete methoden is b0 = 0 y 
       

                                
                            

                                

Als eerste stellen we de fout die onze methode geeft op:

error = sum(a[i]*y(h-i*h) for i in range(3)) - h*sum(b[i]*f(h-i*h) for i in range(3)) error 
       

                                
                            

                                

Vervolgens nemen we hier de Taylorreeks naar $h$ van:

error_series = error.series(h, 7) error_series 
       

                                
                            

                                

Als laatste eisen we dat de eerste paar coefficiënten in deze ontwikkeling $0$ zijn:

# We proberen verschillende waarden voor n uit. Bij n = 3 zien we dat dit stelsel exact 1 oplossing heeft n = 3 solutions = solve([error_series.coefficient(h, i) == 0 for i in range(n)], [a[1], b[1], b[2]], solution_dict=True) solutions 
       

                                
                            

                                

We voeren nog een laatste controle uit om te zien dat deze Adams-Bashforth methode orde 2 heeft:

error_series.subs(solutions[0]) 
       

                                
                            

                                

Probeer deze zelfde techniek ook eens voor de 2 staps Adams-Moulton methode. Dit is een impliciete methode met $a_2 = 0$.

Bepaal daarnaast ook de 3 staps Adams-Bashforth methode. Expliciet met $a_2 = a_3 = 0$.

Controlleer jouw resultaten aan de hand van https://en.wikipedia.org/wiki/Linear_multistep_method#Families_of_multistep_methods.

 
       
 
       

Op dezelfde wijze kunnen we ook methoden opstellen voor het probleem $y''(t) = f(t)$.

Hierbij neemt de methode dan de vorm aan (met $a_0 = 1$):

$$ \sum_{i=0}^s a_i y_{n+1-i} = h^2\sum_{i=0}^s b_i f(t_{n}+h-ih) $$

Stel de Numerovmethode op. Deze is impliciet met maximale orde. Gebruik hiervoor dezelfde stappen als daarnet, merk op dat nu $f(t) = y''(t)$ en de error krijgt een factor $h^2$.

 
       
 
       

De Poissonvergelijking

Beschouw de volgende vergelijking

$$ y''(t) = f(t, y) \text{.}$$

Wanneer $f(t, y) = f(t)$ wordt dit de Poissonvergelijking genoemd.

We zullen een randwaarde probleem bij deze vergelijking oplossen. Als eerste maken we een discretisatie in $n$ intervallen, dus de punten $t_0, t_1, \dots, t_n$.

Met behulp van de Numerov methode kunnen we een verband leggen tussen de onbekenden $y_1, y_2, \dots, y_{n-1}$:

$$y_{k+1} - 2y_k + y_{k-1} = \frac{h^2}{12} \left(f(t_{k+1}, y_{k+1}) + 10 f(t_{k}, y_k) + f(t_{k-1}, y_{k-1})\right)\qquad\text{voor $k=1,2,\dots, n-1$}$$

De oplossing van dit stelsel ($n-1$ onbekenden, $n-1$ vergelijkingen) benadert de oplossing voor het randwaardenprobleem.

 

Stel $f(t, y) = \cos(t)$, $t_0 = 0$ en $t_n = \pi$, als randvoorwaarden stellen we $y_0 = y(t_0) = 0$ en $y_n = y(t_n) = 0$. Vergelijk oplossingen voor verschillende waarden van $n$. De exacte oplossing is $y(t) = 1 -\frac{2 \, t}{\pi} - \cos\left(t\right)$.

Doe dezelfde oefening ook met $f(t, y) = \sin(4t)y$ en $y(0) = y(\pi) = 1$.

def poisson(f, a, b, ya, yb, n): h = (b-a)/n t = [RDF(a + i*h) for i in range(0,n+1)] y = None # TODO f = None # TODO equations = None # TODO solution = solve(equations, y[1:-1], solution_dict=True)[0] return t, [ya] + [RDF(solution[y[i]]) for i in range(1,n)] + [yb] 
       
 
       

Het probleem

$$ y''(t) = \sin(4t)y \qquad \text{met $y(0) = y(\pi) = 1$}$$

kunnen we ook met behulp van shooting oplossen. Hiervoor is het nodig om het probleem te vertalen naar een stelsel eerste orde differentiaalvergelijkingen:

$$\begin{pmatrix}y(t)\\y'(t)\end{pmatrix}' = \begin{pmatrix}y'(t)\\ f(t,y)\end{pmatrix} = \begin{pmatrix}y'(t)\\ \sin(4t)y\end{pmatrix}$$

Dan kan men gebruik maken van de volgende functie om het beginwaarde probleem op te lossen.

def ode_solve(f, t0, t1, y0, n): h = RDF((t1-t0)/n) def F(*args): fy = f(*args) return RDF(fy) if fy in RR else vector(RDF, fy) t = [RDF(t0)] y = [RDF(y0) if y0 in RR else vector(RDF, y0)] for i in range(n): ti = t[-1] yi = y[-1] k1 = h*F(ti, yi) k2 = h*F(ti+h/2, yi+k1/2) k3 = h*F(ti+h/2, yi+k2/2) k4 = h*F(ti+h, yi+k3) y.append(yi + (k1 + 2*(k2+k3) +k4)/6) t.append(ti + h) return t, y 
       

Als beginvoorwaarde neemt men $y(0) = 1$ en voor $y'(0)$ probeert men verschillende waarden tot $y(\pi) = 0$ blijkt.

 
       

Experimenteer met verschillende waarden voor $n$.

 
       
 
       

Schrödingervergelijking

Deze shooting techniek kan men in een gewijzigde vorm ook toepassen op de Schrödingervergelijking:

$$y''(t) = (V(t) - E) y(t)\text{.}$$

Voor een onbekende $E$, een gegeven $V$ en gegeven randvoorwaarden.

 

Kies $V(t) = 2\cos(2t)$ en $y(-\frac{\pi}{2}) = y(\frac{\pi}{2}) = 0$.

Pas shooting toe:

  • Gok een waarde voor $E$
  • Los het beginvoorwaarde probleem op: $y''(t) = (V(t) - E)y(t)$ met $y(-\frac{\pi}{2}) = 0$ en $y'(-\frac{\pi}{2}) = 1$. Gebruik 'ode_solve'.
  • Bereken hoe ver $y(\frac{\pi}{2})$ van $0$ afwijkt en corrigeer $E$.
  • Herhaal tot $y(\frac{\pi}{2}) \approx 0$