Forelesning 8. april
Apr 7, 2016
Elliptic and parabolic partial differential equations
Problem 1: Stationary temperature field in beam cross-section
In this problem we shall solve an elliptical PDE. We are going to find the stationary solution of the temperature field in a quadratic beam cross-section, see Figure 1. The temperature is known on three of the surface of the beam, and the derivative is known on the last. So this is a combination of Dirichlet- and Neumann boundary conditions.
Figure 1: Beam cross-section.
If we assume isotropic heat conduction (heat conduction coefficient is equal in all directions) we can write the stationary heat conduction equation as $$ \begin{equation} \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=0 \label{eq:Eliptic_pde} \end{equation} $$
Discretize the geometry with equal number of nodes in \( x \)-, and \( y \)-direction (choose a suitable step size). Choose a suitable difference scheme, and find the stationary temperature field in the beam cross-section (solution of \eqref{eq:Eliptic_pde}) numerically.
numerical solution:
Using central differences we get the familiar discretized equation: $$ \begin{equation} T_{i+1,j}+T_{i-1,j}+T_{i,j+1}+T_{i,j-1}-4 T_{i,j}=0 \label{eq:stencil} \end{equation} $$
Figure 2: Illustration of the numerical stencil.
We number the unknowns as shown in Fig. 3. Note that on the line \( y=0 \) we also have unknowns.
Figure 3: Illustration of discretized grid.
Applying Eq. \eqref{eq:stencil} on all interior nodes we get the following matrix equation $$ \begin{equation} \left(\begin{array}{ccccccccccccccc} -4& 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1& -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0& 1 & -4 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1& 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0& 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0& 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0& 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\ 0& 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\ 0& 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ 0& 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\ 0& 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\ 0& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 \end{array}\right) \cdot \left(\begin{array}{c} T_{1}\\ T_{2}\\ T_{3}\\ T_{4}\\ T_{5}\\ T_{6}\\ T_{7}\\ T_{8}\\ T_{9}\\ T_{10}\\ T_{11}\\ T_{12} \end{array}\right) = \left(\begin{array}{c} -10\\ 0\\ -10\\ -10\\ 0\\ -10\\ -10\\ 0\\ -10\\ -40\\ -30\\ -40 \end{array}\right) \label{eq:7205} \end{equation} $$
# lec-8-apr/Problem1.py
import numpy as np
import scipy
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg
import matplotlib.pylab as plt
import time
from math import sinh
#import matplotlib.pyplot as plt
# Change some default values to make plots more readable on the screen
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
def plot_SurfaceNeumann_xy(Temp, Ttop, Tright, Tleft, xmax, ymax, Nx, Ny, nxTicks=4, nyTicks=4):
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# surfaceplot:
x = np.linspace(0, xmax, Nx + 2)
y = np.linspace(0, ymax, Ny + 1)
X, Y = np.meshgrid(x, y)
T = np.zeros_like(X)
T[-1,:] = Ttop
T[:,-1] = Tright
T[:,0] = Tleft
k = 1
for j in range(Ny):
T[j,1:-1] = Temp[Nx*(k-1):Nx*k]
k+=1
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, T, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_zlim(0, Ttop)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('T [$^o$C]')
xticks=np.linspace(0.0,xmax,nxTicks+1)
ax.set_xticks(xticks)
yticks=np.linspace(0.0,ymax,nyTicks+1)
ax.set_yticks(yticks)
def setup_LaplaceNeumann_xy(Ttop, Tright, Tleft, nx, ny):
""" Function that returns A matrix and b vector of the laplace Neumann heat problem A*T=b using central differences
and assuming dx=dy, based on numbering with respect to x-dir, e.g:
30 30 30 -4 1 0 2 0 0 0 0 0 0 0 0 -10
10 T10 T11 T12 10 1 -4 1 0 2 0 0 0 0 0 0 0 0
10 T7 T8 T9 10 0 1 -4 0 0 2 0 0 0 0 0 0 -10
10 T4 T5 T6 10 --> A= 1 0 0 -4 1 0 1 0 0 0 0 0 ,b = -10
10 T1 T2 T3 10 0 1 0 1 -4 1 0 1 0 0 0 0 0
T4 T5 T6 0 0 1 0 1 -4 0 0 1 0 0 0 -10
0 0 0 1 0 0 -4 1 0 1 0 0 -10
0 0 0 0 1 0 1 -4 1 0 1 0 -0
0 0 0 0 0 1 0 1 -4 0 0 1 -10
0 0 0 0 0 0 1 0 0 -4 1 0 -40
0 0 0 0 0 0 0 1 0 1 -4 1 -30
0 0 0 0 0 0 0 0 1 0 1 -4 -40
T = [T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12]^T
Args:
nx(int): number of internal elements in each row in the grid, nx=3 in the example above
ny(int): number of internal elements in each column in the grid, ny=4 in the example above
Returns:
A(matrix): Sparse matrix A, in the equation A*T = b
b(array): RHS, of the equation A*t = b
"""
n = (nx)*(ny) # number of unknowns
b = np.zeros(n) # RHS
d0 = -4*np.ones(n) # main diagonal
d1_lower = np.ones(n - 1) # first lower diagonal
d1_upper = d1_lower.copy()
dnx_lower = np.ones(n - nx) # the nx-th lower diagonal
dnx_upper = dnx_lower.copy()
d1_lower[?::?] = 0 # every nx element on first diagonal is zero; starting from the nx-th element
d1_upper[?::?] = 0
dnx_upper[?:?] = 2 # the first nx elements in the nx-th upper diagonal is two;
# This correspond to all equations on border (x, y=0)
b[?:] += -Ttop
b[?::?] += -Tright
b[?::?] += -Tleft
A = scipy.sparse.diags([d0, d1_upper, d1_lower, dnx_upper, dnx_lower], [0, 1, -1, nx, -nx], format='csc')
return A, b
if __name__ == '__main__':
# Main program
# Set temperature at the top
Ttop=30
Tright = 10
Tleft = 10
xmax=1.0
ymax=1.
Nx = 4
h=xmax/Nx
Ny = int(ymax/h)
A, b = setup_LaplaceNeumann_xy(Ttop, Tright, Tleft, Nx - 1, Ny)
Temp = scipy.sparse.linalg.spsolve(A, b)
plot_SurfaceNeumann_xy(Temp, Ttop, Tright, Tleft, xmax, ymax, Nx - 1, Ny)
plt.show()
Problem 2: Transient temperature field in one-dimensional beam
In this problem we shall solve a parabolic PDE. We are going to find the temperature field along a one-dimensional beam as function of time (transient temperature field). The temperature field is known in the initial state \( t=0 \), and the temperature is known at the boundary.
This problem is goverened by the one-dimensional transient heat conduction equation, $$ \begin{equation} \frac{\partial T}{\partial t}=k \frac{\partial^2 T}{\partial x^2} \label{eq:pde} \end{equation} $$ where \( k \) is the heat conduction coefficient. In this exercise we set \( k=1 \). The beam has length \( L=1 \).
Boundary conditions: \( T(0,t)=100 \) and \( T(1,t)=0 \).
Initial conditions: \( T(x,0)=0 \) Find the temperature field along the beam as function of time. Discretize the length with a suitable number of nodes that you think are sufficient to describe the temperature field while not too computationally intensive.
a) Discretize \eqref{eq:pde} with the FTCS scheme and solve it numerically. Program the scheme both with and without the use of numpy in the integration and compare the CPU time for the two methods. Try different values for \( D \) (the numerical Fourier number), e.g. \( D=0.2 \), \( D=0.5 \) and \( D=1.0 \). Plot the temperature field for \( t=0.001 \), \( t=0.025 \) and \( t=0.4 \).
b) Same as a), but now using the \( \theta \)-scheme instead of the FTCS scheme. Use both \( \theta = \frac{1}{2} \) (Crank) and \( \theta = 1 \) (Laasonen).
numerical solution
a) The following difference equation is obtained using FTCS: $$ \begin{equation} \label{eq:FTCS} u_j^{n+1}= D\, u^n_{j-1}+(1-2D) \, u_j^n + D\, u_{j+1}^n \end{equation} $$
and adding the boundaries: $$ \begin{equation} \label{eq:bc} u_0^{n} = 100, \quad u_{-1}^n = 0 \end{equation} $$
b) The following difference equation is obtained using the general \( \theta \)-scheme: $$ \begin{equation} \label{eq:theta} -D \theta \, u^{n+1}_{j-1} + \left(1 + 2D \, \theta \right) u^{n+1}_{j} -D \theta \, u^{n+1}_{j+1}= D(1-\theta)u^n_{j-1} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{j} + D(1-\theta)u^n_{j+1} \end{equation} $$
And writing the equation especially for \( j = 1 \), we get: $$ \begin{equation} \label{eq:theta_0} \left(1 + 2D \, \theta \right) u^{n+1}_{1} -D \theta \, u^{n+1}_{2}= D(1-\theta)u^n_{0} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{1} + D(1-\theta)u^n_{2} + D \theta \, u^{n+1}_{0} \end{equation} $$
where \( u^{n+1}_{0} \) is known.
# lec-8-apr/Problem2.py
# The equation solved is the parabolic equaiton
# du d du
# -- = k -- --
# dt dx dx
# along with boundary conditions
import numpy as np
import scipy.sparse.linalg
def FTCS_python(uOld, D):
"""Solve using FTCS using a Python expression"""
u = np.zeros_like(uOld)
for j in range(1,nx):
u[j] = D*uOld[?] + (1.0 - 2.0*D)*uOld[?] + D*uOld[?]
u[0] = uleft
u[-1] = uright
return u
def FTCS_numpy(uOld, D):
"""Solve using FTCS using a Numpy expression"""
u = np.zeros_like(uOld)
uOld_minus = uOld[?:?]
uOld_mid = uOld[?:?]
uOld_plus = uOld[?:?]
u[1:-1] = D*uOld_minus + (1.0 - 2.0*D)*uOld_mid + D*uOld_plus
u[0] = uleft
u[-1] = uright
return u
def theta_scheme_python(uOld, D):
"""Solve using the theta scheme"""
u = np.zeros_like(uOld)
subDiag = -D*theta*np.ones(nx - 2)
mainDiag = (1 + 2.0*D*theta)*np.ones(nx - 1)
superDiag = -D*theta*np.ones(nx - 2)
RHS = np.zeros_like(mainDiag)
for j in range(1, nx):
RHS[j - 1] = D*(1-theta)*uOld[?] + (1 - 2*D*(1-theta))*uOld[?] + D*(1-theta)*uOld[?]
RHS[0] += D*theta*uleft
A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc') # sparse matrix instance
u[1:-1] = scipy.sparse.linalg.spsolve(A, RHS)
u[0] = uleft
u[-1] = uright
return u
def theta_scheme_numpy(uOld, D):
"""Solve using the theta scheme"""
u = np.zeros_like(uOld)
subDiag = -D*theta*np.ones(nx - 2)
mainDiag = (1 + 2.0*D*theta)*np.ones(nx - 1)
superDiag = -D*theta*np.ones(nx - 2)
uOld_minus = uOld[?:?]
uOld_mid = uOld[?:?]
uOld_plus = uOld[?:?]
RHS = D*(1-theta)*uOld_minus + (1 - 2*D*(1-theta))*uOld_mid + D*(1-theta)*uOld_plus
RHS[0] += D*theta*uleft
A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc') # sparse matrix instance
u[1:-1] = scipy.sparse.linalg.spsolve(A, RHS)
u[0] = uleft
u[-1] = uright
return u
if __name__ == '__main__':
import matplotlib.pylab as plt
import time as timemodule
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
## Main program starts here
k = 1.0
nx = 10 # number of dx-spaces
L = 1.0 # length of beam
dx = L/nx
x = np.linspace(0, L, nx + 1)
tstart, tend = 0, 0.4 # time period
D = 0.2 # numerical Fourier number
dt = D*dx**2/k # compute timestep based on Fourier number, dx and diffusivity
theta = 0.5 # theta parameter for the theta scheme
uleft, uright, uinitial = 100.0, 0.0, 0.0 # boundary and initial conditions
time = np.arange(tstart, tend + dt, dt)
solvers = [FTCS_nump] #, FTCS_python, theta_scheme_numpy, theta_scheme_python]
uSolutions = np.zeros((len (solvers), len(time), nx + 1))
uSolutions[:, 0, 0] = uleft
uSolutions[:, 0, -1] = uright
for k, solver in enumerate(solvers):
cpu_start = timemodule.time()
for n in range(len(time) - 1):
uOld = uSolutions[k, n, :]
uSolutions[k, n + 1, :] = solver(uOld, D)
cpu_end = timemodule.time()
print "cpu-time for " + solver.func_name + " = ", cpu_end - cpu_start
t_plot_list = [0.001, 0.025, 0.4]
l_style=['c', 'ro', 'g^', 'b--']
legendList = []
plt.figure()
for t in t_plot_list:
n = round(t/dt)
for k in range(len(solvers)):
plt.plot(x, uSolutions[k, n, :], l_style[k])
for solver in solvers:
legendList.append(solver.func_name)
plt.legend(legendList, frameon=False)
plt.show()
Problem 3: Oppstart av rørstrømning
Dersom vi transformerer diffusjonsligningen \( \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2} \) til sylinder-koordinater og forlanger at \( u \) bare skal være funksjon av tiden \( t \) og radien \( r \), får vi: $$ \begin{equation*} \frac{\partial u}{\partial t}= \frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial u }{\partial r} \end{equation*} $$
Figure 4:
Figuren viser hastighetsprofilet ved strømning i et rør av en inkompressibel fluid ved et gitt tidspunkt. Vi tenker oss at profilet har utviklet seg fra null ved å sette på en konstant trykkgradient \( \frac{dp}{dz} < 0 \), slik at hastighetsprofilet for det stasjonære tilfellet er det velkjente parabolske profilet for Poiseuille-strømning.
Bevegelsesligning: $$ \begin{equation} \label{eq:56010} \frac{\partial U}{\partial \tau}= -\frac{1}{\rho}\frac{dp}{dz}+\nu \left( \frac{\partial^2 U}{\partial R^2}+\frac{1}{R}\frac{\partial U}{\partial R} \right) \end{equation} $$
der hastighetsprofilet \( U=U(R,\ \tau).\ 0\leq R \leq R_0, \) og \( \tau \) er den fysikalske tiden.
Kan vise at problemet kan skrives om til: $$ \begin{equation} \label{eq:56015a} \frac{\partial \omega}{\partial t}= \frac{\partial^2\omega}{\partial r^2}+ \frac{1}{r} \frac{\partial \omega}{\partial r} \end{equation} $$
Randbetingelser: $$ \begin{equation} \label{eq:56015b} \omega(\pm1,t)=0,\ \frac{\partial \omega}{\partial r}(0,t)=0 \end{equation} $$
Startbetingelse: $$ \begin{equation} \label{eq:56015c} \omega(r,0)=u_s=1-r^2 \end{equation} $$
Det opprinnelige problemet er da: $$ \begin{equation} \label{eq:56016} u(r,t)=1-r^2-\omega(r,t) \end{equation} $$
\( \theta \)-skjemaet for diffusjonslikningen i sylindriske koordinater kan skrives: $$ \begin{equation} \label{eq:theta_syl} \begin{array}{l} -D\theta\cdot\left(1-\frac{1}{2j}\right)\cdot \omega^{n+1}_{j-1}+(1+2D\theta)\cdot \omega_j^{n+1}-D\theta \cdot\left(1+\frac{1}{2j}\right)\cdot \omega^{n+1}_{j+1}\\ = D(1-\theta)\cdot \left(1-\frac{1}{2j}\right)\cdot \omega^n_{j-1}+[1-2D(1-\theta)]\cdot\omega^n_j +D(1-\theta)\cdot\left(1+\frac{1}{2j}\right)\cdot\omega^n_{j+1} \end{array} \end{equation} $$
men ser at vi får ett problem ved \( r=0 \); Leddet \( \dfrac{1}{r}\dfrac{\partial u}{\partial r} \) må behandles spesielt for \( r=0 \).
L'Hopitals regel: $$ \begin{equation} \label{eq:5606} \lim_{r\to 0}\frac{1}{r}\frac{\partial u}{\partial r} = 1 \frac{\partial^2u}{\partial r^2}\to \frac{\partial u}{\partial t}= 2\frac{\partial^2u}{\partial r^2} \text{ for }r=0 \end{equation} $$
Diskretiserer \eqref{eq:5606} for \( r=0 \) og utnytter symmetribetingelsen \( \frac{\partial u}{\partial r}(0)=0 \): $$ \begin{equation} \label{eq:5607} (1+4D\theta)\omega^{n+1}_0 -4D\theta \omega_1^{n+1}=\omega_0^n + 4D(1-\theta)(\omega_1^n - \omega_0^n) \end{equation} $$
# lec-8-apr/startup.py;Visualization.py @ git@lrhgit/tkt4140/src/src-ch5/Visualization.py;Startupfuncs.py @ git@lrhgit/tkt4140/src/src-ch5/Startupfuncs.py;
import scipy
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg
import numpy as np
def thetaScheme_numpy(theta, D, N, wOld):
""" solve problem by array-slicing and scipy sparse solver"""
superDiag = np.zeros(N - 1)
subDiag = np.zeros(N - 1)
mainDiag = np.zeros(N)
RHS = np.zeros(N)
j_array = np.linspace(0, N, N + 1)
tmp = D*(1. - theta)
superDiag[1:] = -D*theta*(1 + 0.5/j_array[?:?])
mainDiag[1:] = np.ones(N - 1)*(1 + 2*D*theta)
subDiag[:] = -D*theta*(1 - 0.5/j_array[?:?])
a = tmp*(1 - 1./(2*j_array[?:-?]))*wOld[?:?]
b = (1 - 2*tmp)*wOld[?:?]
c = tmp*(1 + 1/(2*j_array[?:?]))*wOld[?:?]
RHS[1:] = a + b + c
superDiag[0] = -4*D*theta
mainDiag[0] = 1 + 4*D*theta
RHS[0] = (1 - 4*tmp)*wOld[0] + 4*tmp*wOld[1]
A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc')
wNew = scipy.sparse.linalg.spsolve(A, RHS)
wNew = np.append(wNew, 0)
return wNew
def thetaScheme_python(theta, D, N, wold):
""" solve problem by looping and using tdma"""
tmp1 = D*(1. - theta)
subDiag = np.zeros(N)
mainDiag = np.zeros(N)
superDiag = np.zeros(N)
RHS = np.zeros(N)
for j in range(1, N):
fac = ?
subDiag[j] = -D*theta*(1 - fac)
mainDiag[j] = (1 + 2*D*theta)
superDiag[j] = -D*theta*(1 + fac)
RHS[j] = tmp1*(1 - fac)*wold[j - 1] + (1 - 2*tmp1)*wold[j]+ tmp1*(1 + fac)*wold[j + 1] + (1 - 2*tmp1)*wold[j]
mainDiag[0] = 1 + 4*D*theta
superDiag[0] = -4*D*theta
RHS[0] = wold[0] + 4*tmp1*(wold[1] - wold[0])
A = scipy.sparse.diags([subDiag[1:], mainDiag, superDiag[0:-1]], [-1, 0, 1], format='csc')
wNew = scipy.sparse.linalg.spsolve(A, RHS)
wNew = np.append(wNew, [0]) # add boundary
return wNew
if __name__ == '__main__':
import matplotlib.pylab as plt
import time as timeModule
from Visualization import createAnimation
from Startupfuncs import analyticSolution
# change some default values to make plots more readable
LNWDT=3; FNT=9
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
#### main program starts here ####
N = 25 # No. of parts
h = 1./N # length of r-step
D = 0.25 # numerical diffusion number
dt = D*h**2
nmax = 1500 # No. of time-steps
T = nmax*dt
time = np.linspace(0, T, nmax + 1)
r = np.linspace(0, 1, N + 1)
eps = 2.2204*10**-16 #tollerance to be used in calculations of analytical solution
theta = 0
scheme_name = 'FTCS'
wStart = 1 - r**2 # initial profile is the pouseille profile
wOld = wStart
#initialize solution matrices
uSolutions = np.zeros((1, len(time), N + 1))
uSolutions2 = np.zeros((1, len(time), N + 1))
uAnalytic = np.zeros((len(time), N + 1))
# solve numerically with vectorized solver
tic = timeModule.time()
for i in range(nmax):
wNew = thetaScheme_numpy(theta, D, N, wOld)
uNew = 1 - r**2 - wNew
uSolutions[0, i + 1, :] = uNew
wOld = wNew
toc = timeModule.time()
print "vectorized scipy sparse solver CPU-time: ", toc - tic
tic = timeModule.time()
for i in range(nmax):
wNew = thetaScheme_python(theta, D, N, wOld)
uNew = 1 - r**2 - wNew
uSolutions2[0, i + 1, :] = uNew
wOld = wNew
toc = timeModule.time()
print "scipy sparse loop solver CPU-time: ", toc - tic
animate = True
if animate:
# solve analytically:
tic = timeModule.time()
for i in range(nmax):
t = time[i + 1]
for k, r_variable in enumerate(r):
uAnalytic[i + 1, k] = 1 - r_variable**2 - analyticSolution(r_variable, t, eps)
toc = timeModule.time()
print "analytic solution CPU-time: ", toc - tic
createAnimation(uSolutions, uAnalytic, [scheme_name], r, time, jump = int(round(0.002/dt)))