偏微分方程数值解
8.Solve菜单 Solve PDE 对当前的几何结构实体 CSC、三角形网格和图形解偏微分方 程。
Parameters…打开PDE对话框,输入解 PDE的参数。 Expeort Solution… 输出PDE方程的解 矢量u。如果可行,将计算特征值1输 出到主工作区间。
Parameters…打开对话框,可以输入解
在MATLAB的命令窗口中输入pdetool命令, 然后单击回车键,就将显示PDE图形界面,如图:
PDE Toolbox 菜单
1.File 菜单 New 更新或建立一个新的几何结果实体模型 Constructive Soild Geometry(CSG), 并取名为 Untitled。 Open … 从硬盘装载M文件。 Save 将在GUI内完成的成果储存到一个M文件 中。
PDE Specifficatiom…打开一个对话框,输入偏微分方程
类型和应用参数。参数的维数决定于偏微分方程的维数。
如果选择专业应用模式,那么特殊偏微分方程和参数将 代替标准偏微分方程系数。每一个参数c、a、f和d皆可 作为有效的MATLAB表达式,以作为计算三角形单元质 量中心的参数值。下面的变量是很有用的。
二阶线性偏微分方程的分类
自变量多于一个微分方程称为偏微分方程。 一般的二阶线性方程总可以写成如下的形状:
若二阶偏导数项的系数 满足: 则称方程在点
则称方程在点
则称方程在点
在区域G中某点
为双曲型的; 为抛物型的; 为椭圆型的。
二阶偏微分方程的解法的方法主要有两大类,其中 每一类又包含若干种方法。
第一类是解析法 分为分离变量法、保角变换法、镜像法和格林函 数法等。
前处理
网格剖分
取定沿x轴和y轴方向的步长 和 ,
。作
两族与坐标轴平行的直线:
, y ih2 , j 0,1, 10
两族直线的交点 ih1, jh2 称为网点或节点,记为 i, j
有限差分的网格分割
方程离散
在直角坐标系中,矩形槽中的电位函数满 足拉普拉斯方程,即
根据Taylor展开式 ,有
v1(1,:)=zeros(1,hx); %左右两列的Dirichlet边界条件
for i=1:hy
v1(i,1)=0;
v1(i,hx)=0;
End v2=v1;maxt=1;t=0; %初始化 k=0
while(maxt>1e-6) %由v1迭代,算出v2,迭代精度为
k=k+1
%计算迭代次数
0.000001
式中: ——加速收敛因子
迭代收敛的速度与 有明显关系:
收敛因子( ) 1.0 1.7 1.8 1.83 1.85 1.87 1.9 2.0 迭代次数( N) >1000 269 174 143 122 133 171 发散
最佳收敛因子的经验公式: 正方形场域、正方形网格
矩形场域、正方形网格
迭代收敛的速度与电位初始值的给定及网格剖分精 细有关; 迭代收敛的速度与工程精度要求有关。
end
end
v1=v2
end
subplot(1,2,1),mesh(v2) %画三维曲面 axis([0,10,0,10,0,100]) subplot(1,2,2),contour(v2,15) %画等电位线 hold on x=1:1:hx;y=1:1:hy [xx,yy]=meshgrid(x,y); [Gx,Gy]=gradient(v2,0.6,0.6); axis([-1.5,hx+2.5,-2,13]) plot([1,1,hx,hx,1],[1,hy,hy,1,1],'k') text(hx/2,0.3,'OV','fontsize',11); text(hx/2-0.5,hy+0.5',11); text(-0.5,hy/2,'OV','fontsize',11); text(hx+0.3,hy/2,'OV','fontsize',11); hold off
方程参数。每组参数的选择取决于 PDE的类型。
建立几何模型 设定边界条件 进行三角剖分
用PDE工具箱进行问题求解
x, y yb 100
x, y x0 0
x, y xa 0
x, y y0 0
应用PDE工具箱求解结 果如下图2:
2 xi1, y j 2 xi , y j xi1, y j
x 2
h12
2 xi , y j1 2 xi , y j xi , y j1
y 2
h22
2 2 xi1, y j
①x和y:x和y坐标; ②u:解; ③ux,uy:解的关于x和y的导数; ④t:时间。
7.Mesh 菜单 Mesh Mode 输入网格模式。 Initialize Mesh 建立和显示初始化三角形网格。 Refine Mesh 加密当前三角形网格。 Jiggle Mesh 优化网格。 Undo Mesh Change 退回上一次网格操作。 Display Triangle Quality 用0~1之间数字化颜色显 示三角形网格的质量,大于0.6的网格是可接受的。 Show Mode Labels 显示网格节点标识开关,三角形 网格标识数据是三角形矩阵t的列。 Export Mesh… 输出节点矩阵p,边界矩阵e和三角 形矩阵t到主工作区间。 Parameters… 打开对话框,修改网格生成参数。
Clear 删除已选的实体。 Select All 选择当前几何结构实体造 型CSG中的所有的实体及其边界和子
域。
3.Optionst 菜单 Grid 绘图时栅格的开启和关闭。 Grid Spacing… 调整栅格大小。 Snap 捕捉栅格开启和关闭。 Axis Limits… 改变绘图轴的比例。 Axis Equal 绘图轴打开和关闭。 Turn off Toolbar 关闭工具栏按钮的帮助文档。 Help 帮助 Zoom 图形缩放的开启和关闭。 Application 应用模式选择。 Refresh 重新显示PDE工具箱中所有的图形实体。
借助计算机进行计算时,其程序框图如下:
N
迭代解程序框图
启动 赋边界节点已知电位值 赋予场域内各节点电位初始值
累计迭代次数 N=0
N=N+1
按超松弛法进行一
次迭代,求
(N i, j
1)
Y
所有内点 相邻二次迭代值的最大误差
是否小于W
打印 N,(i, j) 停机
算法的实现
算法: hx=10;hy=10; %设置网格节点数 v1=ones(hy,hx); %设置行列二维数组 %上下两行的Dirichlet边界条件 v1(hy,:)=ones(1,hx)*100;
如果待求N个点的电位,就需解含有N 个方程的线性方程组。若点的数目较 多,用迭代法较为方便。 1.高斯—赛德尔迭代法
高斯——赛德尔迭代法
式中
节点迭代顺序一般按自然顺序排列,即先从左到右,在 从下到上。
迭代过程遇到边界节点时,代入边界值或边界差分
格式,直到所有节点电位满足
为止。
2、逐次超松弛迭代法 高斯—赛德尔迭代法收敛慢,为了加速收敛,把所得的结果 依次代入进行计算的同时,还进行迭代变化量的加权,有:
Ellipse/Circle(centered) 以中心方式画椭圆/圆(ctrl+ 鼠标)。
Polygon 画多边形,按右键可封闭多边形。 Rotate… 旋转已选的物体。 Export Geometry Description,Set Formula,Labels… 将 几何描述矩阵sd、公式设置字符sf和标识空间矩阵nk 输出到主工作空间去。
第二类是数值法 有限差分法(FDM)和有限单元(FEM)等。 其中有限差分法最为成熟的.在有限差分法中, 通常以差分代替微分,用有限求和代替积分,这 样,就将问题转化为求解差分方程或代数方程问 题。
以椭圆型方程为例介绍有限差分法
例:横截面为矩形的无限长槽由三块接地 导体板构成,如图所示,槽的盖板接直流 电压100V,求矩形槽的电位分布。
2 xi , y j
xi1, y j
xi , y j1
2 xi , y j
xi , y j1
x2 y 2
h12
h22
特别取正方形网格:h1 h2 h 则有:
二维场域的边界条件
区域S的边界为曲线l所包围,若二维场域边界上
的电位
, 是一个事先已知的电位函数,则
4.Draw菜单 Draw Mode 进入绘图模式。 Rectangle/Square 以角点方式画矩形/方形(ctrl+鼠 标)。
Rectangle/Square(centered) 以中心方式画矩形/方形 (ctrl+鼠标)。 Ellipse/Circle 以矩形角点方式画椭圆/圆(ctrl+鼠 标)。
Remove Border Subdomain 当图形进行布尔计算时,删除 已选取的子域边界。
Remove AllSubdomain Borders 当图形进行布尔计算时, 删除所有的子域边界。
Export Decomposed Geometry, Boundary Cond’s...将分解 几何矩阵g和边界条件矩阵b,输出到主工作区间。 Specify Boundary Conditions... 定义边界条件。从显示的对话 框,对已选的边界可输入边界条件,共有三种不同的条件类型。
6.PDE菜单 PDE Mode 进入偏微分方程模式。 Show Subdomain Labels 显示子区域标识开关。 PDE Specification…打开对话框以输入PDE参数和类型 Export PDE Coefficients…将当前PDE参数c,a,f和d 输出到主工作区间,其参数变量为字符类型。