偏微分方程Matlab数值解法（补充4）

偏微分方程Matlab数值解法（补充4） 偏微分方程Matlab数值解法(补充4) Matlab可以求解一般的偏微分方程，也可以利用偏微分方程工具箱中给出的函数求解一些偏微分方程。 1 偏微分方程组求解 Matlab语言提供了pdepe()函数，可以直接求解偏微分方程 ,,,,,uuuu,mm (4.1) (,,,)[(,,,)](,,,)Cxtuxxfxtusxtu,， ,,,,,xtxxx 这样，偏微分方程可以编写为以下函数的描述，其入口为 [,,](,,,)cfspdefunxtuu,x cfs,,其中:pdefun为函数名。由给定输入变量可计算出这三个函数。 边界条件可以用下面的函数描述 ,u (4.2) pxtuqxtufxtu(,,)(,,).*(,,,)0，, ,x 这样的边值函数可以写为Matlab函数 [,,,](,,,)pqpqpdebcxtuu,aabbx 初始条件数学描述为，编写一个简单的函数即可 uxtu(,),00 updeicx,() 0 tx还可以选择和的向量，再加上描述这些函数，就可以用pdepe()函数求解次偏微分方程，需要用如下 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 求解 solpdepempdefunpdeicpdebcxt,(,@,@,@,,) 【例1】 试求下列偏微分方程 2,,,uu11,,,0.024()Fuu12,2,,,tx ,2,,uu,22,，,0.17()Fuu122,,,tx, 5.7311.46xx,ux(,0)1,其中:Fxee(),,，且满足初始条件，及边界ux(,1)0,12条件 ,,uu12(0,)0,(0,)0,(1,)1,(1,)0tututt,,,, 21,,xx 1 解:对照给出的偏微分方程和(4.1)，可将原方程改写为 uuxFuu0.024/(),,,,1，，，，，，，，,,1112 .*,，,,,,,,,,10.17/()uuxFuu,,,,,tx，，2212，，，，，， ，且 可知m,0 0.024/(),,,,uxFuu1，，，，，，112 cfs,,,,,,,,,,,10.17/(),,,uxFuu，，212，，，， 编写下面的Matlab函数 function [c,f,s]=c7mpde(x,t,u,du) c=[1;1];y=u(1)-u(2);F=exp(5.73*y)-exp(-11.46*y);s=F*[-1;1]; f=[0.024*du(1);0.17*du(2)]; 套用(4.2)中的边界条件，可以写出如下的边值方程 0u,11000，，，，，，，，，，，，1，右边界 左边界，,.*f，,.*f,,,,,,,,,,,,u001002，，，，，，，，，，，， 编写下面的Matlab函数 function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t) pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0,1]; 另外，描述初值的函数 function u0=c7mpic(x) u0=[1;0]; tx有了这三个函数，选定和向量，则可以由下面的程序直接求此微分方程， uu得出解和。 12 x=0:0.05:1;t=0:0.05:2;m=0; >> sol=pdepe(m,@c7mpde,@c7mpic,@c7mpbc,x,t); >> surf(x,t,sol(:,:,1)) >> figure;surf(x,t,sol(:,:,2)) 2 二阶偏微分方程的属性描述 2.1 椭圆型偏微分方程 椭圆型偏微分方程的一般形式为 ,,，,divcuaufxt()(,) u其中:若，为的梯度，则其定义为 uuxxxtuxt,,(,,,,)(,),u12n ，，,,,,,uu,,, ,,,,,xxx12n，， divv()散度的定义为 2 ,,,,, divvv(),，，，,,,,,xxx12n,, divcu(),这样，可以更明确地 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示为 ，，,,,,,,,,,,,,uuu()divcuccc,,，，， ,,,,,,,,,,,,,,xxxxxx1122nn,,,,,,，， 若为常数，则进一步化简为 c 222,,,,, divcucucu(),,，，，,,,,222,,,xxxn,,12 其中，又称为Laplace算子。这样椭圆型偏微分方程可以简单地写为 , 222,,,,, ,，，，，,cuaufxt(,),,222,,,xxxn,,12 2.2 抛物型偏微分方程 抛物型偏微分方程的一般形式为 ,u ddivcuaufxt,,，,()(,) ,t 根据上面叙述，若c为常数，则该方程可以更简单地写为 222,,,,,,u(,) dcuaufxt,，，，，,,,222,,,,txxxn12,, 2.3双曲型偏微分方程 双曲型偏微分方程的一般形式为 2,uddivcuaufxt,,，,()(,) 2,t c若为常数，则可以将该方程简化为 2222,,,,,,u(,) dcuaufxt,，，，，,,,2222,,,,txxxn12,, tu三类方程的直接的区别在于对的导数的阶次。 t若对没有求导，可以理解为其值为常数，故称为椭圆型的。 tuux若取对时间的一阶导数，则与对的二阶导数直接构成了抛物线关系， 故称为抛物型偏微分方程。 tu若取对时间的二阶导数，称其为双曲型偏微分方程。 3 2.4 特征值型偏微分方程 特征值型偏微分方程为 ,,，,divcuaudu(), 对常数该方程可以化简为 c 222,,,,, cuaufxt,(,),，，，，,,,222xxx,,,n12,, 该方程是椭圆型偏微分方程的特例。 3 偏微分方程求解界面应用举例 3.1 偏微分方程求解程序概述(已讲) 3.2 偏微分方程求解区域绘制(已讲) 3.3 偏微分方程边界条件描述 求解边界在偏微分方程界面下用。在偏微分方程工具箱支持包括,,Dirichlet条件和Neumann条件。 (1) Dirichlet条件。 一般描述为 uu,,,,,,hxtuurxtu,,,,,,, ,,,,,,xx,,,,,, r其中:表示求解区域的边界。假设在边界上满足该方程，则只需要给出和,,h xuux,/,,函数即可，这两个参数可以为常数，也可以为的函数，甚至可以是的 函数，为了方便起见，一般可以令。 h,1 (2) Neumann条件。 其扩展形式为 ,，，()cuqug,，, ,,,,n,，， x其中:为向量法向量的偏导数。 ,,un/ hr,,1,0选择Boundary?Specify Boundary Condition菜单，设置即可。 3.4 偏微分方程求解举例 【10.2】试求双曲型偏微分方程 4 222duuu,,,,，,210u 222dtxy,, 解:(双曲型方程的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 形式为 2222,,,,,,u (,)dcuaufxt,，，，，,,,2222,,,,txxxn12,, cafd,,,,1,2,10,1因此，。 区域绘制:(R1+E1+E2)-E3 网格绘制(网格加密) Boundary?Specify Boundary Condition菜单，设置hr,,1,0。 PDE?PDE Specification? Parabolic选项，输入 cafd,,,,1,2,10,1 也可修改边界条件为(边界上所有的值为5)。 r,5 ?Plot Solution Plot Paramters 3.5 时变动画显示 偏微分方程默认的时间向量，以所给的是最终时刻时的解。t,0:10t,10从双曲型偏微分方程看，方程的解应该是时间t的函数，所以应该用动画来显示出来。 (以例2为例) 可由Solve?Parameters菜单引出的对话框设置时间向量，在time栏内填入，这样求解即为这段时间的解。 0:0.1:4 选择其中的动画选项:Animation 单击Options设置动画的播放速度，默认为6帧/秒，这样就可获得动画。 可以用Plot?Export Movie菜单将动画输出到Matlab工作空间，例如存成变量M，则可以用movie(M)在Matlab图形窗口中播放得出的动画，也可以用movie2avi(M,'myavi.avi')命令将动画存成myavi.avi文件，以备以后播放。 3.6 函数参数的偏微分方程求解 cadf,,,前面介绍的偏微分方程均为常数，而实际遇到的方程中，经常遇到这些量为函数的情况。 ,,uu用分别表示偏导数。 uxuy,, ,,xy 【例2】 假设偏微分方程为 5 22,,122(),，xy divuxyue,,，，,(),,2u1||，,,, 其中边界为0，求解该偏微分方程。 解 观察该方程发现，它是椭圆型偏微分方程。 ,,，,divcuaufxt()(,)(标准型为) 22122(),，xycaxyfe,,，,,,其中 22,,,,uu,,1，，,,,,,,xy,,,, PDE?PDE Specification?Elliptic选项， 1./(1.^2.^2)sqrtuxuy，，在栏输入; c xy.^2.^2)，在a栏输入; exp(.^2.^2),,xyf在栏输入; 打开Solve?parameters对话框，从中选定Use nonlinear solver属性(注意，该属性只适用于椭圆型偏微分方程求解)。 4 偏微分方程应用实例的Matlab实现 细棒的热流问题:设有一细棒由各向均匀材料制成，长度为a。 utx(,) 设表示细棒上点u(t,x) cold hot x=0 x x=a tx处在时刻的温度，有物理 utx(,)定律知满足热传导方 程: 2,,uu,, 2,,tx ,其中，为比热系数，取决于细棒的材料。为了确保唯一性，设细棒的两个端点 绝热，即在端点没有热量流失。从而有边界条件 ututa(,0)0,(,)0,,xx ,,0.005设棒长，比热系数，初始时刻细棒的温度分布为 a,2 0uxxC(0,)3010[1cos()],，,, o温度分布为一正弦曲线，在细棒的中间，温度达到50C，在细棒的两端，温 6 o度低于30C。取时间取决，时间步长，空间步长，t,[0,120],,t2.4,,x0.05 运用6点隐式格式求解数值解。设在第k个时刻t，棒上各点的温度为utx(,)，kjk ,,t记，用6点隐式格式离散偏微分方程后，得线性代数方程组 ,,2,x 22，,,,,,,，，，， ,,,,,，,,2(1)2(1),,,,,,,,,,kk，1,,,,uu,,,,,,，,,2(1)2(1),,,,,,,,,, ,,,,,，,22，，，，,,,, k其中u表示向量的近似值。这是一个三对角方utxutxutx(,),(,),,(,)kknn1211,,程组，可以用追赶法求解。程序如下 function [t,x,U]=heatrod(m,n,beta) dt=120/m; dx=2/n; gamma=beta*dt/(dx*dx); t=zeros(m+1,1); %t网格值 x=zeros(n-1,1); %x网格值 U=zeros(m+1,n-1); %解函数 b=zeros(n-1,1); li=zeros(n-1,1); ui=zeros(n-1,1); t=[0:dt:120]; x=[dx:dx:(n-1)*dx]'; f=inline('30+10*(1-cos(pi*x))','x'); %初值函数 for j=1:n-1 U(1,j)=f(x(j)); end g0=2*(1+gamma);g1=2*(1-gamma); for k=1:n-2 A(k,k+1)=-gamma; A(k+1,k)=-gamma; A(k,k)=g0; end A(1,1)=A(1,1)-gamma; A(n-1,n-1)=g0-gamma; %分解矩阵A ui=A(1,1); 7 for i=2:n-1 li(i)=A(i,i-1)/ui(i-1); ui(i)=A(i,i)-A(i-1,i)*li(i); end for k=1:m b(1)=(g1+gamma)*U(k,1)+gamma*U(k,2); for j=2:n-2 b(j)=gamma*U(k,j-1)+g1*U(k,j)+gamma*U(k,j+1); end b(n-1)=gamma*U(k,n-2)+(g1+gamma)*U(k,n-1); %用追赶法求解线性方程组AU=b y(1)=b(1); for i=2:n-1 y(i)=b(i)-li(i)*y(i-1); end; U(k+1,n-1)=y(n-1)/ui(n-1); for i=n-2:(-1):1 U(k+1,i)=(y(i)+gamma*U(k+1,i+1))/ui(i); end uu=zeros(m+1,n+1); for i=1:m uu(i,1)=U(i,1); uu(i,n+1)=U(i,n-1); for j=1:n-1 uu(i,j+1)=U(i,j); end end end 在Matlab命令窗输入[t,x,U]=heatrod(50,40,0.005)，回车得到运算结果。 5 偏微分方程Matlab常用函数举例 【例4】求L形状区域上的偏微分方程，在边界上，并画出函,,,u1u,0数的解。 解 在Matlab命令窗输入 [p,e,t]=initmesh('lshapeg','Hmax',0.2); >> [p,e,t]=refinemesh('lshapeg',p,e,t); >> u=assempde('lshapeb',p,e,t,1,0,1); >> pdesurf(p,t,u) Assempde是PDE工具箱中最基本的函数，它用有限元的方法生成刚度矩阵 和右端向量，该函数主要求解下列标量形式的偏微分方程 ,,,,，,()cuauf 8 及系统微分方程 ,,,,,，,()cuauf 】求解小方程 【例5 ,,,uu, 在L形区域上的小于100的特征值及相应的特征模态，并演示第16个特征值模 态。 在命令窗输入 解 [p,e,t]=initmesh('lshapeg'); >> [p,e,t]=refinemesh('lshapeg',p,e,t); >> [v,l]=pdeeig('lshapeb',p,e,t,1,0,1,[-Inf 100]); >>l(1) >>pdesurf(p,t,v(:,1)) >> figure >> membrane(1,20,9,9); >> figure >> l(16) >> pdesurf(p,t,v(:,16)) Pdeeig函数解特征值问题，是生成用有限元法求解的特征方程: ,,,,，,()cuaudu, 或 ,,,,,，,()cuaudu, 【例6】在以L形定义的几何层中画方程解的等高线，Dirichlet,,,u1 u,,,0()用边界条件。 解 在命令窗输入 [p,e,t]=initmesh('lshapeg'); >> [p,e,t]=refinemesh('lshapeg',p,e,t); >> u=assempde('lshapeb',p,e,t,1,0,1); >> pdecont(p,t,u) 【例7】解Possion方程，几何条件由L形平面定义。Dirichlet,,,u1 u,,,0()用边界条件，且画出结果。 解 在命令窗输入 [p,e,t]=initmesh('lshapeg'); [p,e,t]=refinemesh('lshapeg',p,e,t); >> u=assempde('lshapeb',p,e,t,1,0,1); >> pdemesh(p,e,t,u) 9 得到Possion方程的解。 【例8】精细化L型膜若干次 解 在命令窗输入 [p,e,t]=initmesh('lshapeg','Hmax',inf); subplot(2,2,1),pdemesh(p,e,t); [p,e,t]=refinemesh('lshapeg',p,e,t); subplot(2,2,2),pdemesh(p,e,t) [p,e,t]=refinemesh('lshapeg',p,e,t); subplot(2,2,3),pdemesh(p,e,t); >> [p,e,t]=refinemesh('lshapeg',p,e,t); >> subplot(2,2,4),pdemesh(p,e,t) >> subplot 【例9】求热传导方程 ,u ,,u ,t 22,,,11yu(0)1,求解区域为正方形，，初始条件:当，;xy，,1,,,11x u(0)0,其它情况。 解 在命令窗输入 [p,e,t]=initmesh('squareg'); >> [p,e,t]=refinemesh('squareg',p,e,t); >> u0=zeros(size(p,2),1); >> ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4); >> u0(ix)=ones(size(ix)); >> tlist=linspace(0,0.1,20); >> u1=parabolic(u0,tlist,'squareb1',p,e,t,1,0,1,1); 执行上述程序段，得到如下结果: 经过96步计算;194次函数赋值;1次求偏导;20次LU分解运算;193组 偏微分方程的数值解。 10