I have the following inhomogeneous parabolic initial/boundary value problem: $$u_{t}(t,x) = (1-x^{2})u_{xx}(t,x)+u(t,x),$$ for $t \in [0,1]$ and $x \in [-1,1]$ $$u(0,x) = \sin(\pi x),$$ for $x \in [-1,1]$ initial condition $$u(t,-1)=u(t,1)=0,$$ $t \in [0,1]$ Dirichlet boundary conditions.
I want to construct a Backward Euler method with $N_{x} =39$ and $Nt=400$
and a Crank Nicolson method with $N_{x} =39$ and $Nt=20$ but I don't know how to put $$a = (1-x^2)$$ inside my method in the script below.Any help?
Backward Euler
def g(x):
return(np.sin(np.pi*x))
Nx = 39
Nt = 400
L = 1
dx = (L - (-L))/(Nx - 1)
t0 = 0
Tf = 1
dt = (Tf - t0)/(Nt - 1)
h = (L - (-L))/(Nx+1)
t = Tf / Nt
m = t/h**2
print("m =", round(m))
x = np.linspace(-L, L, Nx+1)
t = np.linspace(t0, Tf, Nt+1)
a = np.array([1-x**2]).reshape(Nx+1)
u = np.zeros(Nx+1)
u_n = np.zeros(Nx+1)
A = np.zeros((Nx+1, Nx+1))
b = np.zeros(Nx+1)
for i in range(1, Nx):
A[i,i-1] = -m
A[i,i+1] = -m
A[i,i] = 1 + 2*m
A[0,0] = A[Nx,Nx] = 1
A = A*a
#--- initial condition u(x,0) = g(x)
for i in range(0, Nx+1):
u_n[i] = g(x[i])
for n in range(0, Nt):
# Compute b and solve linear system
for i in range(1, Nx):
b[i] = -u_n[i]
b[0] = b[Nx] = 0
u[:] = scipy.linalg.solve(A, b)
# Update u_n before next step
u_n[:] = u
plt.plot(u)
plt.show()
Crank - Nicolson
from scipy.sparse.linalg import spsolve
Nx = 39
Nt = 20
L = 1
dx = (L - (-L))/(Nx - 1)
t0 = 0
Tf = 1
dt = (Tf - t0)/(Nt - 1)
h = (L - (-L))/(Nx+1)
t = Tf / Nt
m = t/h**2
print("m =", round(m))
x = np.linspace(-L, L, Nx+1)
t = np.linspace(t0, Tf, Nt+1)
# Representation of sparse matrix and right-hand side
main = np.zeros(Nx+1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b = np.zeros(Nx+1)
# Precompute sparse matrix
main[:] = 1+m
lower[:] = -1/2*m
upper[:] = -1/2*m
# Insert boundary conditions
main[0] = 0
main[Nx] = 0
A = scipy.sparse.diags(
diagonals=[main, lower, upper],
offsets=[0, -1, 1], shape=(Nx+1, Nx+1),
format='csr')
A = A*a
print(A)
# Set initial condition
for i in range(0,Nx+1):
u_n[i] = g(x[i])
for n in range(0, Nt):
b = u_n
b[0] = b[-1] = 0.0 # boundary conditions
u[:] = scipy.sparse.linalg.spsolve(A, b)
u_n[:] = u
plt.plot(u)
plt.show()
```