matlab pde瞄点,[转载]Matlab 中PDE一些函数的用法

Matlab 中PDE一些函数的用法

Matlab工具箱中任何功能都能用命令行的形式解决，对于一些复杂的问题，gui形式是解决不了的，这时就必须用命令行形式自己编程来解决问题。最近查看了一些pde工具箱中的函数，对两个函数的用法做一些说明。

PDE解决问题流程与一般的fam方法一样，先画出几何模型，再画网格，加边界、载荷，求解。MATLAB中的画几何命令有两种方式，其中一种方式是自己编写m文件，MATLAB中的命令是pdegeom，实际上这不是一个命令，而是要用户自己编写的一个几何m文件。在网上找了半天，也没见有人介绍过这个函数的用法，网上大多数只是列出pde中的函数而已，感叹一下….

这个函数名可以任意，假设为tt，参数有两个(bs,s)如果输入参数为0，1时，函数也必须要能有返回值。函数的意义为：

参数为0个，即没有参数时，返回边界总共段数。

参数为1个是，即只有bs，返回边界指定段的参变量范围。

参数为2个时，返回每段边界长度为s时的坐标。

函数参数意义bs表示几何边界段数，如矩形，边界为四段，三角开边界肯定为三段….s为第bs段时的长度，bs与s可以为向量，但是要一一对应，即bs为几个值，s也得为几个值。

下面具体来讲下函数如何实现当参数为2时的返回值。其实函数中实现的内容是指用一个参变量来控制指定段的边界坐标，也就是说对任意一段边界，上面的坐标设为参变量的函数，假设参变量为t,则对任一边界，坐标为(x(t),y(t))，参变量的范围由用户自己定义。由以上可以看出，参变量其实有无穷多种取法，看自己怎么设定而已。

例如：一个矩形，四条边，返回每一段对应长度的坐标，m文件为：

function [x,y]=tt(bs,s);

ux=1;

uy=1;

nbs=4;  //表示有四段

if nargin==0

x=nbs;

return

end

dl=[0,0,0,0; //参变量t范围，第一行为起始值，第二行为结束值，这两个值理论

3,1,2,3;

//可以任意设置，只要后面能对应求出边界上坐标(x(t),y(t))就行

1,1,1,1;

//第三行指面积在边界左侧

0,0,0,0,];

//第四行指面积在边界右侧

if nargin==1

x=dl(:,bs);

return;

end

x=zeros(size(s));

y=zeros(size(s));

[m,n]=size(bs);

if m==1 & n==1

bs=bs*ones(size(s));

end

ii=find(bs==1);

if ~isempty(ii)

x(ii)=interp1([dl(1,1),dl(2,1)],[0 3],s(ii));

y(ii)=interp1([dl(1,1),dl(2,1)],[0,0],s(ii));

//其实这一段是一条水平直线，

//x=(3-0)/(dl(2,1)-dl(1,1))*s(1)=s(1)

//y=(0-0)/(dl(2,1)-dl(1,1))*s(1)=0,

//如果dl(2,1)由3变成1，那么

//x(ii)=interp1([dl(1,1),dl(2,1)],[0 3],s(ii));结果也是一样的。

end

ii=find(bs==2);

if length(ii)

x(ii)=interp1([dl(1,2),dl(2,2)],[3*ux,3*ux],s(ii));

y(ii)=interp1([dl(1,2),dl(2,2)],[0,1],s(ii));

end

ii=find(bs==3);

if length(ii)

x(ii)=interp1([dl(1,3),dl(2,3)],[3*ux,-2*ux+2],s(ii));

y(ii)=interp1([dl(1,3),dl(2,3)],[uy,uy],s(ii));

end

ii=find(bs==4);

if ~isempty(ii)

x(ii)=interp1([dl(1,4),dl(2,4)],[0 0],s(ii));

y(ii)=interp1([dl(1,4),dl(2,4)],[uy,0],s(ii));

end

end

图形如下所示：

这个函数编制的关键是，函数内边界上的坐标((x(t),y(t))是用参变量表示的，返回值是求得边界任意长度时的坐标(x(t),y(t))值，明白这一点其它就很容易了。

另外，值得注意的是，函数所绘制图形必须是封闭的，否则之后划分网格会出错的。

MATLAB帮助文档中有一个曲线的例子，可以对照着看，之所以不用帮助文档上的例子是因为我在网上找了大半天，绝大多数都是对函数的罗列，还有一些就是对帮助文档的翻译，不知道是觉得太简单不想多说还是本来就没弄明白。

说了这么多也不知道讲清楚了没有，希望能对大家有帮助。

% pdegplot('tt'), axis equal

% [p,e,t]=initmesh('tt');

% pdemesh(p,e,t), axis equal