偏微分方程求解

1、关于偏微分方程求解现在学习的是第一页，共46页基础知识预习 微分方程的求解包含 ：常微分方程的求解（上节课已经讲过）这里不再赘述。 ：偏微分方程的求解（本次教学内容）现在学习的是第二页，共46页偏微分方程概念 偏微分方程（Partial Differential Equation，简称PDE）指含有未知函数及其偏导数的方程。描述自变量、未知函数及其偏导数之间的关系。 偏微分方程分为线性偏微分方程式与非线性偏微分方程式，常常有几个解而且涉及额外的边界条件。常微分方程：在微分方程中，若自变量的个数只有一个的微分方程。偏微分方程：自变量的个数有两个或两个以上的微分方程。现在学习的是第三页，共46页求

2、解偏微分方程的方法 求解偏微分方程的数值方法： 1. 有限元法（Finite Element Method, FEM）- hp-FEM 2. 有限体积法（Finite Volume Method, FVM） 3. 有限差分法（Finite Difference Method, FDM）。其它：广义有限元法（Generalized Finite Element Method, FFEM）、扩展有限元法（eXtended Finite Element Method, XFEM）、无网格有限元法（Meshfree Finite Element Method）、离散迦辽金有限元法（Discontinu

3、ous Galerkin Finite Element Method, DGFEM）等。现在学习的是第四页，共46页MATLAB解偏微分方程 MATLAB提供了两种方法解决PDE 问题：pdepe()函数，它可以求解一般的PDEs，具有较大的通用性，但只支持命令行形式调用。 PDE 工具箱，可以求解特殊PDE 问题，PDEtool 有较大的局限性，比如只能求解二阶PDE 问题，并且不能解决偏微分方程组，但是它提供了GUI界面，从繁杂的编程中解脱出来了，同时还可以通过File-Save As直接生成M代码 使用pdeval()直接计算某个点的函数值？现在学习的是第五页，共46页一般偏微分方程组(

4、PDEs)的MATLAB求解 直接求解一般偏微分方程(组)，它的调用格式为sol=pdepe(m,pdefun,pdeic,pdebc,x,t)(1) ) u,t,s(x,),()u,( cxuxuutxfxxxtuxutxmm，问题描述函数初值条件边界条件输出参数自变量参数现在学习的是第六页，共46页【输入参数】（1） pdefun：是PDE 的问题描述函数，它必须换成下面的标准形式 PDE 就可以编写下面的入口函数 c,f,s=pdefun(x,t,u,du) m,x,t就是对应于(式1)中相关参数和自变量，du是u的一阶导数，由给定的输入变量即可表示出出c,f,s这三个函数(1) ) u

