Eksamen i fag TKT4140 Numeriske berekningar m/datalab

Ved å derivere likning \eqref{eq:ODEsol} får me \( y' = \frac{-1}{x^2} \) og \( y''=\frac{2}{x^3} \) eg ved innsetjing i likning \eqref{eq:ODE} og \eqref{eq:BCs} får me: $$ \begin{align} \frac{2}{x^3} + 2 \left(\frac{1}{x} \frac{-1}{x^2}\right) = 0 \quad y(x_{\min}) = y &= \frac{1}{1}=1 \quad y(x_{\max}) = \frac{1}{2} = 0.5 \label{eq:ODEsol_show} \end{align} $$

Me brukar at \( y_0 = y \) som vidare gir: $$ \begin{align*} y_0^\prime &= y_1 \label{} \\ y_1^\prime &= -2\left(y_0 \cdot y_1 \right) \label{} \end{align*} $$

Merk at programmet nyttar funksjonen rk4 frå modulen ODEschemes til å løyse initialverdiproblemet for kvar startverdi me tippar. Denne funksjonen kan du rekne som kjend og det er ikkje naudsynt å greia ut om denne.

Det som manglar for at me kan løyse dette som eit initialverdiproblem er \( y'(x_{\min}) \) som her blir verdien for s. Frå den analytiske løysinga har me at \( y'=\frac{-1}{x^2} \rightarrow y'(x_{\min}) = \frac{-1}{1^2} = -1 \). Som er verdien av s som me ventar at programmet vil avslutta med.

# Eksamen2016/oppg1c.py


import numpy as np
import matplotlib.pylab as plt
from ODEschemes import euler, heun, rk4


def func(y, x):
    y0 = y[0]
    y1 = y[1]
    
    return np.array([y1, -2*y0*y1])
    

x_0, x_end = 1, 2
y_0, y_end = 1, 0.5
N = 50

X = np.linspace(x_0, x_end, N +1)
Y_analytic = 1/X

s0 = 0
Y_0 =[y_0, s0]

Y_num = euler(func, Y_0, X)
Y0 = Y_num[:, 0]

phi0 = y_end - Y0[-1]

s1 = -0.5

itMax = 10
tol = 10**(-10)
for n in range(itMax):
    Y_0 =[y_0, s1]

    Y_num = euler(func, Y_0, X)
    Y0 = Y_num[:, 0]

    phi1 = y_end - Y0[-1]
    
    s2 = s1 - phi1*(s1 - s0)/(phi1 - phi0)
    
    s0, s1, phi0 = s1, s2, phi1

    if abs(s1 - s0) < tol:
        break   

print n, s1, phi1

plt.figur()
plt.plot(X, Y0)
plt.plot(X, Y_analytic, "k--")
plt.show()

Bruk at $$ \begin{align} F\left(y_{i-1}, y_{i}, y_{i+1}\right)_{m+1} = F_m &+ \frac{\partial F_m}{\partial y_{i-1}^m}\;\delta y_{i-1} + \frac{\partial F_m}{\partial y_{i}^m}\;\delta y_{i} + \frac{\partial F_m}{\partial y_{i+1}^m}\;\delta y_{i+1} \label{eq:newton} \end{align} $$ der \( F_m = F\left(y_{i-1}^m, y_{i}^m, y_{i+1}^m\right) \) er det ikkje-lineære leddet, og $$ \begin{align*} \delta y_{i-1} = y_{i-1}^{m+1} - y_{i-1}^{m}, \qquad \delta y_{i} = y_{i}^{m+1} - y_{i}^{m}, \qquad \delta y_{i+1} = y_{i+1}^{m+1} - y_{i+1}^{m} \end{align*} $$ der me har innført sub-iterasjonsindeksen \( m \).

Me brukar sentraldifferansar slik at \( y^{\prime\prime} \approx \frac{y_{i-1} - 2 \, y_i + y_{i+1}}{{\Delta x}^2} \) og \( y^{\prime} \approx \frac{y_{i+1}-y_{i-1}}{2 \, \Delta x} \), som insatt i likning \eqref{eq:ODE} gir: $$ \begin{align} \frac{y_{i-1} - 2 \, y_i + y_{i+1}}{{\Delta x}^2} + 2 \, y_i \frac{y_{i+1}-y_{i-1}}{2 \, \Delta x} = 0 \label{_auto1}\\ \rightarrow y_{i-1} - 2 \, y_i + y_{i+1} + \Delta x \, y_i \left( y_{i+1}-y_{i-1} \right) = 0 \label{eq:discretized} \end{align} $$

Det ikkje-lineære leddet blir då \( F_{m+1} = y_i^{m+1} \left( y_{i+1}^{m+1}-y_{i-1}^{m+1} \right) \), og me bruk av likning \eqref{eq:newton} får me: $$ \begin{align} F_{m+1} \approx y_i^{m} \left( y_{i+1}^{m}-y_{i-1}^{m} \right) - \left( y_{i-1}^{m+1}-y_{i-1}^m \right) y_{i}^m + \left( y_{i}^{m+1}-y_{i}^m\right) \left( y_{i+1}^{m}-y_{i-1}^{m}\right) + \left( y_{i+1}^{m+1}-y_{i+1}^m\right) y_{i}^m \nonumber \\ = - y_i^m \, y_{i-1}^{m+1} + \left( y_{i+1}^m - y_{i-1}^m \right) \, y_i^{m+1} + y_i^m \, y_{i+1}^{m+1} - y_i^m \left( y_{i+1}^m - y_{i-1}^m\right) \label{_auto2}\\ = -\alpha_i \, y_{i-1}^{m+1} + \beta_i \, y_i^{m+1} + \alpha_i \, y_{i+1}^{m+1} - \alpha_i \cdot \beta_i \nonumber \end{align} $$ Der me har brukt at \( \alpha_i = y_i^m \) og \( \beta_i = y_{i+1}^m - y_{i-1}^m \), som dessutan gir at \( \alpha_i \cdot \beta_i = F_m \). Vidare innsetjing i \eqref{eq:discretized}: $$ \begin{align} \left( 1 - \Delta x \, \alpha_i \right) y_{i-1}^{m+1} + \left(\Delta x \, \beta_i -2 \right) y_i^{m+1} + \left( 1 + \Delta x \, \alpha_i\right) y_{i+1}^{m+1} = \Delta x \, \alpha_i \cdot \beta_i \label{_auto3} \end{align} $$

For \( i=1 \) får me: $$ \begin{align} \left(\Delta x \, \beta_1 -2 \right) y_1^{m+1} + \left( 1 + \Delta x \, \alpha_1\right) y_{2}^{m+1} = \Delta x \, \alpha_1 \cdot \beta_1 - \left( 1 - \Delta x \, \alpha_1 \right) y(x_{\min}) \label{_auto4} \end{align} $$

og for \( i=N-1 \): $$ \begin{align} \left( 1 - \Delta x \, \alpha_{N-1} \right) y_{N-2}^{m+1} + \left(\Delta x \, \beta_{N-1} -2 \right) y_{N-1}^{m+1} = \Delta x \, \alpha_{N-1} \cdot \beta_{N-1} - \left( 1 + \Delta x \, \alpha_{N-1}\right) y(x_{\max}) \label{_auto5} \end{align} $$

med \( C = \frac{a_0 \Delta t}{\Delta x} \)

Skjemaet er \( O(\Delta t) \) og \( O(\Delta x^2) \). Forover-differansar gir første orden i tid og sentral-differansar gir andre orden i rom.

Skjemaet kan skrivast om til: $$ \begin{equation} u_j^{n+1}=\frac{1 - C}{2} u_{j + 1}^n + \frac{1 + C}{2} u_{j - 1}^n \label{eq:lax_pk} \end{equation} $$

Ser at koeffisientane summeres opp til 1; \( \frac{1 - C}{2} + \frac{1 + C}{2} = 1 \), og at andre ledd på høgre side alltid er positivt medan det første leddet er positiv dersom \( C \leq 1 \), som gir eit tilstrekkeleg stabilitets-krav for Lax-Friedrich. Dersom \( C < 0 \) har me at det andre leddet er positivt dersom \( C \geq -1 \).

Skjemaet kan skrivast om til: $$ \begin{equation} u_j^{n+1}=\frac{1}{2} \left(u_{j + 1}^n + u_{j - 1}^n\right) - \frac{C}{2} \left(u_{j + 1}^n - u_{j - 1}^n\right) \label{eq:lax_Von} \end{equation} $$ $$ \begin{equation} G^{n+1} \cdot e^{i \beta y_j}=\frac{1}{2} \left(G^{n} \cdot e^{i \beta y_{j+1}} + G^{n} \cdot e^{i \beta y_{j-1}}\right) - \frac{C}{2} \left(G^{n}\cdot e^{i \beta y_{j+1}} - G^{n}\cdot e^{i \beta y_{j-1}}\right) \label{eq:lax_Von1} \end{equation} $$

Ved å introdusera \( y_{j \pm 1} = y_j \pm h \) and \( \delta = \beta \cdot h \) then dividing by \( G^{n} \cdot e^{i y_j} \) får me: $$ \begin{equation} G =\frac{1}{2} \left(e^{i \delta} + e^{- i \delta}\right) - \frac{C}{2} \left(e^{i \delta} - e^{-i \delta}\right) \label{eq:lax_Von2} \end{equation} $$

med hjelp av følgjande relasjonar: $$ \begin{equation*} e^{i \delta} = cos \delta + i sin \delta \\ e^{-i \delta} = cos \delta - i sin \delta \\ e^{i \delta} + e^{- i \delta} = 2cos \delta \\ e^{i \delta} - e^{- i \delta} = 2 i sin \delta \nonumber \end{equation*} $$

får me: $$ \begin{equation} G = cos \delta - i C sin \label{eq:lax_Von3} \end{equation} $$

og $$ \begin{equation} \lvert G \rvert = \sqrt{cos^2 \delta + C^2 sin^2 \delta} \label{eq:lax_Von4} \end{equation} $$

kor me har brukt at modulusen til eit komplekst tal \( z = a + i b \), er \( \lvert z\rvert = \sqrt{a^2 + b^2} \)

Me har at \( \lvert G \rvert \leq 1 \) (\( \frac{\lvert G \rvert}{\lvert G_a \rvert } \leq 1 \) med \( \lvert G_a \rvert = 1 \)). Som vidare gir at \( -1 \leq C \leq 1 \).

Von Neumann gir eit naudsynt og dermed også tilstrekkeleg kriteria. PK gir eit tilstrekkeleg kriteria, som dermed kan vera meir strengt enn Von-Neumann. Her ser med derimot at dei gir like kriteria.

Frå førre oppgåve har me at $$ \begin{equation*} \lvert G \rvert = \sqrt{cos^2 \delta + C^2 sin^2 \delta} \nonumber \end{equation*} $$

som gir: $$ \begin{equation*} \delta = 0 \rightarrow \lvert G \rvert = 1 \\ \delta = \frac{\pi}{2} \; rad = 90^o \rightarrow \lvert G \rvert = C \\ \delta = \pi \; rad = 180^o \rightarrow \lvert G \rvert = 1 \nonumber \end{equation*} $$

som samsvarar med den blå linja. Lax-Wendroff er eit andre ordens skjema i både tid og rom og er dermed mykje mindre diffusivt en dei andre to (for relativt små verdier av \( \delta \)), og samsvarer derfor med den raue linja.

Frå figur 2 har me at \( \lvert G \rvert_{LW} > \lvert G \rvert_{FTBS} > \lvert G \rvert_{LF} \). Dette gir at tap av amplitude for Lax Friedrich (LF) er større enn hos FTBS, som igjen er større enn hos Lax Wendroff (LW). I tilegg veit me at Lax Wendroff som er eit andre ordens skjema kan gi oskilasjonar i nærleiken av løysingar med store gradientar (lilla linje).

I likning \eqref{eq:oppg3_disk} er den eindimensjonale diffusjonslikninga: $$ \begin{equation} \label{eq:5211} \frac{\partial u}{\partial t}=\alpha\frac{\partial^2u}{\partial x^2} \end{equation} $$

diskretisert. Eksemplar på problem er varmeoverføring, oppstart av rør-straumning og straumning i porøse media.

skjemaet i likning \eqref{eq:oppg3_disk} er eksplisitt for \( \theta=0 \), og elles implisitt. Ein seier at skjemaet er eksplisitt dersom \( u \) kan skrivast på ein formel der høgre side berre inneheld kjende størrelsar, altså dersom \( u_j^{n+1} \) kan uttrykkjast som ein funksjon av \( u^n \). Dersom det er fleire ukjente på høgre side (\( u_{j-1}^{n+1} \) og \( u_{j+1}^{n+1} \) til dømes), er skjemaet implisitt, og må løysast som eit system av likningar.

Med \( \theta=1/2 \) blir både \( \frac{\partial u}{\partial t} \) og \( \frac{\partial^2u}{\partial x^2} \) diskretisert med sentraldifferansar, og me får Crank-Nicholson skjemaet som er av 2. orden.

Uttrykkjet i \eqref{eq:order} forventast å nærme seg skjemaet si teoretiske orden; 2 for heuns metode som er av 2. orden.