用隐式欧拉求解python中的PDE错误输出

Question

我会尽力解释到底是怎么回事和我的问题。

这是有点马虎，所以不支持乳胶，所以可悲的是，我不得不求助于图像。我希望这样可以。

​

​

​

​

我不知道为什么会倒转，很抱歉。无论如何，这是一个线性系统Ax = b，在这里我们知道A和b，所以我们可以找到x，这是我们在下一步的近似。我们继续这样做，直到时间t_final。

这是密码

import numpy as np

tau = 2 * np.pi
tau2 = tau * tau
i = complex(0,1)

def solution_f(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) + np.exp(tau * i * x) * np.exp((tau2 + 4) * i * t))

def solution_g(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) - np.exp(tau * i * x) * np.exp((tau2 + 4) * i * t))

for l in range(2, 12):
    N = 2 ** l #number of grid points
    dx = 1.0 / N #space between grid points
    dx2 = dx * dx
    dt = dx #time step
    t_final = 1
    approximate_f = np.zeros((N, 1), dtype = np.complex)
    approximate_g = np.zeros((N, 1), dtype = np.complex)

    #Insert initial conditions
    for k in range(N):
        approximate_f[k, 0] = np.cos(tau * k * dx)
        approximate_g[k, 0] = -i * np.sin(tau * k * dx)

    #Create coefficient matrix
    A = np.zeros((2 * N, 2 * N), dtype = np.complex)

    #First row is special
    A[0, 0] = 1 -3*i*dt
    A[0, N] = ((2 * dt / dx2) + dt) * i
    A[0, N + 1] = (-dt / dx2) * i
    A[0, -1] = (-dt / dx2) * i

    #Last row is special
    A[N - 1, N - 1] = 1 - (3 * dt) * i
    A[N - 1, N] = (-dt / dx2) * i
    A[N - 1, -2] = (-dt / dx2) * i
    A[N - 1, -1] = ((2 * dt / dx2) + dt) * i

    #middle
    for k in range(1, N - 1):
        A[k, k] = 1 - (3 * dt) * i
        A[k, k + N - 1] = (-dt / dx2) * i
        A[k, k + N] = ((2 * dt / dx2) + dt) * i
        A[k, k + N + 1] = (-dt / dx2) * i

    #Bottom half
    A[N :, :N] = A[:N, N:]
    A[N:, N:] = A[:N, :N]

    Ainv = np.linalg.inv(A)

    #Advance through time
    time = 0
    while time < t_final:
        b = np.concatenate((approximate_f, approximate_g), axis = 0)
        x = np.dot(Ainv, b) #Solve Ax = b
        approximate_f = x[:N]
        approximate_g = x[N:]
        time += dt
    approximate_solution = np.concatenate((approximate_f, approximate_g), axis=0)

    #Calculate the actual solution
    actual_f = np.zeros((N, 1), dtype = np.complex)
    actual_g = np.zeros((N, 1), dtype = np.complex)
    for k in range(N):
        actual_f[k, 0] = solution_f(t_final, k * dx)
        actual_g[k, 0] = solution_g(t_final, k * dx)
    actual_solution = np.concatenate((actual_f, actual_g), axis = 0)

    print(np.sqrt(dx) * np.linalg.norm(actual_solution - approximate_solution))

它不起作用。至少不是在一开始，它不应该开始这么慢。我应该无条件地稳定，并收敛到正确的答案。

这里出什么问题了？

Answer 1

L2-范数可能是测试收敛性的有用指标，但在调试时并不理想，因为它不能解释问题是什么。尽管您的解决方案应该是无条件稳定的，但向后Euler不一定会收敛到正确的答案。正如前欧拉是出了名的不稳定(反耗散)，向后欧拉是著名的耗散。你的解决方案证明了这一点。数值解收敛到零。对于下一阶近似，克兰克-尼科尔森是一个合理的候选人。下面的代码包含了更通用的theta-方法，这样您就可以调优解决方案的隐式性。theta=0.5给CN，theta=1给BE，theta=0给FE。我还调整了几件事：

