MATLAB求解PDE问题.docx
- 文档编号:27549920
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:12
- 大小:396.72KB
MATLAB求解PDE问题.docx
《MATLAB求解PDE问题.docx》由会员分享,可在线阅读,更多相关《MATLAB求解PDE问题.docx(12页珍藏版)》请在冰豆网上搜索。
MATLAB求解PDE问题
MATLAB求解PDE问题
(1)-—概述、例子(转)
(2011—07-2016:
48:
45)
MATLABPDEToolbox提供利用有限元方法求解偏微分方程的GUI以及相应的命令行函数。
利用该工具箱可以求解椭圆型方程、抛物型方程、双曲型方程、特征值方程以及非线性方程.PDEToolbox的功能非常强大,网上有许多利用PDEToolbox解决各种物理问题的论文,还有专门介绍工具箱的参考书。
网上的例子虽然很多,但是大部分是介绍PDE工具箱自带的一些例子,这些例子中解的区域,边界条件是PDE工具箱已经编写好的,直接调用就可以。
对于该如何自己设定求解区域及边界条件,却很少有人涉及.网上搜索发现只有刘平在博客中详细介绍过求解区域的设定。
下面以一个椭圆型方程的例子来详细说明求解的各个步骤,希望对大家能有所帮助。
设要求如下形式的椭圆方程的解:
按照PDE的要求,将方程化为标准形式
求解后的图像如下,第一幅图是解的图像,第二幅是计算误差。
从第二幅图可以看到,计算的最大误差是10—3方量级。
通过这个例子我们可以基本掌握PDE求解偏微分方程的步骤和方法,后面我将详细介绍如何设置区域及边界条件。
掌握了区域和边界条件的设定,就可以轻松求解遇到的偏微分方程了。
图后是附带的matlab命令以及注释,并提供m文件附件下载,下载后解压即可.希望能对大家有所帮助。
下面是编写的求解上述方程的matlab语句及说明:
g='mygeom';
b=’mybound';
定义区域,边界条件。
mygeom是定义区域的子函数名,函数名可根据自己的需要取定,区域的确定规则由pdegeom函数说明,注意pdegeom函数只是说明如何定义区域,它并不直接确定区域;mybound是定义边界条件的子函数名,与区域类似,边界的确定规则由函数pdebound确定。
后面我会详细介绍区域和边界的取法。
[p,e,t]=initmesh(g);
网格初始化,此处也可以写成[p,e,t]=initmesh(’mygeom');这样可以省略上面的语句
[p,e,t]=refinemesh(g,p,e,t);
[p,e,t]=refinemesh(g,p,e,t);
加密网格两次,需要加密几次重复几次即可,根据具体问题确定加密次数
U=assempde(b,p,e,t,1,0,'2*(x+y)-4');
调用assempde函数计算方程的数值解,assempde函数的详细用法可以参考MATH网站或者PDE的使用指南。
常用的用法是[u,res]=assempde(b,p,e,t,c,a,f),其中b为边界条件,此处也可以写为’mybound’,p,e,t,为网格参数,c,a,f,为方程的参数,后面也可以加猜测值以及各种属性。
pdesurf(p,t,U)
gridon;
xlabel(’x');ylabel(’y');zlabel('u’)
colorbar
view([6030])
画出解的图形.注意,为了让结果更直观一些,使用view函数调整了视点位置。
大家可以自行调整视角,满意即可.
exact=p(1,:
)。
^2+p(2,:
).^2-p(1,:
)。
*p(2,:
).*(p(1,:
)+p(2,:
));
exact=exact';
figure
pdesurf(p,t,U—exact)
gridon
xlabel('x’);ylabel(’y');zlabel('error’)
colorbar
view([6030])
由于方程有解析解,我们可以比较数值计算的误差.如果能求得解析解,我们也不会设计各种方法求数值解了,因此,这一步在大多数情况下是用不上的,这里只是为了比较计算结果,验证计算的精度。
MATLAB求解PDE问题
(2)—-确定几何区域
(2012—06—1416:
20:
38)
前一篇介绍了如何利用Matlab求解椭圆型方程,下面介绍如何确定求解的几何区域。
PDEToolbox中规定几何区域的m文件是pdegeom.m。
但是pdegeom并不是一个可以调用的函数,它只是规定了应该何如定义区域,具体的区域则要根据研究的问题来决定。
函数pdegeom释义如下:
参数为0个时,即没有参数时,返回边界的段数;
参数为1个时,即只有bs,返回输出区域边界的参变量范围矩阵d;
参数为2个时,返回每段边界长度为s时的坐标.
函数参数意义bs表示指定的边缘线段,如矩形边界为四段,三角开边界肯定为三段…。
s为第bs段线段弧长的近似(估计)值,bs与s可以为向量,但是要一一对应,即bs为几个值,s也得为几个值.输出变量[x,y]是每条线段起点和终点所对应的坐标。
这个函数编制的关键是,函数内边界上的坐标((x(t),y(t))是用参变量t表示的,返回值是求得边界任意长度时的坐标(x(t),y(t))值,参量可以有很多种选法.
回到前一篇中,给定的方程的求解区域是[0,1,0,1]的一个正方形,我们将它命名为mygeom。
下面我们来看下mygeom是怎么编写的。
function[x,y]=mygeom(bs,s)
nbs=4;%表示边界的段数
ifnargin==0,
x=nbs;%不给定输入变量时,输出表示几何区域边界的线段数
return
end
d=[
0000%表示每条线段起始值的参量值t(四条边界,所以有四列)
1111%表示每条线段的终点值参量值t
0000%沿线段方向左边区域的标识值,如果是1,表示选定左边区域,如果
%是0,表示不选定左边区域
1111%沿线段方向右边区域的标识值,规则同上。
];
bs1=bs(:
)’;
iffind(bs1<1|bs1〉nbs),
error('PDE:
squareg:
InvalidBs’,’Nonexistentboundarysegmentnumber.')
end
ifnargin==1,
x=d(:
bs1);%给定一个输入变量时,输出区域边界数据的矩阵
return
end
x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
ifm==1&&n==1,
bs=bs*ones(size(s));%expandbs
elseifm~=size(s,1)||n~=size(s,2),
error('PDE:
squareg:
SizeBs’,’bsmustbescalarorofsamesizeass。
');
end
if~isempty(s),本文为互联网收集,请勿用作商业用途
%第一段边界
ii=find(bs==1);
iflength(ii)
x(ii)=interp1([d(1,1),d(2,1)],[01],s(ii));%通过参量来确定边界上的值
y(ii)=interp1([d(1,1),d(2,1)],[11],s(ii));%
end
%第二段边界
ii=find(bs==2);
iflength(ii)
x(ii)=interp1([d(1,2),d(2,2)],[11],s(ii));
y(ii)=interp1([d(1,2),d(2,2)],[10],s(ii));
end
%第三段边界
ii=find(bs==3);
iflength(ii)
x(ii)=interp1([d(1,3),d(2,3)],[10],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[00],s(ii));
end
%第四段边界
ii=find(bs==4);
iflength(ii)
x(ii)=interp1([d(1,4),d(2,4)],[00],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[01],s(ii));
end
end
编辑完该函数后,可以利用pdeplot函数来测试设置的几何区域是否满足要求。
在matlab命令窗口输入:
pdegplot(’mygeom’),axisequal
[p,e,t]=initmesh(’mygeom');
pdemesh(p,e,t),axisequal
结果如下:
mygeom的m文件放在附件里,大家可以下载。
本篇写作时参考了刘平的博客,已在前面指出,另参考了《偏微分方程的MATLAB解法》及《MATLABPDEToolboxUserGuide》,有兴趣的读者可以参阅。
MATLAB求解PDE问题(3)—-确定边界条件
前一篇给出了如何确定几何区域,本篇接着给出如何确定边界条件,在几何区域和边界条件都确定好之后,就可以利用PDEToolbox对给定的PDE进行计算了。
先回忆下前面的边界条件:
上面的四个边界条件中,前两个是Dirichlet边界,后两个是Neumann边界。
我们先了解下PDEToolbox中规定有三种边界条件,一是Dirichlet边界条件,一是广义Neumann条件,一是混合边界条件(只用于方程组)。
这和传统的定义有所不同,按照传统的观点,PDEToolbox中的广义Neumann条件应该是混合边界条件(Robin边界条件),而传统的Neumann边界条件是指边界处的法向导数为零。
我们先不考虑方程组,Dirichlet边界条件和广义Neumann条件写成如下形式:
其中
是外法线方向,
是边界法向量与x轴的夹角。
规定边界条件的m文件是pdebound函数,该函数并不提供边界条件,而是要求用户按照规定的方法给出边界条件,用户可以自己给定函数的名字。
我们计算中给出的边界条件函数是mybound。
调用格式为[q,g,h,r]=mybound(p,e,u,time)。
pdebound函数中最主要的内容是bl矩阵,bl包括了所有的边界信息。
bl的每一列都是边界条件矩阵的对应列,每一列都必须满足下面的规则:
第一行是方程的维数N,一个方程N=1,两个方程构成的方程组,N=2,…;
第二行是Dirichlet边界条件数M,Dirichlet边界条件时M=N,广义Neumann边界条件时M=0,如果是方程组,0 第三行到第3+N2—1行是表示q的字符串的长度,这个长度按与q有关的列方向的次序储存; 第3+N2行到第3+N2+N-1是表示g的字符串的长度; 第3+N2+N行到第3+N2+N+MN—1是表示h的字符串的长度,这个长度按与h有关的列方向的次序储存; 第3+N2+N+MN行到第3+N2+N+MN+M—1是表示r的字符串的长度; 接下来的行是包括MATLAB文本表达式所表示的真实边界条件函数。 文本表达式是将实际的表达式通过double函数转化为ASCII码后的表达式。 文本字符串的长度与前面四行规定的长度是一致的,两个字符串之间没有分隔符。 可以插入含有下列变量的文本表达式: 二维坐标x和y; 边界线段参数s,弧长的比例。 在边界线段的起始处s=0,向着线段的终点处增加到s=1; 外法向量分量nx,ny。 如果要用到切向量,可以用tx和ty表示,其中tx=—ny,ty=nx; 解u(除非已指定输入参数u); 时间t(除非已指定输入参数time)。 注意,在只有一个方程的时候,即N=1,如果M=0时,即是广义Neumann边界条件,此时不需要表示表示Dirichlet边界条件的h,r两行,因此表示q和g的字符串从第5行开始(前四行分别为N,M,q,g的指示行);如果M=1,此时表示h,r的字符串从第9行开始(前六行分别是N,M,q,g,h,r的指示行,第七八行表示q,g,它们均为0,ASCII码为48) 下面看下我们计算你的例子中边界mybound函数的内容: function[q,g,h,r]=mybound(p,e,u,time) %upper=double('2-2*x—x.^2');%y=1上边界条件函数 %down=double('x。 ^2’);%y=0下边界条件函数 %left=double(’y。 ^2’);%x=0左边界条件函数 %right=double(’2—2*y-y。 ^2');%x=1右边界条件函数 %upper=upper’;%转化为列向量 %down=down'; %left=left'; %right=right'; % bl=[ 1111%N表示方程维数的指示行 0011%M表示Dirichlet边界条件的指示行,M=0为广义Neumann边界条件 1111%q表示q的字符串的长度,根据实际长度确定, 101011%g表示g的字符串的长度,如1表示字符串的长度为1 484811%h表示h的字符串的长度,这是针对M=1的列,本例子是三四列 505044%r表示r的字符串的长度,规则同上。 45454848%本行以下是边界条件函数的字符串表达式,字符串的长度与上面 50504848%规定的一致。 注意,这些都是将表达式转化为ASCII码后的数 42424949%例如x的ASCII码是120,ASCII码的0表示null(空) 120121120121%需要注意的还有,M=0的列,字符表达式从g的下一行开始, 45454646%即从第5行开始。 M=1的列,从r的下一行开始,即7,8两行是 1201219494%表示q,g的字符串的长度,此时认为q,g均为0,0的ASCII 46465050%码是48,因此7,8两行均为48。 在往下,按照h,r的长度排列。 949400%由此造成的行数不匹配,可以由任意ASCII码补齐,一般以0 505000%(表示null,空)表示,不容易引起歧义。 ]; ifany(size(u)) [q,g,h,r]=pdeexpd(p,e,u,time,bl); else [q,g,h,r]=pdeexpd(p,e,time,bl); end本文为互联网收集,请勿用作商业用途 规定边界条件的矩阵bl较为复杂,详细的解释可以使大家省不少时间,但是也容易把大家整晕,我也是在这方面花费了不少功夫。 下面具体分析下bl中的两列。 第一列对应着上边界条件(Neumann边界条件): 与边界条件的标准形式比较,有,q=0,g=2—2*x-x.^2,于是第一列可以写为 [10110'0’'2-2*x—x.^3’]’ 其中1表示N=1,即一个方程,0表示Neumann边界条件;1表示q的长度;10表示g的长度,'0’表示q=0,长度为1;'2—2*x—x.^3’表示g=2—2*x—x。 ^2,转化为ASCII码后长度为10。 因此第一列最终写为 [10110485045504212045120469450]' 第三列对应着下边界条件,是Dirichlet边界条件: 与标准形式比较,有,h=1,r=x.^2,于是第三列可以写为 [111114’0’'0’'1''x。 ^2']’ 按顺序依次为,1表示N=1,一个方程;1表示Dirichlet边界条件;1表示q的长度,此时必为1;1表示g的长度,此时也必为1;1表示h的长度是1;4表示r的长度为4;'0'表示q的值为0(0的字符串长度是1);'0’表示g的值是0;'1’表示h=1(1的ASCII码是49);’x.^2'表示r=x。 ^2,转化为ASCII码后长度是4。 于是第三列可以写为 [111114484849120469450]' 此时第三列的长度小于第一列的长度,缺几个长度就可以添加几个ASCII码的0,表示空,这样第三列就是mybound中的第三列了,如下 [11111448484912046945000]’ 最后附上mybound的m文件,有需要的可以下载
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 求解 PDE 问题
![提示](https://static.bdocx.com/images/bang_tan.gif)