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

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

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.

• Thank you for your precious and to the point answer as always.But how about the next term u(t,x)? – user37062 Nov 16 '20 at 23:45
• I just edited my answer. Basically, you need to add $u_i$ to each row as you can see. It means that you're adding the vector $(u_i)_i$ to your system @Mr.Podilatis – VoB Nov 16 '20 at 23:54
• Ok.So my code is completely wrong.Thanks again for your assistance.You are awesome mate. – user37062 Nov 17 '20 at 0:14
• You're welcome @Mr.Podilatis. Actually, I had no look at your code. To find out if you did things properly, you could try to write your $i-$th row and see whether or not it's equal to the one I wrote above. Also, you could check the order of convergence of your method: if it has the order you expect, you implemented it correctly. Otherwise, there's something wrong. – VoB Nov 17 '20 at 0:25
• No worries.If you find time it’s ok.If no again ok.Cheers mate – user37062 Nov 17 '20 at 0:41