How to insert a(x) function in non homogeneous parabolic pde for implicit method in Python?

Question

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()
```

Answer 1

I'm not going to debug your code, but since $$(1-x)u_{xx}$$ can be discretized, using second order, centered, finite differences with:

$$(1-x_i) \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$$ where $x_i$ is one of the grid points, $h$ your discretization parameter, this means that the $i$-th row of the second derivative matrix is multiplied by the factor $1-x_i$

The $i-$th row of the system is therefore$$(1-x_i) \frac{u_{i+1}-2u_i+u_{i-1}}{h^2} + u_i$$

where $u_i$ is actually $t \mapsto u_i(t)$, and hence can be integrated in time using a suitable time integration method.