• 我选择了一个更合适的时间步长为dt = ( dx **2)/2，而不是dt =dx。后者不能通过CN收敛到正确的解决方案。
• 这是一个次要的注意事项，但由于t_final不一定是dt的倍数，所以您没有在同一时间步骤中比较解决方案。
• 关于你的评论，它是缓慢的:随着空间分辨率的提高，你的时间分辨率也需要增加。即使在使用dt=dx的情况下，也必须执行(1024x1024)* 1024矩阵乘法1024次。我没有发现这个在我的机器上花了很长时间。为了加快速度，我删除了一些不必要的连接，但不幸的是，将时间步骤更改为dt = (dx**2)/2会使事情陷入困境。如果你关心速度，你可以试着用Numba编译。

所有这些，我并没有找到巨大的成功与一致性的CN。我必须设置N=2^6才能在t_final=1上得到任何东西。增加t_final会使情况更糟，减少t_final会使情况更好。根据您的需要，您可以考虑实现TR-BDF2 2或其他线性多步方法来改进这一点。

有情节的代码如下：

import numpy as np
import matplotlib.pyplot as plt

tau = 2 * np.pi
tau2 = tau * tau
i = complex(0,1)

def solution_f(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) + np.exp(tau * i * x) * np.exp((tau2 + 4) * i * t))

def solution_g(t, x):
    return 0.5 * (np.exp(-tau * i * x) * np.exp((2 - tau2) * i * t) - np.exp(tau * i * x) * 
np.exp((tau2 + 4) * i * t))

l=6
N = 2 ** l 
dx = 1.0 / N 
dx2 = dx * dx
dt = dx2/2
t_final = 1.
x_arr = np.arange(0,1,dx)

approximate_f = np.cos(tau*x_arr)
approximate_g = -i*np.sin(tau*x_arr)

H = np.zeros([2*N,2*N], dtype=np.complex)
for k in range(N):
    H[k,k] = -3*i*dt
    H[k,k+N] = (2/dx2+1)*i*dt    
    if k==0:
        H[k,N+1] = -i/dx2*dt
        H[k,-1] = -i/dx2*dt     
    elif k==N-1:
        H[N-1,N] = -i/dx2*dt
        H[N-1,-2] = -i/dx2*dt    
    else:
        H[k,k+N-1] = -i/dx2*dt
        H[k,k+N+1] = -i/dx2*dt
### Bottom half
H[N :, :N] = H[:N, N:]
H[N:, N:] = H[:N, :N]

### Theta method. 0.5 -> Crank Nicolson
theta=0.5
A = np.eye(2*N)+H*theta
B = np.eye(2*N)-H*(1-theta)

### Precompute for faster computations
mat = np.linalg.inv(A)@B

t = 0
b = np.concatenate((approximate_f, approximate_g))
while t < t_final:
    t += dt
    b = mat@b

approximate_f = b[:N]
approximate_g = b[N:]
approximate_solution = np.concatenate((approximate_f, approximate_g))

#Calculate the actual solution
actual_f = solution_f(t,np.arange(0,1,dx))
actual_g = solution_g(t,np.arange(0,1,dx))
actual_solution = np.concatenate((actual_f, actual_g))

plt.figure(figsize=(7,5))
plt.plot(x_arr,actual_f.real,c="C0",label=r"$Re(f_\mathrm{true})$")
plt.plot(x_arr,actual_f.imag,c="C1",label=r"$Im(f_\mathrm{true})$")
plt.plot(x_arr,approximate_f.real,c="C0",ls="--",label=r"$Re(f_\mathrm{num})$")
plt.plot(x_arr,approximate_f.imag,c="C1",ls="--",label=r"$Im(f_\mathrm{num})$")
plt.legend(loc=3,fontsize=12)
plt.xlabel("x")

plt.savefig("num_approx.png",dpi=150)

​

​