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

        拉普拉斯方程(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)
        得到如下图所示的位场图:


  • 16
    点赞
  • 82
    收藏
    觉得还不错? 一键收藏
  • 9
    评论
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值