0
$\begingroup$

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()
```
$\endgroup$
1
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Thank you for your precious and to the point answer as always.But how about the next term u(t,x)? $\endgroup$ – Mr.Podilatis Nov 16 '20 at 23:45
  • $\begingroup$ 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 $\endgroup$ – VoB Nov 16 '20 at 23:54
  • $\begingroup$ Ok.So my code is completely wrong.Thanks again for your assistance.You are awesome mate. $\endgroup$ – Mr.Podilatis Nov 17 '20 at 0:14
  • $\begingroup$ 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. $\endgroup$ – VoB Nov 17 '20 at 0:25
  • $\begingroup$ No worries.If you find time it’s ok.If no again ok.Cheers mate $\endgroup$ – Mr.Podilatis Nov 17 '20 at 0:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.