5、,t,s(x,),()u,( cxuxuutxfxxxtuxutxmm，现在学习的是第七页，共46页【输入参数】（2） pdeic：是PDE 的初值条件，必须化为下面的形式 我们使用下面的简单的函数来描述为u0=pdeic(x)00),(uutx现在学习的是第八页，共46页【输入参数】（3） pdebc：是PDE的边界条件描述函数，必须先化为下面的形式 于是边值条件可以编写下面函数描述为pa,qa,pb,qb=pdebc(x,t,u,du)其中a 表示下边界，b 表示下边界0 ) xu u,t,(x, f *u).t,q(x, u) t,p(x,现在学习的是第九页，共46页【输入参数】（4）

6、m：就是对应于(式1)中相关参数 x,t：就是对应于(式1)中自变量(1) ) u,t,s(x,),()u,( cxuxuutxfxxxtuxutxmm，现在学习的是第十页，共46页【输出参数】 sol：是一个三维数组，sol(:,:,i)表示ui的解，换句话说uk对应x(i)和t(j)时的解为sol(i,j,k)现在学习的是第十一页，共46页实例讲解（题目） 例：现在学习的是第十二页，共46页初值条件边界条件现在学习的是第十三页，共46页实例讲解（解法） 【解】第一步根据（1）对照给出的偏微分方程，则原方程可以改写为现在学习的是第十四页，共46页输入参数（1）目标PDE函数 % 目标PDE函

7、数 function c,f,s=pdefun (x,t,u,du) c=1;1; f=0.024*du(1);0.17*du(2); temp=u(1)-u(2); s=-1;1.*(exp(5.73*temp)-exp(-11.46*temp);现在学习的是第十五页，共46页输入参数（2）初值条件 初值条件改写为 % 初值条件函数 function u0=pdeic(x) u0=1;0;现在学习的是第十六页，共46页输入参数（3）边界条件 边界条件改写为 % 边界条件函数 function pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t) %a表示左边界，b表示右边界 pa

8、=0;ua(2);qa=1;0; pb=ub(1)-1;0;qb=0;1;现在学习的是第十七页，共46页（4）主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t);figure(numbertitle,off,name,PDE Demoby Matlabsky)%创建个窗口，窗口名字是name后边的名字NumberTitle,off是关掉默认显示名字。subplot(211)surf(x,t,sol(:,:,1)%sol(:,:,i)表示ui的解title(The Solution of u_1)xlabel(

9、X)ylabel(T)zlabel(U)subplot(212)surf(x,t,sol(:,:,2)%sol(:,:,i)表示ui的解title(The Solution of u_2)xlabel(X)ylabel(T)zlabel(U)现在学习的是第十八页，共46页现在学习的是第十九页，共46页PDEtool求解特殊PDE问题 MATLAB的偏微分工具箱(PDE toolbox)可以比较规范的求解各种常见的二阶偏微分方程(特殊二阶的PDE)现在学习的是第二十页，共46页典型偏微分方程的描述 现在学习的是第二十一页，共46页 现在学习的是第二十二页，共46页 （3）双曲线型偏微分方程的一般

10、形式 现在学习的是第二十三页，共46页 （4）特征值型偏微分方程的一般形式，注 意它是（1）的变形，不能算独立的一类 现在学习的是第二十四页，共46页MATLAB 采用有限元的方法求解各种PDE MATLAB 为我们提供一个pdetool (在command window 中键输pdetool打开)的交互界面，可以求解二元偏微分u(x1,x2)(注意只能求解二元)。方程的参数由a、c、d和f确定，求解域由图形确定，求解域确定好后，需要对求解域进行栅格化(这个是自动)。现在学习的是第二十五页，共46页 现在学习的是第二十六页，共46页 现在学习的是第二十七页，共46页偏微分方程边界条件的描述 D

11、irichlet(狄利克莱)条件 Neumann(纽曼)条件现在学习的是第二十八页，共46页 现在学习的是第二十九页，共46页 现在学习的是第三十页，共46页求解实例 现在学习的是第三十一页，共46页 【解】由给定的PDE，可以得出d=1,c=1,a=2,f=10现在学习的是第三十二页，共46页step1:点击工具栏的【PDE】按钮，如下输入PDE的参数，注意选择Hyperbolic 现在学习的是第三十三页，共46页step2:绘制求解域对坐标轴的操作可以在【Options】主菜单中操作，包括设置网格、坐标系范围等 (1)【Options】-Axis Limits设置如下现在学习的是第三十四页

12、，共46页 现在学习的是第三十五页，共46页 (2)点击工具栏上的第三个按钮【绘制椭圆】，任意绘制一个椭圆，双击椭圆，设置如下 现在学习的是第三十六页，共46页重复上面的操作，参数如下 现在学习的是第三十七页，共46页得到 现在学习的是第三十八页，共46页 (3)在set formula 中如下输入，“+”表示求并集，“-”表示求差集，注意没有直接求交接的操作符现在学习的是第三十九页，共46页step3:边界条件和初值条件 初值条件可以通过【Solve】-【Parameters】设置 边值条件设置如下 (1)点击工具栏的第6 个按钮【区域边界】，显示如下现在学习的是第四十页，共46页 (2)【

13、Boundary】-【Remove All Subdomain Borders】移除所有子域的边界，将得到所有子域合并成一个求解域 (3) 【Boundary】-【Secify Boundary Conditons】设置边界如下，注意我们这里只有Dirichlet条件现在学习的是第四十一页，共46页step4:生成使用有限元方法求解方程所需的栅格 点击工具栏的第8/9 个按钮，对求解域生成栅格，多次点击可以在原来基础上继续细化栅格，直到自己觉得满意 为止，当然可以通过【Mesh】主菜单进行精确控制现在学习的是第四十二页，共46页 step5:求解方程点解工具栏的第10 个按钮“=”【求解方程】 step6:求解结果绘图 点击第11 个按钮【绘制图形】，里面的选项很丰富，可以绘制等高线等好多，甚至播放动画，具体大家可以自己慢慢摸索现在学习的是第四十三页，共46页动画播放设置：(1)【Solve】-【Parameters】设置合适的时间向量Time(2)【Plot】-【Parameters】选中【Animation】，点击后面的【Options】，设置