有限差分法解泊松方程

更新：代码中ii和jj的顺序反了

import numpy as np
from matplotlib import pyplot as plt
from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import spsolve
def FDM_poison(div_f, x, y):
    ii, jj = div_f.shape
    ii2 = ii - 2
    jj2 = jj - 2

    hx = np.mean(x[1:]-x[:-1])
    hy = np.mean(y[1:]-y[:-1]) 
    print('hx=', hx, 'hy=', hy)

    A = lil_matrix((ii2 * jj2, ii2 * jj2))

    diagonals = [2 * (1/hx**2 + 1/hy**2) * np.ones(jj2),
            -1/hy**2 * np.ones(jj2-1),
            -1/hy**2 * np.ones(jj2-1)]
    Dmatrix = diags(diagonals, [0, 1, -1])

    for k in range(ii2):
        p = slice(k * jj2, (k + 1) * jj2)
        A[p, p] = Dmatrix
        if k < ii2 - 1:
            next_p = slice((k + 1) * jj2, (k + 2) * jj2)
            A[p, next_p] = diags([-1/hx**2 * np.ones(jj2)], [0])
            A[next_p, p] = diags([-1/hx**2 * np.ones(jj2)], [0])

    A = A.tocsr()

    phi = np.zeros((ii, jj))
    # solver A*phi=b
    b = -div_f[1:-1, 1:-1].reshape(ii2 * jj2)
    tmp_phi = spsolve(A, b)
    phi[1:-1, 1:-1] = tmp_phi.reshape((ii2, jj2))
    return phi

这样运行就不会出现下述情况。

---

我在用有限差分方法求解泊松方程，在实验过程中，我发现改变网格密度和区间大小会对求解结果产生显著影响。

我使用的基本方程是二维泊松方程，边界条件设为零，设置了一组解析解，以下是我进行的一些实验设置：

实验一：网格为40x40，区间为[0, 1]x[0, 1]。

实验二：增加网格密度到80x80，区间不变。

实验三：x方向网格为80，y方向为40，区间不变。

实验四：x方向区间扩大到[0, 2]，网格点为40x40。

实验五：区间为[0, 2]x[0, 1]，但是网格点调整为x和y方向网格距离相同。

实验中，加密网格能减小误差，但单方向加密时，误差反而增加。
增大x方向的区间，即使调整网格保持网格距离一致，误差依然较大。

为什么即使在x和y方向网格距离一致的情况下，增大区间仍会导致误差增加？是我的代码有问题吗？

---

import numpy as np
from matplotlib import pyplot as plt
from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import spsolve
def FDM_poison(div_f, x, y):
    ii, jj = div_f.shape
    ii2 = ii - 2
    jj2 = jj - 2

    hx = np.mean(x[1:]-x[:-1])
    hy = np.mean(y[1:]-y[:-1]) 
    print('hx=', hx, 'hy=', hy)

    A = lil_matrix((ii2 * jj2, ii2 * jj2))

    diagonals = [2 * (1/hx**2 + 1/hy**2) * np.ones(ii2),
            -1/hy**2 * np.ones(ii2-1),
            -1/hy**2 * np.ones(ii2-1)]
    Dmatrix = diags(diagonals, [0, 1, -1])

    for k in range(jj2):
        p = slice(k * ii2, (k + 1) * ii2)
        A[p, p] = Dmatrix
        if k < jj2 - 1:
            next_p = slice((k + 1) * ii2, (k + 2) * ii2)
            A[p, next_p] = diags([-1/hx**2 * np.ones(ii2)], [0])
            A[next_p, p] = diags([-1/hx**2 * np.ones(ii2)], [0])

    A = A.tocsr()

    phi = np.zeros((ii, jj))
    # solver A*phi=b
    b = -div_f[1:-1, 1:-1].reshape(ii2 * jj2)
    tmp_phi = spsolve(A, b)
    phi[1:-1, 1:-1] = tmp_phi.reshape((ii2, jj2))
    return phi

x = np.linspace(0, 1, 41)
y = np.linspace(0, 1, 41)
Y, X =np.meshgrid(y, x)
f = -2*np.pi**2*np.sin(np.pi*X)*np.sin(np.pi*Y)
phi_exact = np.sin(np.pi*X)*np.sin(np.pi*Y)
phi = FDM_poison(f, x, y)

fig, ax = plt.subplots(1,3,figsize=(12, 3))
mesh1 = ax[0].pcolormesh(X, Y, phi, shading='auto', cmap='RdBu_r')
fig.colorbar(mesh1, ax=ax[0])
ax[0].set_title('numeric solution', fontsize=15, fontweight='bold')

mesh2 = ax[1].pcolormesh(X, Y, phi_exact, shading='auto', cmap='RdBu_r')
fig.colorbar(mesh2, ax=ax[1])
ax[1].set_title('exact solution', fontsize=15, fontweight='bold')

mesh3 = ax[2].pcolormesh(X, Y, phi-phi_exact, vmin=-1e-3, vmax=1e-3, shading='auto', cmap='RdBu_r')
fig.colorbar(mesh3, ax=ax[2])
ax[2].set_title('error', fontsize=15, fontweight='bold')
plt.show()

---