二维拉普拉斯方程的数值解法

        拉普拉斯方程（Laplace's equation）又称调和方程、位势方程，是一种偏微分方程，因由法国数学家拉普拉斯首先提出而得名。

        两个自变量的拉普拉斯方程具有以下形式：

        对于这样一个偏微分方程，我们想要通过计算机求它的数值解：对于一个函数f(x)在x处的一阶导数我们可以用它的中心差分来表示：f'(x)=(f(x+dx/2)-f(x-dx/2))/dx.令dx=1可得f(x)步长为1的一阶导数f'(x)=f(x+0.5)-f(x-0.5).同理f(x)的二阶导数f"(x)可以看作f(x+0.5)-f(x-0.5)的一阶导数，再使用一次中心差分可得f"(x)=f(x+1)-2f(x)+f(x-1)。

        对于二元函数F(x,y),F(x,y)对自变量x的二阶偏导为F(x+1,y)-2F(x,y)+F(x-1,y)，对自变量y的二阶偏导为F(x,y+1)-2F(x,y)+F(x,y-1)。对拉普拉斯方程可以得到F(x+1,y)-2F(x,y)+F(x-1,y)+F(x,y+1)-2F(x,y)+F(x,y-1)=0，即F(x,y)=(F(x+1,y)+F(x-1,y)+F(x,y+1)+F(x,y-1))/4。在已知边界条件和初始条件的情况下可以迭代求和得到满足精度要求的数值解F(x,y)。

        对于以下案例：

        将求解区域划分成100X100的网格，定义初始条件和边界条件之后，使用F(x,y)=(F(x+1,y)+F(x-1,y)+F(x,y+1)+F(x,y-1))/4公式进行迭代计算。以下为python代码我让求解过程迭代10000次：

import numpy as np
import matplotlib.pyplot as plt
m=101
n=101
U=np.zeros([m,n])
F=np.zeros([m,n])
#定义边界条件
F[m-1,:]=100
e=0.001
count=0

for i in range(10000):
    count+=1
    F[:,0]=0
    F[:,n-1]=0
    F[0,:]=0
    F[m-1,:]=100
    U[1:-1,1:-1]=(F[1:-1,:-2]+F[1:-1,2:]+F[:-2,1:-1]+F[2:,1:-1])/4
    e=abs(U-F).max() #求助老师我这个值为何总是0
    F=U
       
x=np.arange(0,m)
y=np.arange(0,n)
X,Y=np.meshgrid(x,y)
C=plt.contour(X,Y,U)
plt.clabel(C, inline = True, fontsize = 10)

        得到如下位场等值线图：

        我将相邻两次迭代的结果定义为误差。我将矩阵各元素的最大误差控制在0.001内，MATLAB代码如下：

clc
clear
m=101;
n=101;
U=zeros(m,n);%电势
F=U;
count=0 %计数器
e=0.001; %要求误差精度
while e>=0.001
    count=count+1;
    %初始条件
    F(1,:)=0;
    F(:,1)=0;
    F(:,n)=0;
    F(m,:)=100;
    U(2:m-1,2:n-1)=(F(2:m-1,1:n-2)+F(2:m-1,3:n)+F(1:m-2,2:n-1)+F(3:n,2:n-1))/4;
    e=abs(max(max(U-F)));
    F=U;
end

x=0:m-1;
y=0:n-1;
[X,Y]=meshgrid(x,y);
figure
[c,h]=contour(X,Y,U,[0:10:100]);
clabel(c,h)

        得到如下图所示的位场图：