超松弛迭代法求解二维电磁场有限差分方程（附Matlab代码）

FTT不会FFT

已于 2022-04-19 21:48:09 修改

阅读量5.1k

点赞数 16

文章标签： matlab

于 2022-04-19 20:24:23 首次发布

本文链接：https://blog.csdn.net/m0_53560590/article/details/124274858

版权

二维电磁场泊松方程差分格式

由泰勒公式：

$f(x+h)=f(x)+h\frac{df(x)}{dx}+\frac{h^{2}}{2!}\frac{d^{2}f(x)}{dx^{2}}+\cdot \cdot \cdot$

以及：

$f(x-h)=f(x)-h\frac{df(x)}{dx}+\frac{h^{2}}{2!}\frac{d^{2}f(x)}{dx^{2}}+\cdot \cdot \cdot$

两式做和，截断于 $h^{2}\frac{d^{2}f(x)}{dx^{2}}$ 项，得到二阶差商：

$\frac{d ^{2}f(x)}{d x^{2}}\approx \frac{ f(x+h)+ f(x-h)-2f(x)}{h^2}$

二维场域内泊松方程为：

$\Delta \varphi =\frac{\partial ^{2}\varphi }{\partial x^2}+\frac{\partial ^{2}\varphi }{\partial y^2}=f$

于是在等距五点格式下，点 $\varphi _{i,j}$ 的差分方程为

$\frac{1}{h^2}(\varphi _{i+1,j}-2\varphi _{i,j}+\varphi _{i-1,j})+ \frac{1}{h^2}(\varphi _{i,j+1}-2\varphi _{i,j}+\varphi _{i,j-1})=f_{i,j}$

对于拉普拉斯方程，则有：

$\varphi _{i+1,j}+\varphi _{i-1,j}+ \varphi _{i,j+1}+\varphi _{i,j-1}-4\varphi _{i,j}=0$

差分方程组的解法

1直接迭代

由于 $\varphi _{i,j}=\frac{1}{4}(\varphi _{i+1,j}+\varphi _{i-1,j}+ \varphi _{i,j+1}+\varphi _{i,j-1}-h^2f_{i,j})$

于是任意给出各节点处函数值 $\varphi ^0_{i,j}$ 然后带入上式右端，得到第一次函数近似值 $\varphi ^1_{i,j}$ .

依次循环，以第n次迭代值求得第n+1次近似值，即：

$\varphi ^{(n+1)}_{i,j}=\frac{1}{4}(\varphi^{(n)} _{i+1,j}+\varphi^{(n)} _{i-1,j}+ \varphi^{(n)} _{i,j+1}+\varphi^{(n)} _{i,j-1}-h^2f_{i,j})$

由于此法收敛慢，又消耗两套储存单元，占用内存大，实际中并不采用。

2高斯—赛德尔迭代法

这一方法的基本思想是：在第n+1次迭代时，如果某些相关节点上的第n+1次迭代值已经得到，便可直接将这些新值代入进行运算。因此，高斯—赛德尔迭代法的公式与计算顺序是有关的。在从左到右，从上到下的逐点计算顺序中，高斯—赛德尔公式如下：

$\varphi ^{(n+1)}_{i,j}=\frac{1}{4}(\varphi^{(n)} _{i+1,j}+\varphi^{(n+1)} _{i-1,j}+ \varphi^{(n)} _{i,j+1}+\varphi^{(n+1)} _{i,j-1}-h^2f_{i,j})$

可以证明，这一方法迭代次数近似与 $h^2$ 成反比。但网格节点数目很大时，这一方法收敛速度仍很慢。

3超松弛迭代法

为了加快收敛，引入松弛因子 $\omega$ ，把高斯—赛德尔迭代法的值作为一个中间结果

$\overline\varphi _{i,j}=\frac{1}{4}(\varphi^{(n)} _{i+1,j}+\varphi^{(n+1)} _{i-1,j}+ \varphi^{(n)} _{i,j+1}+\varphi^{(n+1)} _{i,j-1}-h^2f_{i,j})$

取n+1次迭代为 $\overline\varphi _{i,j}$ 和上次迭代值 $\varphi^{(n)} _{i,j}$ 的加权平均，即

$\varphi ^{(n+1)}_{i,j}=\varphi ^{(n)}_{i,j}+\omega (\overline\varphi _{i,j}-\varphi ^{(n)}_{i,j})$

整理得：

$\varphi ^{(n+1)}_{i,j}=\varphi ^{(n)}_{i,j}+\frac{\omega}{4} (\varphi ^{(n)}_{i+1,j}+\varphi ^{(n)}_{i,j+1}+\varphi ^{(n+1)}_{i-1,j}+\varphi ^{(n+1)}_{i,j-1}-h^2f_{i,j}-4\varphi ^{(n)}_{i,j})$

$\omega$ 取值范围一般为[1,2)，其具体值依靠经验。当 $\omega$ =1时，还原为高斯—赛德尔迭代法。

应用示例：

对于这一题，编写matlab函数实现功能：输入长、宽节点数以及松弛因子，输出最终迭代结果与迭代次数。这是一个简单矩形无源空间的迪利克雷边界条件问题。不难根据上述原理写出matlab代码如下：

function [Phi2,n] = DSOR(a,b,w)

%a列数
%b行数
%w松弛因子

Phi = [100*ones(1,a);zeros(b-1,a)]; %生产点阵并配置边界条件
Phi(2:round(b/2-0.5),2:a-1) = 100;  %稍微填充一下内部，减少迭代次数

Phi2 = Phi;
Phi1=Phi2;  %生成两个矩阵迭代

n=0;    %迭代次数
p=1;    %判断条件，前后两次迭代各元素之差的绝对值的最大值，初始化为1

while p>10E-10  %判断条件，p大于10E-10时继续迭代
   
    for i = 2:b-1
        for j = 2:a-1
            Phi2(i,j) = Phi1(i,j)+w*0.25*(Phi1(i+1,j)+Phi1(i,j+1)+Phi2(i-1,j)+Phi2(i,j-1)-4*Phi1(i,j));
        end
    end     %完成一次迭代
    n=n+1;
    p=max(max(abs(Phi2-Phi1))); %判断条件，前后两次迭代各元素之差的绝对值的最大值
    Phi1=Phi2;
end

如下图参数调用此函数