Forelesning del 1

Fredrik Eikeland Fossan

Feb 19, 2016

Example 1: bruk av MMS på Heun's skjemat

La oss bruke MMS på følgende differensialligning: $$ \begin{align} \label{eq:simple_linear_ode} u'' + u^2 = g \end{align} $$

Denne ligningen kan omgjøres til et system av førsteordens differensialligninger (ved bruk av $ u_0 = u $): $$ \begin{align*} u_0'&= u_1& \\ u_1'&= g - u_0^2& \end{align*} $$

Predictor: $$ \begin{align*} u_0^p&=(u_0)_n+h\cdot (u_1)_n& \\ u_1^p&=(u_1)_n+h\cdot \left[g(t_n) - (u_0)_n^2\right]& \end{align*} $$

Corrector: $$ \begin{align*} (u_0)_{n+1}&=(u_0)_n+0.5h\cdot [(u_1)_n+(u_1^p)]&\\ (u_1)_{n+1}&=(u_1)_n+0.5h\cdot \left[\left(g(t_n) - \left((u_0)_n\right)^2\right) + \left(g(t_n+ h) - \left(u_0^p\right)^2\right)\right]& \end{align*} $$

eller ved å innføre vektorer: $$ \begin{equation} \label{eq:3102} \boldsymbol{U}' = \begin{bmatrix} u_0' \\ u_1' \end{bmatrix} = \begin{bmatrix} u_1 \\ g - u_0^2 \end{bmatrix} \end{equation} $$

predictor: $$ \begin{equation} \boldsymbol{U}_{n+1}^p = \begin{bmatrix} u_0^p \\ u_1^p \end{bmatrix} = \boldsymbol{U}_n + h \cdot \boldsymbol{U}_{n}' \end{equation} $$

corrector: $$ \begin{equation} \boldsymbol{U}_{n+1} = \begin{bmatrix} (u_0)_{n+1} \\ (u_1)_{n+1} \end{bmatrix} = \boldsymbol{U}_n + \frac{h}{2} \cdot \left[\boldsymbol{U}_{n}' + \boldsymbol{U'}_{n+1}^{p}\right] \end{equation} $$

MMS prosedyre:

Velger $ u_m $ = $ sin(t^2) $.
Regner ut $ g = u_{m}''+u_{m}^2 $.
Velger starttid (og sluttid), og regner ut initial-verdier fra $ u_m $. $ t_0 = 0 $, $ t_{end} = 2\pi $.
$ \boldsymbol{U}_0 $ = $ [u_m(t_0), u_{m}'(t_0) ] $
Løser problemet numerisk.

Kan løse for flere $ h $ verdier for å regne ut den observerte ordenen: $$ \begin{align} \label{eq:discretization__error_order} p = \frac{log\left(\frac{{\epsilon}_1}{{\epsilon}_2} \right)}{log\left(r \right)}, \quad \epsilon = max \left( abs(Ua_i- Un_i)\right) \nonumber \end{align} $$

# lectures/lec-16.02/example1.py

import numpy as np
import sympy as sp
import matplotlib.pylab as plt
from math import pi, sqrt

def dUfunc(U, tau):
    
    u0 = U[0]
    u1 = U[1]
    
    du0 = u1
    du1 = g_py(tau) - u0**2

    return np.array([du0, du1])

# main program starts here:
t = sp.symbols('t')
um = sp.sin(t**2)

du = sp.diff(um, t)
ddu = sp.diff(du, t)

g = ddu + um**2

um_py = sp.lambdify(t, um, np)
du_py = sp.lambdify(t, du, np)
g_py = sp.lambdify(t, g, np)

t0, tend = 0, 2*pi
N = 200*2
time = np.linspace(t0, tend, N)
h = time[1]-time[0]

U_0 = [um_py(t0), du_py(0)]

U = np.zeros((N, 2))

U[0][0]= U_0[0]
U[0][1]= U_0[1]

for n, tau in enumerate(time[0: -1]):
    Un = U[n, :]
    Up = Un + h*dUfunc(Un, tau)
    Un1 = Un + 0.5*h*(dUfunc(Un, tau) + dUfunc(Up, tau + h))
    U[n+1, :] = Un1
    
u_num = U[:, 0]
u_man = um_py(time)

rms = np.max(np.abs(u_num-u_man))
print "rms error: ", rms

plt.figure()
plt.plot(time, u_man, 'k')
plt.plot(time, u_num, 'r--')
plt.legend(['manufactured', 'numeric'], frameon=False)
plt.show()

Generalisert Heuns metode

$$ \begin{align*} \boldsymbol{U}_{n+1}&=\boldsymbol{U}_n+h\cdot \left[\frac{1}{2}k1 + \frac{1}{2}k2 \right]& \end{align*} $$ $$ \begin{align*} k1 = \boldsymbol{U}'\left(\boldsymbol{U}_n, t_n\right), \quad k2 = \boldsymbol{U}'\left(\boldsymbol{U}_n + h \cdot k1, t_n + h\right) \end{align*} $$

# lectures/lec-16.02/example2.py

import numpy as np
import sympy as sp
import matplotlib.pylab as plt
from math import pi, sqrt

def dUfunc(U, tau):
    
    u0 = U[0]
    u1 = U[1]
    
    du0 = u1
    du1 = g_py(tau) - u0**2

    return np.array([du0, du1])

def heun_gen(dUfunc, U_0, time):
    
    U = np.zeros((len(time), len(U_0)))
    U[0,:] = U_0
    for n, tau in enumerate(time[0:-1]):
        h = time[n + 1] - time[n]
        Un = U[n,:]
        
        k1 = dUfunc(Un, tau)
        k2 = dUfunc(Un + h*k1, tau + h)
        
        U[n + 1,:] = Un + h*(0.5*k1 + 0.5*k2)
    
    return U

# main program starts here:
t = sp.symbols('t')
um = sp.sin(t**2)

du = sp.diff(um, t)
ddu = sp.diff(du, t)

g = ddu + um**2

um_py = sp.lambdify(t, um, np)
du_py = sp.lambdify(t, du, np)
g_py = sp.lambdify(t, g, np)

t0, tend = 0, 2*pi
N = 20
time = np.linspace(t0, tend, N)

U_0 = [um_py(t0), du_py(0)]

Ntds = 7
orderList = []

for i in range(Ntds):
    
    U = heun_gen(dUfunc, U_0, time)
    u_num = U[:, 0]
    u_man = um_py(time)

    error = np.max(np.abs(u_num-u_man))
    if i>0:
        O = np.log(prev_error/error)/np.log(2)
        orderList.append(O)
    prev_error = error
    N *= 2
    time = np.linspace(t0, tend, N)
    

print orderList

plt.figure()
plt.plot(orderList)
plt.show()