Noregs teknisk-naturvitenskapelege universitet
Institutt for konstruksjonsteknikk
Kontaktperson: Leif Rune Hellevik
Tlf: 73594535/98283895
Eksamen i fag TKT4140 Numeriske berekningar m/datalab
Dato: 19. mai 2015
Oppgåve 1
Figur 1: Skjematisk framstilling av diskretisering av domenet.
I denne oppgåva skal me sjå på følgjande ikkje-lineære ODE: $$ \begin{align} y^{\prime\prime} + \left(y^\prime\right)^2 + y &= \ln x \label{eq:ODE} \end{align} $$ som har analytisk løysning $$ \begin{align*} y &= \ln x \label{} \end{align*} $$ La \( x \) gå frå \( x=1 \) til \( x=2 \). Grensevilkåra til likning \eqref{eq:ODE} er gitt som $$ \begin{align*} y(1) = 0 \qquad \text{og} \qquad y(2) = \ln 2 \label{} \end{align*} $$
Me diskretiserar domenet \( x \) med \( N=50 \) steg med lengd \( \Delta x \). Sidan verdiane av \( y \) er kjend på rendene har me dermed \( N-1 \) ukjende, som vist i figur 1.
\newpage
a) Set opp likning \eqref{eq:ODE} som eit system av 1.ordens ODE-ar.
b) Bruk program-skjelettet i slutten av denne oppgåva til å løysa likningssystemet i oppg a) med skyteteknikk. Fyll inn kommandoar som trengst.
Merk at programmet nyttar funksjonen rk4 frå modulen ODEschemes til å løysa initialverdiproblemet for kvar startverdi me tippar. Denne funksjonen kan du gå ut frå er kjend, det er ikkje naudsynt å greia ut om denne.
c) Diskretiser likning \eqref{eq:ODE} med andreordens sentraldifferansar for \( y^{\prime \prime} \) og \( y^{\prime} \). Kvifor er det resulterande likningssystemet ikkje-lineært?
Sentraldifferansar i likning \eqref{eq:ODE} gir $$ \begin{align*} \frac{y_{i+1}-2y_i+y_{i-1}}{\Delta x^2} + \left(\frac{y_{i+1}-y_{i-1}}{2\Delta x}\right)^2 + y_i &= \ln x_i \label{} \end{align*} $$ som omordna kan skrivast som $$ \begin{align*} y_{i+1}-2y_i+y_{i-1} + \frac{1}{4}\left(y_{i+1}-y_{i-1}\right)^2 + \Delta x^2 y_i &= \Delta x^2 \ln x_i \label{} \end{align*} $$ eller $$ \begin{align} y_{i+1}-(2-\Delta x^2) y_i+y_{i-1} + \frac{1}{4}\left(y_{i+1}-y_{i-1}\right)^2 &= \Delta x^2 \ln x_i \label{eq:difference_eq} \end{align} $$ Systemet i \eqref{eq:difference_eq} er ikkje-lineært fordi leddet \( \frac{1}{4}\left(y_{i+1}-y_{i-1}\right)^2 \) er ikkje-lineært. Det kan dermed ikkje løysast direkte.
d) Lineariser det ikkje-lineære leddet i den diskretiserte likninga ved hjelp av Newton-linearisering. 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 \).\\
Vis at likningssystemet kan skrivast på følgjande tri-diagonale form: $$ \begin{align} \left(1-\alpha_i\right) y_{i-1}^{m+1} - \left(2-\Delta x^2 \right)y_{i}^{m+1} + \left(1+\alpha_i\right) y_{i+1}^{m+1} &= \Delta x^2 \ln x_i + \alpha_i^2 \label{eq:tridiagonal} \end{align} $$ og skriv ut uttrykket for \( \alpha_i \).\\
Skriv ut likning \eqref{eq:tridiagonal} for rendene (\( i=1 \) og \( i=N-1 \)).
Leddet som skal lineariserast er $$ \begin{align} F\left(y_{i-1}, y_{i}, y_{i+1}\right) &= F\left(y_{i-1}, y_{i+1}\right) = \frac{1}{4}\left(y_{i+1}-y_{i-1}\right)^2 \label{eq:F} \end{align} $$
Me har dermed at $$ \begin{align} \frac{\partial F_m}{\partial y_{i-1}^m}\delta y_{i-1} &= -\frac{1}{2}\left(y_{i+1}^m-y_{i-1}^m\right)\delta y_{i-1} = -\alpha_i \delta y_{i-1} \label{eq:partial1}\\ \frac{\partial F_m}{\partial y_{i+1}^m}\delta y_{i+1} &= \frac{1}{2}\left(y_{i+1}^m-y_{i-1}^m\right)\delta y_{i+1} = \alpha_i \delta y_{i+1} \label{eq:partial2} \end{align} $$
Bruk av \eqref{eq:F}, \eqref{eq:partial1} og \eqref{eq:partial2} i \eqref{eq:newton} gir det lineariserte leddet som $$ \begin{align*} F_{m+1} &= \frac{1}{4}\left(y_{i+1}^m-y_{i-1}^m\right)^2 + \alpha_i\left(\delta y_{i+1} - \delta y_{i-1}\right) = \alpha_i^2 + \alpha_i\left(\delta y_{i+1} - \delta y_{i-1}\right) \label{}\\ &= \alpha_i^2 + \alpha_i\left(y_{i+1}^{m+1} -y_{i-1}^{m+1}\right) - \alpha_i\left(y_{i+1}^{m} -y_{i-1}^{m}\right)\\ &= \alpha_i^2 + \alpha_i\left(y_{i+1}^{m+1} -y_{i-1}^{m+1}\right) - 2\alpha_i^2\\ &= \alpha_i\left(y_{i+1}^{m+1} -y_{i-1}^{m+1}\right) - \alpha_i^2 \end{align*} $$ Me har her brukt at \( \alpha_i = \frac{1}{2}\left(y_{i+1}^m-y_{i-1}^m\right) \).
Likning \eqref{eq:difference_eq} kan no skrivast som $$ \begin{align} \left(1-\alpha_i\right) y_{i-1}^{m+1} - \left(2-\Delta x^2 \right)y_{i}^{m+1} + \left(1+\alpha_i\right) y_{i+1}^{m+1} &= \Delta x^2 \ln x_i + \alpha_i^2 \label{eq:tridiagonal_sol} \end{align} $$
Likning \eqref{eq:tridiagonal_sol} utskriven for venstre rand vert $$ \begin{align*} - \left(2-\Delta x^2 \right)y_{1}^{m+1} + \left(1+\alpha_1\right) y_{2}^{m+1} &= \Delta x^2 \ln x_1 + \alpha_1^2 - \left(1-\alpha_1\right) y_{0} \end{align*} $$ Sidan \( y_0=y(1)=0 \) har me at $$ \begin{align*} - \left(2-\Delta x^2 \right)y_{1}^{m+1} + \left(1+\alpha_1\right) y_{2}^{m+1} &= \Delta x^2 \ln x_1 + \alpha_1^2 \end{align*} $$
Likning \eqref{eq:tridiagonal_sol} utskriven for høgre rand vert $$ \begin{align*} \left(1-\alpha_{N-1}\right) y_{N-2}^{m+1} - \left(2-\Delta x^2 \right)y_{N-1}^{m+1} &= \Delta x^2 \ln x_{N-1} + \alpha_{N-1}^2 - \left(1+\alpha_{N-1}\right) y_{\text{max}} \end{align*} $$ der \( y_{\text{max}}=y(2)=\ln 2 \).
e) Bruk program-skjelettet i slutten av denne oppgåva til å løysa likning \eqref{eq:ODE} som eit randverdiproblem ved å bruka ein direkte løysar på likning \eqref{eq:tridiagonal}. Fyll inn kommandoar som trengst.
\newpage
Program-skjelett
Følgjande skjelett kan brukast til å svara på oppgåve b) og e). Lever sidene saman med resten av svaret.
import matplotlib.pyplot as plt
import scipy.sparse.linalg
import numpy as np
import numpy.linalg as la
from ODEschemes import rk4
import plotstyles # Set default linewidth and fontsize
def frhs(y, x):
"""ODE-system."""
########### Fill in here #############
return
######################################
def dsfunction(phi0, phi1, s0, s1):
if (abs(phi1-phi0) > 0.0):
########### Fill in here #############
return
######################################
else:
print 'Problem. Division by zero.'
return 0.0
xmin, xmax = 1.0, 2.0
N = 50 # no dx-values => N-1 unknowns for bvp
x, dx = np.linspace(xmin, xmax, N+1, retstep=True) # x includes bndrys
# From initial plots of the phi-function we have two initial guesses
s0, s1 = 1.0, 5.0
# Boundary value for x = (xmin, xmax)
ymin, ymax = 0.0, np.log(2.0)
########## Solve with shooting method ################################
# Compute phi0
y0 = [ymin, s0]
yi = rk4(frhs, y0, x)
########### Fill in here #############
phi0 =
######################################
# Iterate to find solution
nmax = 20 # maximum number of iterations
eps = 1.0e-8 # tolerance
for n in range(nmax):
y0 = [ymin, s1]
yi = rk4(frhs, y0, x)
########### Fill in here #############
phi1 =
######################################
ds = dsfunction(phi0, phi1, s0, s1)
s0 = s1
s1 += ds
phi0 = phi1
if ((abs(ds) <= eps) and (abs(yi[-1,0]-ymax)/ymax <= eps)):
print 'Shooting solution converged.'
y_shoot = yi
break
########## Solve with direct non-linear method #######################
# Crete rhs array for uknowns (i.e. excluding boundaries)
d = np.zeros(N-1)
# Create subiteration array for variable
ym = np.zeros(N+1)*ymax # including boundaries
ym[0]= ymin; ym[-1] = ymax # impose boundary values
alpha = np.zeros(N-1)
# Create matrix for sparse solver
diagonals = np.zeros((3,N-1))
diagonals[1,:] = -(2.0-dx**2)
# Iterate to find solution
nmax = 20 # maximum number of iterations
eps = 1.0e-6 # tolerance
for m in range(nmax):
########### Fill in here #############
alpha[:] =
######################################
# Update diagonals
diagonals[0,:-1] = 1-alpha[1:]
diagonals[2,1:] = 1+alpha[:-1]
# Update sparse matrix
As = scipy.sparse.spdiags(diagonals, [-1,0,1], N-1, N-1,format='csc')
# Update rhs
########### Fill in here #############
d[:] =
######################################
# Calculate solution for current iteration
y_direct = scipy.sparse.linalg.spsolve(As,d)
if (la.norm(y_direct-ym[1:-1]) < eps):
print 'Convergence of direct nonlinear solver achieved.'
ym[1:-1] = y_direct
break
ym[1:-1] = y_direct
y_direct = np.append(y_direct,ymax) # append right boundary value
y_direct = np.insert(y_direct,0,ymin) # insert left boundary value
#Plot the solutions
plt.hold('on')
plt.plot(x,y_shoot[:,0],'-')
plt.plot(x,y_direct,'-.')
plt.plot(x,np.log(x),':')
plt.xlabel('x')
plt.ylabel('y')
legends=[]
legends.append('shooting')
legends.append('direct')
legends.append('analytical')
plt.legend(legends,loc='best')
plt.show()
\newpage
\newpage
import matplotlib.pyplot as plt
import scipy.sparse.linalg
import numpy as np
import numpy.linalg as la
from ODEschemes import rk4
import plotstyles # Set default linewidth and fontsize
def frhs(y, x):
"""ODE-system."""
return [y[1], -y[1]**2-y[0]+np.log(x)]
def dsfunction(phi0, phi1, s0, s1):
if (abs(phi1-phi0) > 0.0):
return -phi1*(s1-s0)/float(phi1-phi0)
else:
print 'Problem. Division by zero.'
return 0.0
xmin, xmax = 1.0, 2.0
N = 50 # no dx-values => N-1 unknowns for bvp
x, dx = np.linspace(xmin, xmax, N+1, retstep=True) # x includes bndrys
# From initial plots of the phi-function we have two initial guesses
s0, s1 = 1.0, 5.0
# Boundary value for x = (xmin, xmax)
ymin, ymax = 0.0, np.log(2.0)
########## Solve with shooting method ################################
# Compute phi0
y0 = [ymin, s0]
yi = rk4(frhs, y0, x)
phi0 = yi[-1,0] - ymax
# Iterate to find solution
nmax = 20 # maximum number of iterations
eps = 1.0e-8 # tolerance
for n in range(nmax):
y0 = [ymin, s1]
yi = rk4(frhs, y0, x)
phi1 = yi[-1,0] - ymax
ds = dsfunction(phi0, phi1, s0, s1)
s0 = s1
s1 += ds
phi0 = phi1
if ((abs(ds) <= eps) and (abs(yi[-1,0]-ymax)/ymax <= eps)):
print 'Shooting solution converged.'
y_shoot = yi
break
########## Solve with direct non-linear method ########################
# Crete rhs array for uknowns (i.e. excluding boundaries)
d = np.zeros(N-1)
# Create subiteration array for variable
ym = np.zeros(N+1)*ymax # including boundaries
ym[0]= ymin; ym[-1] = ymax # impose boundary values
alpha = np.zeros(N-1)
# Create matrix for sparse solver
diagonals = np.zeros((3,N-1))
diagonals[1,:] = -(2.0-dx**2)
# Iterate to find solution
nmax = 20 # maximum number of iterations
eps = 1.0e-6 # tolerance
for m in range(nmax):
alpha[:] = 0.5*(ym[2:]-ym[0:-2])
# Update diagonals
diagonals[0,:-1] = 1-alpha[1:]
diagonals[2,1:] = 1+alpha[:-1]
# Update sparse matrix
As = scipy.sparse.spdiags(diagonals, [-1,0,1], N-1, N-1,format='csc')
# Update rhs
d[:] = np.log(x[1:-1])*dx**2
d[:] += alpha[:]**2
d[0] -= (1-alpha[0])*ymin
d[-1] -= (1+alpha[-1])*ymax
# Calculate solution for current iteration
y_direct = scipy.sparse.linalg.spsolve(As,d)
if (la.norm(y_direct-ym[1:-1]) < eps):
print 'Convergence of direct nonlinear solver achieved.'
ym[1:-1] = y_direct
break
ym[1:-1] = y_direct
y_direct = np.append(y_direct,ymax) # append right boundary value
y_direct = np.insert(y_direct,0,ymin) # insert left boundary value
#Plot the solutions
plt.hold('on')
plt.plot(x,y_shoot[:,0],'-')
plt.plot(x,y_direct,'-.')
plt.plot(x,np.log(x),':')
plt.xlabel('x')
plt.ylabel('y')
legends=[]
legends.append('shooting')
legends.append('direct')
legends.append('analytical')
plt.legend(legends,loc='best')
plt.show()
Oppgåve 2
Adveksjonslikninga er gitt ved: $$ \begin{equation} \frac{\partial u}{\partial t} + a_0 \, \frac{\partial u}{\partial x} = 0 \label{eq:advection} \end{equation} $$
a) Bruk sentral-differansar og Taylor-utvikling til å vise at eit Lax-Wendroff-skjema for \eqref{eq:advection} kan skrivast på forma: $$ \begin{equation} u_j^{n+1}=u_j^n-\frac{C}{2}\left( u_{j+1}^n-u_{j-1}^n \right) + \frac{C^2}{2}\left( u_{j+1}^n - 2u_j^n+u_{j-1}^n \right) \label{eq:lw} \end{equation} $$ og finn uttrykket for \( C \).
b) Kva for orden har Lax-Wendroff-skjemaet i likning \eqref{eq:lw}?
Skjemaet er \( O(\Delta t^2) \) og \( O(\Delta x^2) \). To ledd i Taylor-utviklinga gir andre orden i tid og sentral-differansar gir andre orden i rom.
c) Bruk kriteriet om positive koeffisientar (PK-kriteriet) til å vurdera stabiliteten til skjemaet i likning \eqref{eq:lw}.
d) Ved å nytte von Neumann-analyse på likning \eqref{eq:lw} får ein $$ \begin{equation} | G |^2 = 1 - 4C^2 \, (1-C^2 ) \sin^4 \left(\frac{\delta}{2}\right) \label{eq:vonNeumann} \end{equation} $$ Bruk \eqref{eq:vonNeumann} til å vurdera stabiliteten til Lax-Wendroff-skjemaet i \eqref{eq:lw}.
e) Den modifiserte likninga til \eqref{eq:lw} kan skrivast som $$ \begin{equation} \frac{\partial u}{\partial t} + a_0 \, \frac{\partial u}{\partial x} = -a_0 \frac{\Delta x^2}{6} \left (1 - C^2 \right ) \frac{\partial^3 u}{\partial x^3} + \ldots \label{eq:mod} \end{equation} $$ Korleis finn ein den modifiserte likninga i \eqref{eq:mod}?
f) Kva fortel den modifiserte likninga om kva typar feil ein kan vente med skjemaet i likning \eqref{eq:lw}?
g) Den modifiserte likninga for oppstraums-skjemaet for adveksjonslikninga \eqref{eq:advection} er gitt ved: $$ \begin{equation} \frac{\partial u}{\partial t} + a_0 \, \frac{\partial u}{\partial x} = \frac{a_0 \Delta x}{2} \left (1 - C \right ) \, \frac{\partial^2 u}{\partial x^2}+ \ldots \label{eq:mod_upstream} \end{equation} $$ Gi ei fysisk tolking av \( \frac{a_0 \Delta x}{2} \left (1 - C \right ) \).
h) I figur 2 ser du numeriske løysingar av adveksjonslikninga generert med eit Lax-Wendroff-skjema og eit oppstraums-skjema ved eit gitt tidspunkt. Den analytiske trappefunksjonsløysinga er også illustrert i gult. Marker på figurn kva for løysing som tilhøyrer kva for skjema. Lever figurn inn saman med resten av svaret ditt.
Figur 2: Løysingar frå to ulike numeriske skjema for adveksjonslikninga. Analytisk løysing i gult.
