有限差分法及matlab实现

= 6Ui , j , kρ
这就是泊松方程的有限差分形式，以下估计该方程的精度： 由泰勒公式，易知有以下结果：
U
xh
,
y
,
z =U
x
,
y
,
z
∂U ∂x
h
1 2
∂2U ∂ x2
h2
1 6
∂3 U ∂ x3
h3⋯
U
x
,
yk,
z =U
x
,
y
,
z
∂U ∂y
k
1 2
∂2U ∂ y2
k2
1 6
∂3 U ∂ y3
虑二维情况。)
由一个简单的例子出发先讨论：
矩形截面由四块导体板围成，其中一块有 100v，其他三块全部接地（电势为 0），求解平面
各处电势：
解： 基于有限差分法，我们得到的差分形式事实上是线性方程组。
由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有 15×9=135 个，而对于
每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均 值），
f
xh− h
f
x
− h
f
x
−
f h
x−h
=
f
xh
f
x−h−2f h2
x
设 U x , y , z为空间电势的函数 。
泊松方程：
∇2U=ρ
使用二阶差分代替泊松方程中的拉普拉斯算符，有：
∑ ∇
2
U
=
∂2 U ∂ x2
∂2 U ∂ y2
∂2 U ∂ z2
=
f xh, y , z f x−h , y , z−2f x , y , z h2
这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使 之易于在计算机上实现。
注：有时我们需要解二维静电场，则方程退化为：
U i1,
jU i−1,
jU i ,
j1U i ,
j−1=4Ui ,
j
∂2U ∂ x2
∂2 U ∂ y2
即
U
i1,
jU i−1,
jU i ,
出数值解,在足够的精度上进行逼进,这就引出了有限差分法。
1.1 有限差分法:
有限差分法：
微分：
f
' x=
f
xh− h
f
x
h 0=
dy dx
用有限的 h 代替，使得
f
'
x≈
△ △
y x
差分的种类：
一阶差分：
f xh− f x h
或者
f x −f x−h h
或者
f xh− f x−h 2h
二阶差分:
j , k1=U i ,
j
,
k−
∂U ∂z
12
∂2 U ∂ z2
−
1 6
∂3 U ∂ z3
⋯
上述六式相加
U i1, j , kU i−1, j , kUi , j1, kU i , j−1, kU i , j , k1U i , j , k−1
=
6Ui ,
j
,
k
∂2 U ∂ x2
∂2 U ∂ y2
提前使用新值：因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对 v2（i，j）赋值时，它 周围四个点其中有 2 个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使
用新算出的数值进行迭代，其中：
只要把上述代码中
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4
样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了 模型,讨论了几个与文[1]所持的不同的观点。
正文：
经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而
很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求
∑ 表示分别对三个变元求差分之和，以下相同
矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不 仅使方程化为简单的有限差分形式，而且这样的数据类型在
计算机中易于储存和运算。那么 h=k=l=1，并且令 f(x,y,z)=u(x,y,z)
就有
U i1, j , kU i−1, j , kU i , j1, kU i , j−1,kU i , j , k1U i , j , k−1
时间复杂度和空间复杂度都比较大，而且对于边界条件较复杂的情况，这样的矩阵会使得程序的编写十
分不便，由于我们给出的差分公式本身就是近似公式，求出该方程的精确解没有实际意义，我们只需要
一个符合精度要求的近似解即可，我们所以我们应该考虑一个效率更高的，更容易实现的算法，那么可
以考虑使用迭代法。
迭代解法:
程序如下：
当 m=1.3 迭代次数 逐行扫描 122 次 隔行扫描：120 对称扫描：104
当 m=1.4 迭代次数 逐行扫描 96 次 隔行扫描：94 对称扫描：78 次
当 m=1.5 迭代次数 逐行扫描 71 次 隔行扫描：69 对称扫描：52 次
当 m=1.6 迭代次数 逐行扫描 43 次 隔行扫描：40 对称扫描：42 次
改为
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4 解释执行后，k=222;
迭代次数接近原来的一半，可见这样的算法优化是有效的。
为了能使迭代次数进一步减少，考虑让每次迭代获得更多的增量，
那么将：v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4; 变形为：v2(i,j)=v1(i,j)+(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4; 再考虑：v2(i,j)=v1(i,j)+m*(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4 适当改变 m 的值是否能够减少迭代次数？ 我们做了如下试验:
⋯
U
i,
j
, k1=U
i
,
j
, k
∂U ∂z
12
∂2 U ∂ z2
1 6
∂3U ∂ z3
⋯
U
i−1,
j , k=U i ,
j
,
k−
∂U ∂x
1 2
∂2 U ∂ x2
− 16

∂3 U ∂ x3
⋯
U
i,
j1,k=U i ,
j
,
k−
∂U ∂y
1 2
∂2 U ∂ y2
− 16
∂3 U ∂ y3
⋯
U i ,
lx=17;ly=11;
%定义矩阵维数
v1=zeros(ly,lx); for j=2:lx-1 v1(ly,j)=100;
%建立一个矩阵
end
v2=v1;maxt=1;t=0; k=0;
%设置边界条件
while(maxt>1e-6)
%精度要求，达到精度要求跳出循环
k=k+1
maxt=0;
for i=2:ly-1
end end b=a^(-1); c=zeros(135,1); for i=121:135 c(i,1)=25;end d=b*c;s=zeros(11,17);for i=2:16 s(11,i)=100; end for i=1:9
for j=1:15; s(i+1,j+1)=d(15*(i-1)+j,1);
j1U i ,
j−1=4Ui ,
j
∂2U ∂ x2
∂2 U ∂ y2
1.2 算法选择:
下面我们对算法再进行一些讨论。
MATLAB 的 M 语言本身带有矩阵的数据类型，且 MATLAB 具有高效的数值计算功能。所以我 们选择通过 MATLAB 的 M 语言去实现这种算法。
(可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考
有限差分法解静电场的边值问题的算法实现及相关问题讨论：
王宁远
中国科学技术大学 09 级物理 2 班 E-mail wny@
摘要：
本文用 MATLAB 实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程
组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这
for i=2:14 a(i,i-1)=-0.25;a(i,i+1)=-0.25; a(i,i+15)=-0.25; end for i=122:134 a(i,i-1)=-0.25;a(i,i+1)=-0.25; a(i,i-15)=-0.25; end for i=1:7
for j=2:14; a(15*i+j,15*i+j-1)=-0.25; a(15*i+j,15*i+j+1)=-0.25; a(15*i+j,15*i+j+15)=-0.25; a(15*i+j,15*i+j-15)=-0.25;
从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个 135 阶的方阵，只