拉普拉斯方程的数值解
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1
金文璨 2007201230
利用边值，对（*）式进行迭代，迭代足够的次数后，结果收敛，得到稳定的数 值解。
三． Fortran 程序
PROGRAM LAPLACE IMPLICIT NONE REAL u(21,21) REAL x,y INTEGER i,j,k REAL h=0.1 OPEN(1,file='LAPLACE.txt') DO i=1,21 DOj=1,21 u(i,j)=0 ENDDO ENDDO DO i=1,21 u(i,1)=0 u(i,21)=0 u(1,i)=0 u(21,i)=0 ENDDO DO i=8,14 DO j=8,14 u(i,j)=1.0 ENDDO ENDDO
2
金文璨 2007201230
DO i=8,20 DO j=15,20 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO DO i=15,20 DO j=8,14 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO ENDDO DO i=1,21 DO j=1,21 x=-1+(i-1)*h y=1-(j-1)*h WRITE(1,*)x,y,u(i,j) ENDDO ENDDO CLOSE(1) END PROGRAM
四． 计算结果可视化
(1) 电势的分布图：3D 图中，X,Y 为横纵坐标，Φ(x,y)为电势。
3
金文璨 2007201230
(2) 二维等势线图：
�� (3) 电场满足 E = −∇φ ， 即电场沿电势降落的梯度方向， 故电场线垂直于等势线， 并由高电势指向低电势，如下图：
4
金文璨 2007201230
5
! 中间区域的初值，不妨设为 0
! 外壁的边值
! 金属棒的边值
DO k=1,10000 ! 迭代次数 DO i=2,7 DO j=2,20 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO DO i=8,20 DO j=2,7 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO
用差分法表示：
ui, j =
u i +1, j + u i −1, j + u i , j −1 + u i , j +1
4源自文库
（*）
u 0 , j = 0, u N , j = 0, u i ,0 = 0, u i , N = 0
u i , j = 1; x i ∈ [ − 0 .3, 0 .3 ], y i ∈ [ − 0 .3, 0 .3 ]
金文璨 2007201230
作业 6.1
一． 题目
有一个中空的无限长棱柱，其横截面如图所示，棱柱壁是金属的，棱柱中是金属 棒，在金属棒和外壁上分别加电压，计算金属棒和棱柱壁之间的电场，画出等势 线，电势和电场的分布。
二． 算法
金属棒为等势体，金属棒和外壁之间的区域内，电势满足 Laplace 方程。解电势
φ ( x, y ) 的 Laplace 边值问题：
∂ 2φ ∂ 2φ ∆φ = 2 + 2 = 0 ∂x ∂y
φ(x = ±1, y) = 0 ； φ ( x , y = ± 1) = 0 φ ( x , y ) = 1; x ∈ [ − 0 .3 .0 .3], y ∈ [ − 0 .3 .0 .3]
五． 讨论
(1) 迭代之前，金属棒和外壁之间区域的初值可以任意给定，但方便起见，不妨 设为零； (2) 该计算中，因为划分网格的步长比较短（h=0.1） ，格点比较少，所以迭代次 数不用太多（10000 次） ，即可得到稳定的结果。如果格点比较多的话，可以 对程序做一改进， 设定一个阈值， 当前后两次迭代的结果的差值小于阈值时 ， 迭代停止； (3) 由对称性分析可得，电势分布关于 x 轴、y 轴轴对称，同时也关于坐标原点 中心对称，所以如果格点比较多，计算量大的话，可以只计算第一象限的电 势值，然后通过坐标变换得到其他三个象限的电势值； (4) 从理论上分析，从金属棒到外壁，电势均匀降落，等势线应该是同心的圆角 正方形环，电场线应该是从金属棒向外壁辐射，并垂直于等势线； (5) Fortran 程序计算出电势的分布，在结果可视化的过程中，需要 MATLAB 对 电势进行处理，用 contour(x,y,u)语句绘制等势线，用 gradient(u)绘制电场矢 量的分布。
金文璨 2007201230
利用边值，对（*）式进行迭代，迭代足够的次数后，结果收敛，得到稳定的数 值解。
三． Fortran 程序
PROGRAM LAPLACE IMPLICIT NONE REAL u(21,21) REAL x,y INTEGER i,j,k REAL h=0.1 OPEN(1,file='LAPLACE.txt') DO i=1,21 DOj=1,21 u(i,j)=0 ENDDO ENDDO DO i=1,21 u(i,1)=0 u(i,21)=0 u(1,i)=0 u(21,i)=0 ENDDO DO i=8,14 DO j=8,14 u(i,j)=1.0 ENDDO ENDDO
2
金文璨 2007201230
DO i=8,20 DO j=15,20 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO DO i=15,20 DO j=8,14 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO ENDDO DO i=1,21 DO j=1,21 x=-1+(i-1)*h y=1-(j-1)*h WRITE(1,*)x,y,u(i,j) ENDDO ENDDO CLOSE(1) END PROGRAM
四． 计算结果可视化
(1) 电势的分布图：3D 图中，X,Y 为横纵坐标，Φ(x,y)为电势。
3
金文璨 2007201230
(2) 二维等势线图：
�� (3) 电场满足 E = −∇φ ， 即电场沿电势降落的梯度方向， 故电场线垂直于等势线， 并由高电势指向低电势，如下图：
4
金文璨 2007201230
5
! 中间区域的初值，不妨设为 0
! 外壁的边值
! 金属棒的边值
DO k=1,10000 ! 迭代次数 DO i=2,7 DO j=2,20 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO DO i=8,20 DO j=2,7 u(i,j)=0.25*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)) ENDDO ENDDO
用差分法表示：
ui, j =
u i +1, j + u i −1, j + u i , j −1 + u i , j +1
4源自文库
（*）
u 0 , j = 0, u N , j = 0, u i ,0 = 0, u i , N = 0
u i , j = 1; x i ∈ [ − 0 .3, 0 .3 ], y i ∈ [ − 0 .3, 0 .3 ]
金文璨 2007201230
作业 6.1
一． 题目
有一个中空的无限长棱柱，其横截面如图所示，棱柱壁是金属的，棱柱中是金属 棒，在金属棒和外壁上分别加电压，计算金属棒和棱柱壁之间的电场，画出等势 线，电势和电场的分布。
二． 算法
金属棒为等势体，金属棒和外壁之间的区域内，电势满足 Laplace 方程。解电势
φ ( x, y ) 的 Laplace 边值问题：
∂ 2φ ∂ 2φ ∆φ = 2 + 2 = 0 ∂x ∂y
φ(x = ±1, y) = 0 ； φ ( x , y = ± 1) = 0 φ ( x , y ) = 1; x ∈ [ − 0 .3 .0 .3], y ∈ [ − 0 .3 .0 .3]
五． 讨论
(1) 迭代之前，金属棒和外壁之间区域的初值可以任意给定，但方便起见，不妨 设为零； (2) 该计算中，因为划分网格的步长比较短（h=0.1） ，格点比较少，所以迭代次 数不用太多（10000 次） ，即可得到稳定的结果。如果格点比较多的话，可以 对程序做一改进， 设定一个阈值， 当前后两次迭代的结果的差值小于阈值时 ， 迭代停止； (3) 由对称性分析可得，电势分布关于 x 轴、y 轴轴对称，同时也关于坐标原点 中心对称，所以如果格点比较多，计算量大的话，可以只计算第一象限的电 势值，然后通过坐标变换得到其他三个象限的电势值； (4) 从理论上分析，从金属棒到外壁，电势均匀降落，等势线应该是同心的圆角 正方形环，电场线应该是从金属棒向外壁辐射，并垂直于等势线； (5) Fortran 程序计算出电势的分布，在结果可视化的过程中，需要 MATLAB 对 电势进行处理，用 contour(x,y,u)语句绘制等势线，用 gradient(u)绘制电场矢 量的分布。