为什么不同dx的二维扩散模拟结果不同？

Question

我正在用python做数值2D扩散模拟代码，我有一个大问题

当我增加空间间隔值（dx）时，模拟结果会改变

为什么会这样？我认为这不是稳定问题

你经历过吗？请告诉我有什么问题

我把我做的一个例子代码

• 二维扩散方程

  dC/dt=D*（D^2C/dx^2+D^2C/dy^2）

• 将上述方程离散化。 （n：时间，i:x轴，j:y轴，dx=dy）

  C[n+1，i，j]=C[n，i，j]+D*（dt/dx^2）*（C[n，i+1，j]+C[n，i-1，j]+C[n，i，j+1]+C[n，i，j-1]-4C[n，i，j]）

import numpy as np
import matplotlib.pyplot as plt

T0 = 1e-6 ## 1 us
L0 = 1e-7 ## 100 nm

## Matrix shift
def matrix_shift(CC,ixmax,iymax):
    CC_pn = np.zeros([ixmax,iymax]);CC_mn = np.zeros([ixmax,iymax])
    CC_np = np.zeros([ixmax,iymax]);CC_nm = np.zeros([ixmax,iymax])
    CC_pn[0:ixmax-1,:] = CC[1:,:]; CC_pn[ixmax-1,:] = CC[ixmax-1,:]
    CC_mn[1:,:] = CC[0:ixmax-1,:]; CC_mn[0,:] = CC[0,:]
    CC_np[:,0:iymax-1] = CC[:,1:]; CC_np[:,iymax-1] = CC[:,iymax-1]
    CC_nm[:,1:] = CC[:,0:iymax-1]; CC_nm[:,0] = CC[:,0]
    CC_sum = CC_pn+CC_mn+CC_np+CC_nm
    return CC_sum

dx = 2e-8/L0 ## 20 nm
itmax = int(1e3); ixmax = int(25*(0.2/dx)+1);iymax = int(39*(0.2/dx)+1)      
d7 = 2.2e-10*(T0/(L0**2));dt = (dx**2)/(6.*d7);uu = (dt)/(dx**2)

time_array= np.linspace(0,itmax-1,int(itmax))*dt*1.e3*T0
xarray = np.linspace(0,ixmax-1,ixmax)*dx*1e9*L0

## Initial conditions
CC7 = np.zeros([ixmax,iymax])
ii = int(2./dx); jj = int(4./dx); CC7[ii,jj]=1.
result_ca = np.zeros(itmax)

for it in range(int(itmax)):
    CC7 = CC7 + d7*uu*(matrix_shift(CC7,ixmax,iymax)-4.*CC7)
    result_ca[it] = CC7[ii,jj]

plt.plot(time_array,result_ca); plt.show()

我希望不同dx的结果是一样的

Answer 1

这更多的是一个数学/数值问题，而不是代码问题。改变空间网格改变结果并不奇怪。实际上，在这里你要描述一个连续的数学方程，这会产生一个误差，这个误差随着网格的粗化而增大。所以dx越小，结果就越精确（时间离散化也是如此）

关于数值不稳定性检查，例如Courant–Friedrichs–Lewy condition，以获得更多关于显式格式发生的原因和时间的信息