有限差分法及matlab实现
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
上述六式相加
2
3
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 2 2 2 ∂ U ∂ U ∂ U ⋯ = 6U i , j , k ∂ x2 ∂ y 2 ∂ z 2 2 代入 ∇ U = ρ U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 = 6U i , j , k∇ 2 U ⋯=6U i , j ,k ρ ⋯
即
∂2 U ∂2 U ∂ x2 ∂ y2 ∂2 U ∂2 U ∂ x2 ∂ y2
U i1, j U i−1, j U i , j 1 U i , j −1 =4U i , j
1.2 算法选择:
下面我们对算法再进行一些讨论。
MATLAB 的 M 语言本身带有矩阵的数据类型，且 MATLAB 具有高效的数值计算功能。所以我 们选择通过 MATLAB 的 M 语言去实现这种算法。 (可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考 虑二维情况。)
若考虑离散的点：
∂U 1 ∂ U 1 ∂ U ⋯ ∂ x 2 ∂ x2 6 ∂ x3 ∂ U 1 ∂2 U 1 ∂3 U U i , j 1, k =U i , j ,k ⋯ ∂ y 2 ∂ y2 6 ∂ y3 ∂ U 1 ∂2 U 1 ∂ 3 U U i , j , k1 =U i , j ,k ⋯ ∂ z 2 ∂ z2 6 ∂ z3 ∂ U 1 ∂2 U 1 ∂3 U U i−1, j , k =U i , j , k− − ⋯ ∂ x 2 ∂ x2 6 ∂ x3 ∂ U 1 ∂2 U 1 ∂3 U U i , j 1, k =U i , j , k− − ⋯ ∂ y 2 ∂ y2 6 ∂ y3 ∂ U 1 ∂2 U 1 ∂ 3 U U i , j , k1 =U i , j , k− − ⋯ ∂ z 2 ∂ z2 6 ∂ z3 U i1, j , k =U i , j ,k
迭代解法 :
程序如下：
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 for j=2:lx-1; v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; % 进行迭代计算 t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end k=419 % 输出迭代次数 % 定义矩阵维数 % 建立一个矩阵 % 设置边界条件
由于所有的奇次项被抵消，所以方程：
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 = 6U i , j , k ρ ----- (1) 式
= 6U i , j , k ρ
这就是泊松方程的有限差分形式，以下估计该方程的精度： 由泰勒公式，易知有以下结果：
∂U 1 ∂2 U 2 1 ∂3 U 3 h h h ⋯ ∂x 2 ∂ x2 6 ∂ x3 ∂U 1 ∂ 2 U 2 1 ∂3 U 3 U x , y k, z =U x , y , z k k k ⋯ ∂y 2 ∂ y2 6 ∂ y3 2 3 ∂U 1∂ U 2 1∂ U 3 U x , y , z l =U x , y , z l l l ⋯ ∂z 2 ∂ z2 6 ∂ z3 U x h , y , z =U x , y , z
h
dx
用有限的 h 代替，使得
f ' x ≈
△y △x
差分的种类：
f x h − f x f x − f x −h f x h − f x −h 或者 或者 h 2h h f x h − f x f x − f x −h − h h f x h f x −h −2f x 二阶差分: = h h2 设 U x , y , z 为空间电势的函数 。
有限差分法解静电场的边值问题的算法实现及相关问题讨论：
王宁远
中国科学技术大学 09 级物理 2 班 E-mail wny@mail.ustc.edu.cn
摘要： 本文用 MATLAB 实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程 组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这 样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了 模型,讨论了几个与文[1]所持的不同的观点。
改为
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 的值是否能够减少迭代次数？ 我们做了如下试验:
传统 解法 :
程序如下：
a=zeros(135,135); for i=1:135 a(i,i)=1;end; for i=1:7 a(15*i+1,15*i+2)=-0.25;a(15*i+1,15*i+16)=-0.25;a(15*i+1,15*i-14)=-0.25; end for i=1:7 a(15*i+15,15*i+14)=-0.25;a(15*i+15,15*i+30)=-0.25;a(15*i+15,15*i)=-0.25; end a(1,2)=-0.25;a(1,16)=-0.25; a(121,122)=-0.25;a(121,106)=-0.25; a(135,134)=-0.25;a(135,120)=-0.25; a(15,14)=-0.25;a(15,30)=-0.25;
一阶差分： 泊松方程：
∇2 U = ρ
使用二阶差分代替泊松方程中的拉普拉斯算符，有：
∇ U=
2
f x h, y , z f x −h , y , z −2f x , y , z ∂2 U ∂2 U ∂2 U =∑ 2 2 2 ∂x ∂y ∂z h2
∑
表示分别对三个变元求差分之和，以下相同
正文： 经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而 很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求 出数值解,在足够的精度上进行逼进,这就引出了有限差分法。
1.1 有限差分法:
有限差分法： 微分： f ' x = f x h − f x h 0 = dy
%精度要求，达到精度要求跳出循环
1.3 算法改进
提前使用新值：因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对 v2（i，j）赋值时，它 周围四个点其中有 2 个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使 用新算出的数值进行迭代，其中： 只要把上述代码中
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4
矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不 仅使方程化为简单的有限差分形式，而且这样的数据类型在 计算机中易于储存和运算。那么 h=k=l=1，并且令 f(x,y,z)=u(x,y,z) 就有
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1
此即拉普拉斯方程的有限差分形式。
----- (2)式
ຫໍສະໝຸດ Baidu
这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使 之易于在计算机上实现。
注：有时我们需要解二维静电场，则方程退化为：
U i1, j U i−1, j U i , j 1 U i , j −1 =4U i , j
计算出的矩阵：
绘图：
subplot(1,2,1),mesh(v2) axis([0,17,0,11,0,100]) subplot(1,2,2),contour(v2,32)
有限差分法求出的电势的分布图像：
在这个例子中，我们对网格的划分并不细致（17×11），但却要去计算一个 135 阶方阵的逆，算法的 时间复杂度和空间复杂度都比较大，而且对于边界条件较复杂的情况，这样的矩阵会使得程序的编写十 分不便，由于我们给出的差分公式本身就是近似公式，求出该方程的精确解没有实际意义，我们只需要 一个符合精度要求的近似解即可，我们所以我们应该考虑一个效率更高的，更容易实现的算法，那么可 以考虑使用迭代法。
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; 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); end
的精度为三阶，四阶及更高阶项被略去。 若满足 ∇ U =0 有
2
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1
= 6U i , j , k
由一个简单的例子出发先讨论： 矩形截面由四块导体板围成，其中一块有 100v，其他三块全部接地（电 势 为 0），求解 平 面 各处 电 势 ：
解： 基于有限差分法，我们得到的差分形式事实上是线性方程组。 由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有 15×9=135 个，而对于 每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均 值）， 从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个 135 阶的方阵，只 要求出其逆就可以了。
2
3
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 2 2 2 ∂ U ∂ U ∂ U ⋯ = 6U i , j , k ∂ x2 ∂ y 2 ∂ z 2 2 代入 ∇ U = ρ U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 = 6U i , j , k∇ 2 U ⋯=6U i , j ,k ρ ⋯
即
∂2 U ∂2 U ∂ x2 ∂ y2 ∂2 U ∂2 U ∂ x2 ∂ y2
U i1, j U i−1, j U i , j 1 U i , j −1 =4U i , j
1.2 算法选择:
下面我们对算法再进行一些讨论。
MATLAB 的 M 语言本身带有矩阵的数据类型，且 MATLAB 具有高效的数值计算功能。所以我 们选择通过 MATLAB 的 M 语言去实现这种算法。 (可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考 虑二维情况。)
若考虑离散的点：
∂U 1 ∂ U 1 ∂ U ⋯ ∂ x 2 ∂ x2 6 ∂ x3 ∂ U 1 ∂2 U 1 ∂3 U U i , j 1, k =U i , j ,k ⋯ ∂ y 2 ∂ y2 6 ∂ y3 ∂ U 1 ∂2 U 1 ∂ 3 U U i , j , k1 =U i , j ,k ⋯ ∂ z 2 ∂ z2 6 ∂ z3 ∂ U 1 ∂2 U 1 ∂3 U U i−1, j , k =U i , j , k− − ⋯ ∂ x 2 ∂ x2 6 ∂ x3 ∂ U 1 ∂2 U 1 ∂3 U U i , j 1, k =U i , j , k− − ⋯ ∂ y 2 ∂ y2 6 ∂ y3 ∂ U 1 ∂2 U 1 ∂ 3 U U i , j , k1 =U i , j , k− − ⋯ ∂ z 2 ∂ z2 6 ∂ z3 U i1, j , k =U i , j ,k
迭代解法 :
程序如下：
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 for j=2:lx-1; v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; % 进行迭代计算 t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end k=419 % 输出迭代次数 % 定义矩阵维数 % 建立一个矩阵 % 设置边界条件
由于所有的奇次项被抵消，所以方程：
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1 = 6U i , j , k ρ ----- (1) 式
= 6U i , j , k ρ
这就是泊松方程的有限差分形式，以下估计该方程的精度： 由泰勒公式，易知有以下结果：
∂U 1 ∂2 U 2 1 ∂3 U 3 h h h ⋯ ∂x 2 ∂ x2 6 ∂ x3 ∂U 1 ∂ 2 U 2 1 ∂3 U 3 U x , y k, z =U x , y , z k k k ⋯ ∂y 2 ∂ y2 6 ∂ y3 2 3 ∂U 1∂ U 2 1∂ U 3 U x , y , z l =U x , y , z l l l ⋯ ∂z 2 ∂ z2 6 ∂ z3 U x h , y , z =U x , y , z
h
dx
用有限的 h 代替，使得
f ' x ≈
△y △x
差分的种类：
f x h − f x f x − f x −h f x h − f x −h 或者 或者 h 2h h f x h − f x f x − f x −h − h h f x h f x −h −2f x 二阶差分: = h h2 设 U x , y , z 为空间电势的函数 。
有限差分法解静电场的边值问题的算法实现及相关问题讨论：
王宁远
中国科学技术大学 09 级物理 2 班 E-mail wny@mail.ustc.edu.cn
摘要： 本文用 MATLAB 实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程 组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这 样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了 模型,讨论了几个与文[1]所持的不同的观点。
改为
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 的值是否能够减少迭代次数？ 我们做了如下试验:
传统 解法 :
程序如下：
a=zeros(135,135); for i=1:135 a(i,i)=1;end; for i=1:7 a(15*i+1,15*i+2)=-0.25;a(15*i+1,15*i+16)=-0.25;a(15*i+1,15*i-14)=-0.25; end for i=1:7 a(15*i+15,15*i+14)=-0.25;a(15*i+15,15*i+30)=-0.25;a(15*i+15,15*i)=-0.25; end a(1,2)=-0.25;a(1,16)=-0.25; a(121,122)=-0.25;a(121,106)=-0.25; a(135,134)=-0.25;a(135,120)=-0.25; a(15,14)=-0.25;a(15,30)=-0.25;
一阶差分： 泊松方程：
∇2 U = ρ
使用二阶差分代替泊松方程中的拉普拉斯算符，有：
∇ U=
2
f x h, y , z f x −h , y , z −2f x , y , z ∂2 U ∂2 U ∂2 U =∑ 2 2 2 ∂x ∂y ∂z h2
∑
表示分别对三个变元求差分之和，以下相同
正文： 经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而 很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求 出数值解,在足够的精度上进行逼进,这就引出了有限差分法。
1.1 有限差分法:
有限差分法： 微分： f ' x = f x h − f x h 0 = dy
%精度要求，达到精度要求跳出循环
1.3 算法改进
提前使用新值：因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对 v2（i，j）赋值时，它 周围四个点其中有 2 个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使 用新算出的数值进行迭代，其中： 只要把上述代码中
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4
矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不 仅使方程化为简单的有限差分形式，而且这样的数据类型在 计算机中易于储存和运算。那么 h=k=l=1，并且令 f(x,y,z)=u(x,y,z) 就有
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1
此即拉普拉斯方程的有限差分形式。
----- (2)式
ຫໍສະໝຸດ Baidu
这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使 之易于在计算机上实现。
注：有时我们需要解二维静电场，则方程退化为：
U i1, j U i−1, j U i , j 1 U i , j −1 =4U i , j
计算出的矩阵：
绘图：
subplot(1,2,1),mesh(v2) axis([0,17,0,11,0,100]) subplot(1,2,2),contour(v2,32)
有限差分法求出的电势的分布图像：
在这个例子中，我们对网格的划分并不细致（17×11），但却要去计算一个 135 阶方阵的逆，算法的 时间复杂度和空间复杂度都比较大，而且对于边界条件较复杂的情况，这样的矩阵会使得程序的编写十 分不便，由于我们给出的差分公式本身就是近似公式，求出该方程的精确解没有实际意义，我们只需要 一个符合精度要求的近似解即可，我们所以我们应该考虑一个效率更高的，更容易实现的算法，那么可 以考虑使用迭代法。
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; 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); end
的精度为三阶，四阶及更高阶项被略去。 若满足 ∇ U =0 有
2
U i1, j , k U i−1, j , k U i , j 1, kU i , j −1, kU i , j , k1 U i , j , k−1
= 6U i , j , k
由一个简单的例子出发先讨论： 矩形截面由四块导体板围成，其中一块有 100v，其他三块全部接地（电 势 为 0），求解 平 面 各处 电 势 ：
解： 基于有限差分法，我们得到的差分形式事实上是线性方程组。 由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有 15×9=135 个，而对于 每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均 值）， 从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个 135 阶的方阵，只 要求出其逆就可以了。