简单的一维色散方程数值解

Question

我是编码新手，正在尝试求解简单的一维色散方程。

方程及边界条件：

a

dC/dx = b

d^2C/dx^2

x = hj, C = C0; x = - inf, C =0

解析解为

C = C0 * exp (a/b*(x-hj))

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt

C0 = 3.5         # Concentration at injection point
h_j = 0.7        # Injection point
L = -100         # Location of closed boundary

alpha = 2.597
beta = 1.7

N=10000 # number of segments
x=np.linspace(L,h_j,N+1) # grid points
# determine dx
dx = x[1]-x[0]

# diagonal elements
a1 = -(alpha*dx)/(2*beta)
a2 = 2
a3 = alpha*dx/(2*beta)-1

# construct A
A = np.zeros((N+2, N+2))
A[0,0] = 1
A[N+1,N+1] = 1
for i in np.arange(1,N+1):
    A[i,i-1:i+2] = [a1,a2,a3]
# construct B
B = np.zeros((N+2))
B[0] = 10^(-10)
B[N+1] = C0

xs=np.linspace(0,h_j,100)
Exact = C0*np.exp(alpha/beta*(xs-h_j))

#solve linear system Au=B

u = np.linalg.solve(A,B)

fig = plt.figure()
ax = fig.add_subplot(111)
res = ax.imshow(A, cmap = "jet", interpolation = 'nearest') 
cb = fig.colorbar(res)

# drop the frist point at x_0, as it is outside the domain
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, u[1:N+2], label= "numerical solution")
ax.plot(xs,Exact, label= 'Analytical solution')
ax.legend()
ax.set_xlim(0,h_j)
ax.set_xlabel("z")
ax.set_ylabel("C(z)")
plt.show()

但是求解的数值解与解析解并不相同。我找不到代码或数学中的错误。

Answer 1

您错误地构建了

a

矩阵项。如果您通过添加一阶项和二阶项的贡献来构建矩阵

#first order contributions
a1_first = 0
a2_first = -dx/alpha
a3_first = +dx/alpha

#second order contributions
a1_second = -(dx*dx)/(2*beta)
a2_second = 0
a3_second = -(dx*dx)/(2*beta)

#summed
a1 = a1_first + a1_second
a2 = a2_first + a2_second
a3 = a3_first + a3_second

来源更清晰，情节看起来也更好了