1、3关于偏微分方程的 MATLAB 函数与例题pdepe解决具有一维时间 t，和空间变量 x 的抛物线，或椭圆偏微分方程的初值与边界值的系统。函数形式：solpdepe（m，pdefun，icfun，bcfun，xmesh，tspan）solpdepe（m，pdefun，icfun，bcfun，xmesh，tspan，options）solpdepe（m，pdefun，icfun，bcfun，xmesh，tspan，options，1，2）其中，m，表示解决系统的类型，m 可以取（表示板问题），（表示圆柱问题），（表示球 形问题）。pdefun，定义所要解决的偏微分方程。形式是：c，f，s=pd
2、efun（x，t，u，dudx）icfun，定义系统的初始条件。形式：u=icfun（x）bcfun，定义系统的边界条件。形式;pl,ql,pr,qr=bcfun(xl,ul,xr,ur,t)xmesh，tspan，对于任何一个偏微分方程，我们都可以将其变为如下形式：该偏微分方程的条件为：t 0 t tf 和 a x b。其中a,b的间隔必须是有界的，m 可以取0，1，2如果 m 0，a 必须 0。在上面的方程中，f（x,t,u, ）是一个通量项，s（x,t,u, ）是源项。xu/ xu/对于 t=0 时，该偏微分方程满足的初始条件的形式是：u（x，t 0）=u 0（x）对于所有的 t，和当
3、x=a 或者 x=b 时，该偏微分方程满足的边界条件形式是：p（x，t，u）+q（x，t）f（x,t,u, ）=0xu/q 的值或者是 0，或者不是 0，注意边界条件用通量 f 表示胜于用 表示，p 是 u 来决定的.xu/例 1一块厚度为 13mm 的钢板，其初始温度各处均匀为 35C。现突然将其置于某物质中，使两端面温度骤然升至 146C，比维持此温度不变，是求算钢板中心面温度上升至 145C 时所经历的时间。已知钢板的导温系数为 a=2.60 10-4m2/h解：该题为平板的两个端面维持恒温的不稳态传热的问题，其热传导方程及相应的初始条件和边界条 件如下：at2xt初始条件：， ， 边界
4、条件：，x l， ，x， xtMATLAB 解题方法如下：1 我们将上述热传导方程改写成 MATLAB 函数 pdepe 解偏微分方程时所需的标准形式，=x0 +0a1txt从上式我们可以看出：m=0，c = ,f = ,s =0;tux,tu,xxtu,2 在 MATLAB 中编写 pdefun 子程序，如下：function c，f，s=pdex1pdew（x，t，u，dudx）c=1/2.6;f=dudx;s=03 在 MATLAB 中编写问题的初始条件 icfun 子程序，如下：function u0=pdex1icw(x)u0=35;4 在 MATLAB 中编写问题的边界条件 bcf
5、un 子程序，如下：为了写出边界条件，必须把已知边界条件变为如下形式：当 x=l 时，t(l)-t0+0 =0xt当 x=-l 时，t(-l)-t0+) 0 =0t所以边界条件的子程序如下：function pl,ql,pr,qr=pdex1bcw(xl,ul,xr,ur,t)pl=ul-146;ql=0;pr=ur-146;qr=0;5.在 MATLAB 命令窗口键入如下命令：x=0：1：13；t=0：1：50；m=0；sol=pdepe（m,pdex1pdew,pdex1icw,pdex1bcw,x,t）6.运行结果：u=sol（：，：，1）；surf（x，t，u）图形如下：plot（x，
6、u（end，：）划出最终曲线图。由上图我们可以看出当钢板中心面升至 145 度所经历的时间为 31 秒。例 2一厚度为 46.2mm，温度为 278K 的方块奶油由冷藏室移至 298K 的环境中，奶油盛于容器中，除顶面与环境直接接触外，各侧面和底面军包在容器之内。设容器为绝热体。试求算 5 小时后奶油顶面、中心面及底面处的温度。已知奶油的导热系数 k、比热 c、密度 分别为 0.197W/m K、2300J/kg K、998kg/m 3，奶油表面与环境之间的膜系数 h 为 8.52W/m2 K。解： 本题属于平板以为不稳态导热问题，在此情况下的热传导方程式为：at2xt初始条件： =0,t=t
7、0边界条件：x=l，-k =hts( )-tbxtx=-l,k =hts( )-tb边界条件中的 ts( )为任一瞬时平板表面的温度，此温度随时间而变；t b为流体介质的主体温度，假定为恒定值。MATLAB 解题方法如下：5 我们将上述热传导方程改写成 MATLAB 函数 pdepe 解偏微分方程时所需的标准形式，=x0 +0a1txt从上式我们可以看出：m=0，c = ,f = ,s =0;xtu,a1xtu, xtu,6 在 MATLAB 中编写 pdefun 子程序，如下：function c，f，s=pde2pdew（x，t，u，dudx）c=1/8.58e-8;f=dudx;s=07
8、 在 MATLAB 中编写问题的初始条件 icfun 子程序，如下：function function u0=pde2icw(x)u0=278;8 在 MATLAB 中编写问题的边界条件 bcfun 子程序，如下：为了写出边界条件，必须把已知边界条件变为如下形式：当 x=l 时，- hts( )-tb-k =0xt当 x=-l 时，-hts( )-tb+k =0t所以边界条件的子程序如下：function pl，ql，pr，qr=pde2bcw(xl,ul,xr,ur)h=8.52;k=0.197;ub=298;pl=-h*(ul-ub);ql=k;pr=h*(ur-ub);qr=k;5.在
9、MATLAB 命令窗口键入如下命令：x=-0.0462：0.0002：0.0462；t=0：100：5 3600；m=0；sol=pdepe（m,pde2pdew,pde2icw,pde2bcw,x,t）6.运行结果：u=sol（：，：，1）；surf（x，t，u）图形如下：plot（x，u（end，：）划出最终曲线图。对于奶油顶面的温度：uout,DoutDx=pdeval(m,x,u(:,:),0)uout =291.3900DoutDx =81.0628所以奶油顶面的温度为：291.39K对于奶油中面的温度：uout,DoutDx=pdeval(m,x,u(:,:),-0.0231)uo
10、ut =289.1503DoutDx =126.4672所以奶油中面的温度为：289.15K对于奶油底面的温度：uout,DoutDx=pdeval(m,x,u(:,:),0.0462)uout =288.5792DoutDx =161.4001所以奶油底面的温度为：288.58K4关于常微分方程的 MATLAB 函数与例题1Ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb解决初始条件下的常微分方程。函数形式：T，Y=solver（odefun，tspan，y0）T，Y=solver（odefun，tspan，y0，options）其中，odefu
11、n 定义所要解决的常微分方程，必须以如下形式定义：y =f(t,y)或者 M(t,y)=f(t,y)tspan 该常微分方程的时间变量。y0，为该常微分方程的初始条件。普兰德边界层方程的精确解。根据平板边界层的特点，已经证明在 x 方向上的压力梯度为零（即 dp/dx =0）。故普兰德边界层方程最终可化为：u + u =vxyx2x经过一系列化简之后，详见化工传递过程基础p163，164 页可得一仅仅为 函数的三阶常微分方程如下：+2 =0ff它的边界条件为：=0， =0f=0， =0， =1f满足边界条件的解为一穷级数如下： 8168452 0427.10972.1093.4160. f解：
12、Matlab 解题方法如下（其中 t= ）：已知： 当 =0 时， =0， =0，当 ， =1tftf所以，令 = ， = ， ；1y2 3fy则， = ， = ， = 15.0y我们已经知道了当 ， ， ，但 的值我们不知道，t)(y0)(2)(3我们可以从 时等于 1，可以采取初步估计一个 的值，计算出 ，当2y)0(3y)(2y=1 时，就得到一个合适的 的值。为此，我们采用 ，求)(2y )0(3y 1,xabsw，为此编写如下程序：该程序调用一函数 fminsearch，如下：)0(3minwxfunction f=p164find（x）t,y = ode45(p164,t,0;0;
13、x);n1=n(1,1);m1=y(n1,2); %m1=y2(infinite)f=abs(m1-1);在 Matlab 命令窗口键入如下命令：x0=0.1;t=0:0.2:10; size(t)ans = 1 51将 n1=51 代入 fminsearch 程序，运行， x=fminsearch(p164find,x0)x =0.3320在编写 p164 函数求出 与 的关系，程序如下：ffunction dydt = p164(t,y)dydt=y(2);y(3);-y(1)*y(3)/2;在 Matlab 命令窗口输入：t,y = ode45(p164,t,0;0;x)结果如下：y =
14、 0 0 0.33200.0066 0.0664 0.33200.0266 0.1328 0.33140.0597 0.1989 0.33010.1061 0.2647 0.32740.1656 0.3298 0.32300.2379 0.3937 0.31660.3229 0.4562 0.30780.4203 0.5167 0.29660.5294 0.5747 0.28290.6500 0.6297 0.26670.7811 0.6813 0.24830.9222 0.7289 0.22811.0724 0.7724 0.20641.2309 0.8115 0.18401.3967 0.
15、8460 0.16141.5690 0.8760 0.13911.7468 0.9017 0.11791.9294 0.9233 0.09802.1159 0.9411 0.08012.3056 0.9555 0.06422.4979 0.9669 0.05062.6922 0.9758 0.03902.8881 0.9826 0.02953.0852 0.9877 0.02193.2831 0.9915 0.01593.4817 0.9942 0.01143.6807 0.9961 0.00793.8801 0.9975 0.00544.0797 0.9983 0.00374.2794 0.
16、9989 0.00244.4793 0.9993 0.00164.6791 0.9996 0.00104.8791 0.9997 0.00065.0790 0.9998 0.00045.2790 0.9999 0.00025.4790 0.9999 0.00015.6790 1.0000 0.00015.8790 1.0000 0.00006.0790 1.0000 0.00006.2790 1.0000 0.00006.4790 1.0000 0.00006.6790 1.0000 0.00006.8789 1.0000 0.00007.0789 1.0000 0.00007.2789 1.
17、0000 0.00007.4789 1.0000 0.00007.6789 1.0000 0.00007.8789 1.0000 -0.00008.0789 1.0000 0.00008.2789 1.0000 0.0000其中第一列为 的值，第二列为 的值，第三列为 的值。fff例 1 边界层能量方程的精确解，有了速度边界层方程的解，即 8168452 0427.10972.1093.460. f就可以求出 、 与 、 的关系，即速度分布关系，于是据此关系便可对边界层能量xuy方程求解，边界层能量方程可写为： 2ytatuxty边界条件为： ，0y0t，x0t解：Matlab 解法如下：%化
18、工传递过程基础p166 方程(9-29)% f+0.5*pr*ef*f=0% t=0: f=0;% t=infinite f=1% 令 y1=f% y2=f% 则原方程变为：% y2+0.5*pr*ef*y2=0% 则：y1=y2; （2）% y2=-0.5*pr*ef*y2; （2）% t=0: y1(0)=0, y2(0)=?, y1(infinite)=1% 但是从 y1(infinite)=1,可以采取初步估计一个 y2(0)的值，% 计算出 y2(infinite),当 y2(infinite)=1 时，就得到合适的 y2(0)值。% 为此，采用 w(y2(0)=abs(y2(inf
19、inite,x)-1)% 求 x=min w(y2(0).见 p166find.m% 然后，再用所求到的 y2(0)，求解微分方程组（2）% 本函数的调用方法：% t,y = ode45(p166,t,0;x);function dydt = p166(t,y)pr=0.016;%注意计算机的数值表达范围，下面的原数值的表达分开做乘方和限除法，以免溢出。ef=(0.16603*t2)-4.5943*t2/100*t2/100*t+2.4972*t2/100*t2/100*t4/100-1.4277*t3/1000*t3/1000*t3/100*t2;dydt=y(2) ;-0.5*pr*ef*y(2);
