In [1]:
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
Backward Euler for unconditional stability¶
In [2]:
D = 4.
L = 1
M = 50
T = 0.1
N = 1000 # critical value for Euler is 2000 (sigma=0.5)
h = L/M # spatial spacing
k = T/N # time step
sigma = D*k/h**2
# matrix for backward Euler
A = np.zeros((M-1,M-1))
I = np.arange(M-1)
A[I,I] = 1 + 2*sigma # set the main diagonal
A[I[1:],I[:-1]] = -sigma # set the subdiagonal
A[I[:-1],I[1:]] = -sigma # set the superdiagonal
def f(x):
fvals = np.zeros_like(x)
fvals[ x>0.5 ] = 1
return fvals
def l(t):
return 0*t
def r(t):
return np.sin(100*t)
#return 0*t
x = np.linspace(0,L,M+1)
w = f(x) # approx to u(x,0)
wa = np.zeros((M+1,N+1)) # place to store values entire grid
wa[:,0] = w
for j in range(N):
# update w
newt = k*(j+1)
b = w[1:-1].copy()
b[ 0] += sigma*l(newt)
b[-1] += sigma*r(newt)
w[1:-1] = np.linalg.solve(A,b) # linear solve (could use tridiagonal solve for greater efficiency)
w[0] = l(newt)
w[-1] = r(newt)
wa[:,j+1] = w
if np.abs(w).max()> 10: break
def plotframe(j):
plt.figure()
plt.plot(x,wa[:,j])
plt.xlim(0,L)
plt.ylim(-2,2)
plt.plot(1.01,-2+4*j/N,'mo',clip_on=False)
plt.title('Backward Euler, t ='+str(round(k*j,3)))
jmax = j
interact(plotframe,j=(0,jmax));
Crank-Nicolson¶
Done quickly, not carefully checked.
In [7]:
D = 4.
L = 1
M = 50
T = 0.1
N = 500
h = L/M # spatial spacing
k = T/N # time step
sigma = D*k/h**2
# matrix for backward Euler
A = np.zeros((M-1,M-1))
I = np.arange(M-1)
A[I,I] = 2 + 2*sigma # set the main diagonal
A[I[1:],I[:-1]] = -sigma # set the subdiagonal
A[I[:-1],I[1:]] = -sigma # set the superdiagonal
def f(x):
fvals = np.zeros_like(x)
fvals[ x>0.5 ] = 1
return fvals
def l(t):
return 0*t
def r(t):
return 0*t
x = np.linspace(0,L,M+1)
w = f(x) # approx to u(x,0)
wa = np.zeros((M+1,N+1)) # place to store values entire grid
wa[:,0] = w
for j in range(N):
# update w
oldt = k*j
newt = k*(j+1)
b = (2-2*sigma)*w[1:-1]
b[:-1] += sigma*w[2:-1]
b[1:] += sigma*w[1:-2]
b[ 0] += sigma*(l(oldt)+l(newt))
b[-1] += sigma*(r(oldt)+r(newt))
#display(b.shape)
w[1:-1] = np.linalg.solve(A,b)
w[0] = l(newt)
w[-1] = r(newt)
wa[:,j+1] = w
if np.abs(w).max()> 10: break
def plotframe(j):
plt.figure()
plt.plot(x,wa[:,j])
plt.xlim(0,L)
plt.ylim(-2,2)
plt.plot(1.01,-2+4*j/N,'mo',clip_on=False)
plt.title('Backward Euler, t ='+str(round(k*j,3)))
jmax = j
interact(plotframe,j=(0,jmax));
Using a banded solver in Backward Euler for efficiency¶
In [26]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_banded ###########
D = 4.
L = 1.
M = 50
T = .01
N = 20
h = L/M # spatial grid spacing
k = T/N # time step
sigma = D*k/h**2
def f(x):
fvals = 0*x
fvals[ (x>.5) & (x<0.7) ] = 1
return fvals
def l(t): return 0.
def r(t): return 0.
x = np.linspace(0,L,M+1)
w = f(x) # set initial values
# 3 by M-1 matrix instead of M-1 by M-1
A = np.zeros((3,M-1)) # just make space for the 3 diagonals
A[0,1:] = -sigma # subdiagonal
A[1,: ] = 1.+2.*sigma # diagonal
A[2,:-1] = -sigma # superdiagonal
#display(A)
plt.plot(x,w,'k')
for j in range(N):
# first update the boundary values
newt = k*(j+1)
w[ 0] = l(newt)
w[-1] = r(newt)
rhs = np.array(w[1:-1])
rhs[ 0] += sigma*w[ 0]
rhs[-1] += sigma*w[-1]
w[1:-1] = solve_banded( (1,1), A, rhs )
plt.plot(x,w) #p.set_ydata(w) # update the curve in the plot
plt.title('j='+str(j).zfill(4))
display( w )
In [ ]: