偏微分方程的MATLAB解法 數學物理方法

偏微分方程的matlab解法

主要講述如何用MATLAB實現對偏微分方程的仿真求解．MATLAB的偏微分方程工具箱（PDE Toolbox）的出現，為偏微分方程的求解以及定性研究提供了捷徑．主要步驟為：

1、設置PDE的定解問題．即設置二維定解區域、邊界條件以及方程的形式；和系數
2、用有限元法（FEM）求解PDE．即網格的生成、方程的離散以及求出數值解；
3、解的可視化．

PDEToolbox注意事項
只能解決二維模型，一維的擴成二維，三維的縮成二維，時間維不計算在內
公式類型，只能解決部分偏微分方程，由公式類型決定
邊界條件兩種，Dirichlet和Neumann
初始條件




第一題

g='circleg';
b='circleb1';
c=1;
a=0;
f=1;
[p,e,t]=initmesh(g,'hmax',1);
figure;
pdemesh(p,e,t); axis equal
er = Inf;
while er > 0.001
    [p,e,t]=refinemesh(g,p,e,t);
    u=assempde(b,p,e,t,c,a,f);
    exact=(1-p(1,:).^2-p(2,:).^2)'/4;
    er=norm(u-exact,'inf');
    fprintf('Error: %e. Number of nodes: %d\n',er,size(p,2));
end
figure;
pdemesh(p,e,t); axis equal
figure;
pdesurf(p,t,u-exact);
figure;
pdesurf(p,t,u);


第二題

a=0;
b=1;
c=0;
d=1;
r='squareg';
z='squareb3';
[p,e,t]=initmesh('squareg');
figure;
pdemesh(p,e,t); axis equal
x=p(1,:)';
y=p(2,:)';
u0=atan(cos(pi/2*x));
ut0=3*sin(pi*x).*exp(sin(pi/2*y));
n=31;
tlist=linspace(0,5,n);
uu=hyperbolic(u0,ut0,tlist,z,p,e,t,b,a,c,d);
figure; set(gcf,'renderer','zbuffer');
delta=-1:0.1:1;
[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta);
gp=[tn;a2;a3];
newplot;
umax=max(max(uu));
umin=min(min(uu));
for i=1:n
    pdeplot(p,e,t,'xydata',uu(:,i),'zdata',uu(:,i),'zstyle','continuous',...
                  'mesh','off','xygrid','on','gridparam',gp,'colorbar','off');
    axis([-1 1 -1 1 umin umax]); caxis([umin umax]);
    M(i)=getframe;
end
movie(M,1);


第三題

g='squareg';
b='squareb1';
c=1;
a=0;
f=1;
d=1;
[p,e,t]=initmesh(g);
figure;
pdemesh(p,e,t); axis equal
u0=zeros(size(p,2),1);
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);
u0(ix)=ones(size(ix));
nframes=20;
tlist=linspace(0,0.1,nframes);
u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);
figure; set(gcf,'Renderer','zbuffer');
newplot;
colormap(cool);
x=linspace(-1,1,31);y=x;
[~,tn,a2,a3]=tri2grid(p,t,u0,x,y);
umax=max(max(u1));
umin=min(min(u1));
for j=1:nframes,
    u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));
    surf(x,y,u);caxis([umin umax]);
    axis([-1 1 -1 1 0 1]);
    shading interp;
    Mv(j) = getframe;
end

全部資料51hei下載地址：
數學物理方法.rar (1.97 MB, 下載次數: 6)

2019-1-22 20:17 上傳

點擊文件名下載附件
下載積分: 黑幣 -5