拉普拉斯方程的数值解
泊松方程与拉普拉斯
泊松方程与拉普拉斯泊松方程与拉普拉斯方程是数学领域中重要的偏微分方程，它们在物理学、工程学、计算机科学等各个领域有着广泛的应用。
本文将介绍泊松方程和拉普拉斯方程的定义、性质以及它们在实际问题中的应用。
泊松方程是一个二阶偏微分方程，通常用于描述电位、温度、流体静压力分布等问题。
其一般形式可以表示为：∆u = f(x,y,z)其中，u是待求函数，∆表示Laplace算子，f(x,y,z)是已知的函数。
泊松方程的求解过程包括确定边界条件、选择适当的解析方法等。
在一些特殊情况下，泊松方程可以通过分离变量、格林函数等方法精确求解。
拉普拉斯方程是泊松方程的特殊情况，即f(x,y,z)=0。
它表示了没有源项的稳定状态下的物理量分布。
例如，在无电荷的情况下，电势的分布可以由拉普拉斯方程描述。
泊松方程和拉普拉斯方程在实际问题中具有重要的应用。
下面将介绍它们在物理学、工程学和计算机科学中的具体应用。
一、物理学应用：1. 电场分布：根据泊松方程，可以求解电荷分布对电场的影响。
例如，在计算静电场、电容器以及电场中带电粒子的运动等问题时，泊松方程能够提供准确的分析结果。
2. 热传导问题：热传导是物体内部以及不同物体之间的热量传递过程。
泊松方程可以描述温度分布的稳定状态，因此可以求解热传导问题。
例如，在石油勘探中，泊松方程可用于分析地下温度场的分布。
二、工程学应用：1. 结构力学：泊松方程可用于模拟材料的弯曲、拉伸、压缩等受力状态。
例如，在工程结构设计中，可以利用泊松方程分析材料的变形和应力分布。
2. 流体力学：泊松方程可以用于模拟流体流动中的压力分布。
例如，在空气动力学中，可以用泊松方程求解空气流动的速度场和压力场。
三、计算机科学应用：1. 图像处理：在数字图像处理中，拉普拉斯算子可以用于图像边缘检测。
通过计算图像中像素灰度值的二阶导数，可以突出显示图像中的边缘结构。
2. 数值计算：泊松方程和拉普拉斯方程是数值计算领域中常用的方程之一。
有限差分法求解拉普拉斯方程
收稿日期:2009-09-01第一作者简介:贾新民(1956-),男,四川邻水人,新疆昌吉学院计算机工程系,副教授,研究方向:计算机程序设计及其语言教学和理论物理研究。
有限差分法求解拉普拉斯方程贾新民1 严文2(1.昌吉学院计算机工程系新疆昌吉831100;2.昌吉学院物理系 新疆昌吉831100)摘 要:以极板上具有半圆截面沟槽的电容器内的电势分布为例,介绍了综合应用计算机软件利用有限差分法求解复杂边界的拉普拉斯方程数值解的方法。
并利用数值解的结果讨论了沟槽表面的电场分布和电荷分布。
关键词:拉普拉斯方程;有限差分法;五点差分格式中图分类号:O411.2 文献标识码:A 文章编号:1671-6469(2009)05-0105-051 引言无源空间的引力场、静电场、稳定的温度分布等问题都满足拉普拉斯(Laplace )方程 2u (x ,y ,z )=0(1)但由于方程(1)是偏微分方程,只有在问题具有高度对称的情况下,才能求出解析解,而这种情形是极少的。
有些情形看上去很简单,但却求不出解析解。
对于这些情况,只能寻求数值解。
2 计算机数值解法方案文献[1][2][3]给出了拉普拉斯方程数值解的方法———有限差分法。
有限差分法的思想是用差分Δu (x +Δx,y )Δx ,Δu (x ,y +Δy )Δy 代替导数5u 5x ,5u 5y,用网格将求解区域覆盖,对于平面拉普拉斯方程,第i 行第j 列小格的电势由Laplace 方程的五点差分格式给出。
u ij =14(u ij -1+u ij +1+u i -1j +u i +1j )(2)考虑图1所示的具有半圆形截面的槽的电容器内部的电势和电场分布。
为了能够对坑(槽)内部的电场进行比较细致的观察,应该将半径R 取的大些,为了满足无限远的条件,应该使求解区域尽量大些。
我们选Excel 为计算工具,因为Excel 具有不用编写程序和直观的优点。
Excel 的一个单元格代表求解区域的一个网格,单元格的值表示该网格处的电势。
有限差分法解拉普拉斯方程python
有限差分法解拉普拉斯方程python一、引言拉普拉斯方程是一个重要的偏微分方程，它在数学、物理、工程等领域都有广泛的应用。
有限差分法是一种常用的数值求解方法，可以有效地解决拉普拉斯方程。
本文将介绍如何使用Python语言实现有限差分法求解拉普拉斯方程。
二、数学模型拉普拉斯方程可以表示为：∇²u = 0其中，u为未知函数，∇²表示Laplace算子。
在二维情况下，可以将该方程写成：∂²u/∂x² + ∂²u/∂y² = 0三、有限差分法有限差分法是一种常用的数值求解方法，在此不再赘述其原理和推导过程。
对于二维情况下的拉普拉斯方程，我们可以采用五点差分公式进行离散化处理：(u(i+1,j) - 2*u(i,j) + u(i-1,j))/Δx² + (u(i,j+1) - 2*u(i,j) + u(i,j-1))/Δy² = 0其中，Δx和Δy分别表示网格间距。
将上式变形可得：u(i,j) = (u(i+1,j) + u(i-1,j))*Δy² + (u(i,j+1) + u(i,j-1))*Δx² / (2*(Δx² + Δy²))四、Python实现在Python中，我们可以使用numpy库来处理数组和矩阵运算，使用matplotlib库来进行可视化。
首先，我们需要定义网格大小和间距：import numpy as npnx = 101 # 网格点数ny = 101dx = 2/(nx-1) # x方向间距dy = 2/(ny-1) # y方向间距接着，我们需要定义初始条件和边界条件：p = np.zeros((ny, nx)) # 初始条件# 边界条件p[:,0] = 0 # 左边界p[:,-1] = y # 右边界p[0,:] = p[1,:] # 下边界p[-1,:] = p[-2,:] # 上边界其中，左右边界分别为零和y的值，上下边界采用一阶差分法进行处理。
拉普拉斯方程在电场中的应用
拉普拉斯方程在电场中的应用拉普拉斯方程是一种常见的偏微分方程，它在数学和物理学中具有广泛的应用。
其中之一就是在电场中的应用。
在本文中，我们将探讨拉普拉斯方程在电场问题中的重要性和应用。
首先，让我们回顾一下拉普拉斯方程的定义。
拉普拉斯方程可以被表示为：$$\nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} = 0$$其中$\nabla^2$是拉普拉斯算子，$\phi$是电势函数。
这个方程描述了在没有电荷存在的区域内，电势函数满足的条件。
换句话说，该方程描述了一个没有电荷的区域内的平衡电势分布。
在电场问题中，拉普拉斯方程的应用主要体现在求解电势分布。
通过求解该方程，我们可以得到一个系统中各个点的电势。
这在研究电场中物体的电势分布和电场强度分布时非常有用。
一个经典的电场问题是带电导体的电势分布。
假设有一个带有电荷的导体，我们希望确定导体周围的电势分布。
通过求解拉普拉斯方程，我们可以获得关于电势分布的信息。
另一个常见的应用是在存在电荷分布的情况下求解电势分布。
如果一个系统中存在电荷密度$\rho$，我们可以使用拉普拉斯方程来求解电势分布并获得系统中各个点的电势。
这对于研究电场对带电粒子的作用以及电势能的计算都是非常重要的。
拉普拉斯方程的数值求解也是一种常见的方法。
在实际问题中，我们往往无法通过解析的方式求解拉普拉斯方程，而是需要使用数值方法来近似求解。
这可以通过离散化问题和应用数值计算方法来实现。
数值求解拉普拉斯方程在电场问题中具有很高的实用性和广泛的应用。
除了电势分布的计算，拉普拉斯方程还可以用于求解电场导致的其他物理量。
例如，根据拉普拉斯方程，我们可以计算电场强度分布、电势能分布以及电荷运动的轨迹等等。
拉普拉斯方程在柱坐标下的表达式
拉普拉斯方程在柱坐标下的表达式在物理学和工程学中，拉普拉斯方程是一个重要的偏微分方程，用于描述场的分布情况。
在柱坐标系下，拉普拉斯方程的表达式具有特定的形式，下面我们来详细介绍。
柱坐标系简介柱坐标系是三维空间中的一种常见坐标系，通常用$(r, \\theta, z)$表示。
其中r表示径向距离，$\\theta$表示与z轴的夹角，z表示沿z轴方向的距离。
拉普拉斯方程在柱坐标系下的表达式在柱坐标系下，拉普拉斯方程可以表示为：$$\ abla^2 \\phi = \\frac{1}{r} \\frac{\\partial}{\\partial r} \\left(r\\frac{\\partial \\phi}{\\partial r}\\right) + \\frac{1}{r^2} \\frac{\\partial^2\\phi}{\\partial \\theta^2} + \\frac{\\partial^2 \\phi}{\\partial z^2} = 0$$ 这里abla2是拉普拉斯算子，$\\phi(r, \\theta, z)$是待求解的场函数。
拉普拉斯方程的物理意义拉普拉斯方程描述了场在空间中的分布情况，其解决了场在各个方向上的二阶导数之和为零的情况。
在物理学中，拉普拉斯方程常常与介电场、温度场等相关联，求解拉普拉斯方程可以帮助我们理解场在空间中的行为。
柱坐标系下的拉普拉斯方程求解为了求解柱坐标系下的拉普拉斯方程，通常需要将方程进行分离变量，将$\\phi(r, \\theta, z)$表示为$R(r) \\Theta(\\theta) Z(z)$的形式，然后分别对$r,\\theta, z$方向的分量进行求解。
经过适当的分离变量和边界条件的设定，可以得到拉普拉斯方程在柱坐标系下的解析解或数值解。
这将帮助我们更好地理解在柱坐标系下场的分布情况。
总结拉普拉斯方程在柱坐标系下的表达式为$\ abla^2 \\phi = \\frac{1}{r}\\frac{\\partial}{\\partial r} \\left(r \\frac{\\partial \\phi}{\\partial r}\\right) +\\frac{1}{r^2} \\frac{\\partial^2 \\phi}{\\partial \\theta^2} + \\frac{\\partial^2\\phi}{\\partial z^2} = 0$，它描述了场在空间中的分布情况。
拉普拉斯方程
拉普拉斯方程拉普拉斯方程(Laplace equation)拉普拉斯方程表示液面曲率与液体压力之间的关系的公式。
一个弯曲的表面称为曲面，通常用相应的两个曲率半径来描述曲面，即在曲面上某点作垂直于表面的直线，再通过此线作一平面，此平面与曲面的截线为曲线，在该点与曲线相重合的圆半径称为该曲线的曲率半径R1。
通过表面垂线并垂直于第一个平面再作第二个平面并与曲面相交，可得到第二条截线和它的曲率半径R2，用 R1与R2可表示出液体表面的弯曲情况。
若液面是弯曲的，液体内部的压力p1与液体外的压力p2就会不同，在液面两边就会产生压力差?P= P1- P2，其数值与液面曲率大小有关，可表示为:在数理方程中，拉普拉斯方程为:?u=d^2u/dx^2+d^2u/dy^2=0，其中?为拉普拉斯算子，此处的拉普拉斯方程为二阶偏微分方程。
三维情况下，拉普拉斯方程可由下面的形式描述，问题归结为求解对实自变量 x 、 y 、 z 二阶可微的实函数φ :上面的方程常常简写作:或其中div表示矢量场的散度(结果是一个标量场)，grad表示标量场的梯度(结果是一个矢量场)，或者简写作:其中Δ称为拉普拉斯算子 .拉普拉斯方程的解称为调和函数。
如果等号右边是一个给定的函数 f ( x , y , z )，即:则该方程称为泊松方程。
拉普拉斯方程和泊松方程是最简单的椭圆型偏微分方程。
偏微分算子或Δ(可以在任意维空间中定义这样的算子)称为拉普拉斯算子，英文是 Laplace operator 或简称作 Laplacian 。
拉普拉斯方程的狄利克雷问题可归结为求解在区域 D 内定义的函数φ，使得在 D 的边界上等于某给定的函数。
为方便叙述，以下采用拉普拉斯算子应用的其中一个例子——热传导问题作为背景进行介绍:固定区域边界上的温度(是边界上各点位置坐标的函数)，直到区域内部热传导使温度分布达到稳定，这个温度分布场就是相应的狄利克雷问题的解。
拉普拉斯方程的诺伊曼边界条件不直接给出区域 D 边界处的温度函数φ本身，而是φ沿 D 的边界法向的导数。
拉普拉斯方程数值解与解析解的研究
拉普拉斯方程数值解与解析解的研究拉普拉斯方程是一个二阶偏微分方程，通常用于描述稳态的物理现象，如电势、温度分布等。
其解析解在大多数情况下难以求出，因此需要使用数值方法解决。
本文将介绍拉普拉斯方程数值解与解析解的研究。
首先，我们需要了解拉普拉斯方程的基本形式。
对于一个二维的拉普拉斯方程，其一般形式为：∇²u(x, y) = 0其中，u(x, y)表示二维平面上的一个标量场，∇²表示拉普拉斯算子，它是二阶偏导数的和。
接下来，我们需要了解解析解和数值解的概念。
解析解是指能够用解析表达式表示的解，通常是利用数学方法求解得到的结果。
数值解是指通过数值计算方法获得的近似解，通常是使用计算机进行求解。
对于拉普拉斯方程的解析解，通常只有在特定的边界条件下才能求得。
例如，对于一个具有矩形边界的区域，可以使用分离变量法求解得到解析解。
但是对于更加复杂的边界条件或者更加复杂的区域形状，通常需要使用数值方法进行求解。
常用的拉普拉斯方程数值解方法有迭代法、差分法和有限元法等。
其中，迭代法是一种基于逐步逼近的求解方法，通常需要指定初始值和逼近精度来进行求解。
差分法是一种基于求解离散方程的方法，通常将连续的区域离散化为离散的网格，然后利用差分近似来求解方程。
有限元法则是一种将连续区域分割成许多小的有限单元，然后在每个单元上使用局部的近似函数来求解方程的方法。
在进行拉普拉斯方程数值解时，我们需要注意选择合适的方法以及合适的边界条件。
同时，我们也需要考虑数值误差的影响，通常需要进行误差分析来评估计算结果的准确性。
综上所述，拉普拉斯方程的解析解在大多数情况下难以求得，而数值解则是一种常用的求解方法。
在进行数值解时，需要选择合适的方法和边界条件，并考虑数值误差的影响。
偏微分方程的数值解法
偏微分方程的数值解法1.peEllip5 用五点差分格式解拉普拉斯方程function u = peEllip5(nx,minx,maxx,ny,miny,maxy) format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);u0 = zeros(nx,ny);for j=1:nyu0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);endfor j=1:nxu0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);endA = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));b = zeros((nx-2)*(ny-2),1);for i=1:(nx-2)*(ny-2)if mod(i,nx-2) == 1if i==1A(1,2) = 1;A(1,nx-1) = 1;b(1) = - u0(1,2) - u0(2,1);elseif i == (ny-3)*(nx-2)+1A(i,i+1) = 1;A(i,i-nx+2) = 1;b(i) = - u0(ny-1,1) - u0(ny,2);elseA(i,i+1) = 1;A(i,i-nx+2) = 1;A(i,i+nx-2) = 1;b(i) = - u0(floor(i/(nx-2))+2,1);endendelseif mod(i,nx-2) == 0if i == nx-2A(i,i-1) = 1;A(i,i+nx-2) = 1;b(i) = - u0(1,nx-1) - u0(2,nx);elseif i == (ny-2)*(nx-2)A(i,i-1) = 1;A(i,i-nx+2) = 1;b(i) = - u0(ny-1,nx) - u0(ny,nx-1);elseA(i,i-1) = 1;A(i,i-nx+2) = 1;A(i,i+nx-2) = 1;b(i) = - u0(floor(i/(nx-2))+1,nx);endendelseif i>1 && i< nx-2A(i,i-1) = 1;A(i,i+nx-2) = 1;A(i,i+1) = 1;b(i) = - u0(1,i+1);elseif i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)A(i,i-1) = 1;A(i,i-nx+2) = 1;A(i,i+1) = 1;b(i) = - u0(ny,mod(i,(nx-2))+1);elseA(i,i-1) = 1;A(i,i+1) = 1;A(i,i+nx-2) = 1;A(i,i-nx+2) = 1;endendendendendul = A\b;for i=1:(ny-2)for j=1:(nx-2)u(i,j) = ul((i-1)*(nx-2)+j);endendformat short;2.peEllip5m 用工字型差分格式解拉普拉斯方程function u = peEllip5m(nx,minx,maxx,ny,miny,maxy) format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);u0 = zeros(nx,ny);for j=1:nyu0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);endfor j=1:nxu0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);endA = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));b = zeros((nx-2)*(ny-2),1);for i=1:(nx-2)*(ny-2)if mod(i,nx-2) == 1if i==1A(1,nx) = 1;b(1) = - u0(1,1) - u0(3,1) - u0(1,3);elseif i == (ny-3)*(nx-2)+1A(i,i-nx+1) = 1;b(i) = - u0(ny,1) - u0(ny,3) - u0(ny-2,1);elseA(i,i-nx+3) = 1;A(i,i+nx-1) = 1;b(i) = - u0(floor(i/(nx-2))+1,1) - u0(floor(i/(nx-2))+3,1);endendelseif mod(i,nx-2) == 0if i == nx-2A(i,i+nx-1) = 1;b(i) = - u0(1,nx-2) - u0(1,nx) - u0(3,nx);elseif i == (ny-2)*(nx-2)A(i,i-nx+1) = 1;b(i) = - u0(ny,nx) - u0(ny,nx-2) - u0(ny-2,nx);elseA(i,i-nx+1) = 1;A(i,i+nx-3) = 1;b(i) = - u0(floor(i/(nx-2)),nx) - u0(floor(i/(nx-2))+2,nx);endendelseif i>1 && i< nx-2A(i,i+nx-1) = 1;A(i,i+nx-3) = 1;b(i) = - u0(1,i) - u0(1,i+2);elseif i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)A(i,i-nx+3) = 1;A(i,i-nx+1) = 1;b(i) = - u0(ny,mod(i,(nx-2))) -u0(ny,mod(i,(nx-2))+2);elseA(i,i-nx+3) = 1;A(i,i+nx-1) = 1;A(i,i+nx-3) = 1;A(i,i-nx+1) = 1;endendendendendul = A\b;for i=1:(ny-2)for j=1:(nx-2)u(i,j) = ul((i-1)*(nx-2)+j);endendformat short;3.peHypbYF 用迎风格式解对流方程function u = peYF(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);if a>0for j=1:(n+M)u0(j) = IniU(minx+(j-M-1)*h);endelsefor j=1:(n+M)u0(j) = IniU(minx+(j-1)*h);endendu1 = u0;for k=1:Mif a>0for i=(k+1):n+Mu1(i) = -dt*a*(u0(i)-u0(i-1))/h+u0(i);endelsefor i=1:n+M-ku1(i) = -dt*a*(u0(i+1)-u0(i))/h+u0(i);endendu0 = u1;endif a>0u = u1((M+1):M+n);elseu = u1(1:n);endformat long;4.peHypbLax 用拉克斯-弗里德里希斯格式解对流方程function u = peHypbLax(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ku1(i) = -dt*a*(u0(i+1)-u0(i-1))/h/2+(u0(i+1)+u0(i-1))/2;endu0 = u1;endu = u1((M+1):(M+n));format short;4.peHypbLaxW 用拉克斯-温德洛夫格式解对流方程function u = peLaxW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ku1(i) = dt*dt*a*a*(u0(i+1)-2*u0(i)+u0(i-1))/2/h/h - ...dt*a*(u0(i+1)-u0(i-1))/h/2+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;6.peHypbBW 用比姆-沃明格式解对流方程function u = peBW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-2*M-1)*h);endu1 = u0;for k=1:Mfor i=2*k+1:n+2*Mu1(i) = u0(i)-dt*a*(u0(i)-u0(i-1))/h-a*dt*(1-a*dt/h)* ...(u0(i)-2*u0(i-1)+u0(i-2))/2/h;endu0 = u1;endu = u1((2*M+1):(2*M+n));format short;7.peHypbRich 用Richtmyer多步格式解对流方程function u = peRich(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+4*M)u0(j) = IniU(minx+(j-2*M-1)*h);endfor k=1:Mfor i=2*k+1:n+4*M-2*ktmpU1 = -dt*a*(u0(i+2)-u0(i))/h/4+(u0(i+2)+u0(i))/2;tmpU2 = -dt*a*(u0(i)-u0(i-2))/h/4+(u0(i)+u0(i-2))/2;u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+u0(i);endu0 = u1;endu = u1((2*M+1):(2*M+n));format short;8.peHypbMLW 用拉克斯-温德洛夫多步格式解对流方程function u = peMLW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i))/h/2+(u0(i+1)+u0(i))/2;tmpU2 = -dt*a*(u0(i)-u0(i-1))/h/2+(u0(i)+u0(i-1))/2;u1(i) = -dt*a*(tmpU1-tmpU2)/h+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;9.peHypbMC 用MacCormack多步格式解对流方程function u = peMC(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endfor k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i))/h+u0(i);tmpU2 = -dt*a*(u0(i)-u0(i-1))/h+u0(i-1);u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+(u0(i)+tmpU1)/2;endu0 = u1;endu = u1((M+1):(M+n));format short;10.peHypb2LF 用拉克斯-弗里德里希斯格式解二维对流方程的初值问题function u = pe2LF(a,b,dt,nx,minx,maxx,ny,miny,maxy,M)%啦-佛format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);for i=1:nx+2*Mfor j=1:(ny+2*M)u0(i,j) = Ini2U(minx+(i-M-1)*hx,miny+(j-M-1)*hy);endendu1 = u0;for k=1:Mfor i=k+1:nx+2*M-kfor j=k+1:ny+2*M-ku1(i,j) = (u0(i+1,j)+u0(i-1,j)+u0(i,j+1)+u0(i,j-1))/4 ...-a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx ...-b*dt*(u0(i,j+1)-u0(i,j-1))/2/hy;endendu0 = u1;endu = u1((M+1):(M+nx),(M+1):(M+ny));format short;11.peHypb2FL 用拉克斯-弗里德里希斯格式解二维对流方程的初值问题function u = pe2FL(a,b,dt,nx,minx,maxx,ny,miny,maxy,M)format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);for i=1:nx+4*Mfor j=1:(ny+4*M)u0(i,j) = Ini2U(minx+(i-2*M-1)*hx,miny+(j-2*M-1)*hy);endendu1 = u0;for k=1:Mfor i=2*k+1:nx+4*M-2*kfor j=2*k-1:ny+4*M-2*k+2tmpU(i,j) = u0(i,j) - a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx + ...(a*dt/hx)^2*(u0(i+1,j)-2*u0(i,j)+u0(i-1,j))/2;endendfor i=2*k+1:nx+4*M-2*kfor j=2*k+1:nx+4*M-2*ku1(i,j) = tmpU(i,j) - b*dt*(tmpU(i,j+1)-tmpU(i,j-1))/2/hy + ...(b*dt/hy)^2*(tmpU(i,j+1)-2*tmpU(i,j)+tmpU(i,j-1))/2;endendu0 = u1;endu = u1((2*M+1):(2*M+nx),(2*M+1):(2*M+ny));format short;12.peParabExp 用显式格式解扩散方程的初值问题function u = peParabExp(c,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = PrIniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ku1(i) = dt*c*(u0(i+1)-2*u0(i)+u0(i-1))/h/h+u0(i);endendu = u1((M+1):(M+n));format short;13.peParabTD 用跳点格式解扩散方程的初值问题function u = peParabTD(c,dt,n,minx,maxx,M)%跳点format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = PrIniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-kif mod(n+i,2) == 0u1(i) = u0(i) + c*dt*(u0(i+1) - 2*u0(i) + u0(i-1))/h/h;if i > 2u1(i-1) =(u0(i-1) + c*dt*(u1(i) + u1(i-2))/h/h)/(1 + 2*c*dt/h/h);endendendu0 = u1;endu = u1((M+1):(M+n));format short;14.peParabImp 用隐式格式解扩散方程的初边值问题function u = peParabImp(c,dt,n,minx,maxx,lbu,rbu,M)format long;h = (maxx-minx)/(n-1);u0(1) = lbu;u0(n) = rbu;for j=2:n-1u0(j) = PrIniU(minx+(j-1)*h);endu1 = u0;for k=1:MA = zeros(n-2,n-2);cb = - transpose(u0(2:(n-1)));cb(1) = cb(1) - dt*c*lbu/h/h;cb(n-2) = cb(n-2) - dt*c*rbu/h/h;A(1,1) = -2*dt*c/h/h -1;A(1,2) = dt*c/h/h ;for i=2:n-3A(i,i-1) = dt*c/h/h ;A(i,i) = - 2*dt*c/h/h -1 ;A(i,i+1) = dt*c/h/h ;endA(n-2,n-2) = -2*dt*c/h/h -1;A(n-2,n-3) = dt*c/h/h;u1(2:(n-1)) = A\cb;u0 = u1;endu = u1;format short;15.peParabKN 用克拉克-尼科尔森格式解扩散方程的初边值问题function u = peParabKN(c,dt,n,minx,maxx,lbu,rbu,M)format long;h = (maxx-minx)/(n-1);u0(1) = lbu;u0(n) = rbu;for j=2:n-1u0(j) = PrIniU(minx+(j-1)*h);endu1 = u0;for k=1:MA = zeros(n-2,n-2);cb = zeros(n-2,1);cb(1) = -u0(2) -(u0(3)-2*u0(2)+u0(1))*dt*c*lbu/h/h/2 - dt*c*lbu/h/h/2;cb(n-2) = -u0(n-1) -(u0(n)-2*u0(n-1)+u0(n-2))*dt*c*lbu/h/h/2 - dt*c*rbu/h/h/2;for i=2:n-3cb(i) = -u0(i+1) -(u0(i+2)-2*u0(i+1)+u0(i))*dt*c*lbu/h/h/2;endA(1,1) = -dt*c/h/h -1;A(1,2) = dt*c/h/h/2 ;for i=2:n-3A(i,i-1) = dt*c/h/h/2 ;A(i,i) = - dt*c/h/h -1 ;A(i,i+1) = dt*c/h/h/2 ;endA(n-2,n-2) = -dt*c/h/h -1;A(n-2,n-3) = dt*c/h/h/2;u1(2:(n-1)) = A\cb;u0 = u1;endu = u1;format short;16.peParabWegImp 用加权隐式格式解扩散方程的初边值问题function u = peParabWegImp(c,sita,dt,n,minx,maxx,lbu,rbu,M)format long;h = (maxx-minx)/(n-1);u0(1) = lbu;u0(n) = rbu;for j=2:n-1u0(j) = PrIniU(minx+(j-1)*h);endu1 = u0;for k=1:MA = zeros(n-2,n-2);cb = zeros(n-2,1);cb(1) = -u0(2) -(1 - sita)*(u0(3)-2*u0(2)+u0(1))*dt*c*lbu/h/h/2 ...- sita*dt*c*lbu/h/h;cb(n-2) = -u0(n-1) -(1 - sita)*(u0(n)-2*u0(n-1)+u0(n-2))*dt*c*lbu/h/h/2 ...- sita*dt*c*rbu/h/h;for i=2:n-3cb(i) = -u0(i+1) -(1 - sita)*(u0(i+2)-2*u0(i+1)+u0(i))*dt*c*lbu/h/h;endA(1,1) = -2*sita*dt*c/h/h -1;A(1,2) = sita*dt*c/h/h ;for i=2:n-3A(i,i-1) = sita*dt*c/h/h ;A(i,i) = - 2*sita*dt*c/h/h -1 ;A(i,i+1) = sita*dt*c/h/h ;endA(n-2,n-2) = - 2*sita*dt*c/h/h -1;A(n-2,n-3) = sita*dt*c/h/h;u1(2:(n-1)) = A\cb;u0 = u1;endu = u1;format short;17.peDKExp 用指数型格式解对流扩散方程的初值问题function u = peDKExp(a,b,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = DKIniU(minx+(j-M-1)*h);endu1 = u0;coff = (exp(a*h/2/b)+exp(-a*h/2/b))/(exp(a*h/2/b)-exp(-a*h/2/b));coff = dt*coff*a*h/2;for k=1:Mfor i=k+1:n+2*M-ku1(i) = coff*(u0(i+1)-2*u0(i)+u0(i-1))+a*dt*(u0(i+1)-u0(i-1))/h/2+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;18.peDKSam 用萨马尔斯基格式解对流扩散方程的初值问题function u = peDKSam(a,b,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = DKIniU(minx+(j-M-1)*h);endu1 = u0;coff = dt*b/(1+a*h/b/2);for k=1:Mfor i=k+1:n+2*M-ku1(i) = coff*(u0(i+1)-2*u0(i)+u0(i-1))+a*dt*(u0(i)-u0(i-1))/h+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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
φ ( 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]
! 中间区域的初值，不妨设为 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
5
四． 计算结果可视化
(1) 电势的分布图：3D 图中，X,Y 为横纵坐标，Φ(x,y)为电势。
3
金文璨 2007201230
(2) 二维等势线图：
�� (3) 电场满足 E = −∇φ ， 即电场沿电势降落的梯度方向， 故电场线垂直于等势线， 并由高电势指向低电势，如下图：
4
金文璨 2007201230
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 Eu 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 方程。解电势
五． 讨论
(1) 迭代之前，金属棒和外壁之间区域的初值可以任意给定，但方便起见，不妨 设为零； (2) 该计算中，因为划分网格的步长比较短（h=0.1） ，格点比较少，所以迭代次 数不用太多（10000 次） ，即可得到稳定的结果。如果格点比较多的话，可以 对程序做一改进， 设定一个阈值， 当前后两次迭代的结果的差值小于阈值时 ， 迭代停止； (3) 由对称性分析可得，电势分布关于 x 轴、y 轴轴对称，同时也关于坐标原点 中心对称，所以如果格点比较多，计算量大的话，可以只计算第一象限的电 势值，然后通过坐标变换得到其他三个象限的电势值； (4) 从理论上分析，从金属棒到外壁，电势均匀降落，等势线应该是同心的圆角 正方形环，电场线应该是从金属棒向外壁辐射，并垂直于等势线； (5) Fortran 程序计算出电势的分布，在结果可视化的过程中，需要 MATLAB 对 电势进行处理，用 contour(x,y,u)语句绘制等势线，用 gradient(u)绘制电场矢 量的分布。