计算科学 - 如何在 Python 中的隐式方法的非齐次抛物线 pde 中插入 a(x) 函数？ - 吾爱随笔录

我有以下不均匀抛物线初始/边界值问题：

u_{t} (t, x) = (1 - x^{2}) u_{x x} (t, x) + u (t, x),

$u_{t}(t,x) = (1-x^{2})u_{xx}(t,x)+u(t,x),$ 为了

t \in [0, 1]

$t \in [0,1]$ 和

x \in [- 1, 1]

$x \in [-1,1]$

u (0, x) = \sin (π x),

$u(0,x) = \sin(\pi x),$ 为了

x \in [- 1, 1]

$x \in [-1,1]$ 初始条件

u (t, - 1) = u (t, 1) = 0,

$u(t,-1)=u(t,1)=0,$

t \in [0, 1]

$t \in [0,1]$ 狄利克雷边界条件。

我想构造一个反向欧拉方法 $N_{x} =39$ 和 $Nt=400$

和曲柄尼科尔森方法 $N_{x} =39$ 和 $Nt=20$ 但我不知道怎么写

a = (1 - x^{2})

$a = (1-x^2)$ 在下面脚本中的我的方法中。有帮助吗？

后向欧拉


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

曲柄 - 尼科尔森

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