Matlab-求解化工常微分方程和偏微分方程课件.ppt

Matlab 求解化工常微分方程和偏微分方程,方利国,Matlab 求解化工常微分方程和偏微分方程,1、常微分方程（组）求解1.1 问题描述及Matlab调用命令1.2 初值问题求解1.3 边值问题求解1.4 加权问题求解（自学内容）2、偏微分方程（组）求解2.1 问题描述及一维动态PDE方程求解2.2 二维求解,1、常微分方程（组）求解,1.1 问题描述及Matlab调用命令 常微分方程：（初值问题）常微分方程：（两点边值问题）,1、常微分方程（组）求解,1.1 问题描述及Matlab调用命令 微分方程组：,1、常微分方程（组）求解,1.1 问题描述及Matlab调用命令 高级微分方程：,1、常微分方程（组）求解,1.1 问题描述及Matlab调用命令 Matlab 调用命令：ODE45:4-5 阶龙格库塔法（非刚性）ODE23:2-3 阶龙格库塔法（非刚性）ODE113:可变D-B-M法（非刚性）ODE15S:基于数值差分的可变阶方法法（刚性）ODE23S、ODE23t、ODE23tb（刚性）,1、常微分方程（组）求解,通用调用格式：x,y=ode*(odefun,xspan,y0,)X:自变量向量，在实际调用时取名不一定要用x，也可以用其他名称，只要前后一致即可。Y:应变量向量，在实际调用时取名不一定要用Y，也可以用其他名称，只要前后一致即可。*:根据不同的问题调用不同格式，如45，23sodefun:自定义函数的函数名，该函数为Xspan：自变量的积分限，xa,xb,也可以是离散点，x0,x1,x2,xfy0：应变量向量的初值:可以没有该选项，如有，具体应用见下面的实际例子,1.2 初值问题求解,例1：已知某高温物体其温降过程符合以下规律，其中温度T的单位为K，时间 的单位为分钟，零时刻高温物体的温度为2000K，以1分钟作为时间步长，请计算零时刻以后每隔1分钟至170分钟的温度。,单个微分方程,function xODEs%铁球从2000K降温曲线,在7.0 版本上调试通过%由华南理工大学方利国编写，2012年2月29日%欢迎读者调用，如有问题请告知 clear all；clcy0=2000;x1,y1=ode45(f,0:1:170,y0);%0到170分钟，每分钟一个计算点x2,y2=ode23(f,0:1:170,y0);plot(x1,y1,r-)xlabel(时间，M)ylabel(温度，K)hold ondisp(Results by using ode45()：)disp(x y(1)disp(x1 y1)disp(Results by using ode23()：)disp(x y(2)disp(x2 y2)plot(x2,y2,k:)%-function dy=f(x,y)%定义降温速率的微分方程%dy=0.04*y(1)-100;dy=-0.04*exp(0.001*(y(1)-300)*(-300+y(1);,Results by using ode45()：x y(1)1.0e+003*Columns 1 through 13 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 2.0000 0.8788 0.6133 0.4880 0.4181 0.3762 0.3498 0.3328 0.3218 0.3145 0.3097 0.3065 0.3043 Columns 14 through 18 0.1300 0.1400 0.1500 0.1600 0.1700 0.3029 0.3019 0.3013 0.3009 0.3006Results by using ode23()：x y(2)1.0e+003*Columns 1 through 13 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 2.0000 0.8779 0.6124 0.4873 0.4176 0.3756 0.3493 0.3323 0.3213 0.3140 0.3092 0.3061 0.3041 Columns 14 through 18 0.1300 0.1400 0.1500 0.1600 0.1700 0.3027 0.3018 0.3012 0.3008 0.3005,1.2 初值问题求解,该问题相当与一个自变量，两个应变量问题，已知初值及微分表达式，可以利用ODE45求解。,微分方程组求解,程序代码,function uvDEs%微生物消亡问题计算,在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月12日%欢迎读者调用，如有问题请告知 clear all;Clcy0=1.6 1.2;x1,y1=ode45(f,0:0.1:10,y0);%0到3分钟，每0.1分钟一个计算点u=y1(:,1);v=y1(:,2);plot(x1,u,r-)xlabel(时间，M)ylabel(微生物浓度)hold onplot(x1,v,k:)disp(Results by using ode45()：)disp(x u v)disp(x1 y1)%-function dy=f(x,y)%定义降温速率的微分方程f1=0.09*y(1)*(1-y(1)/20)-0.45*y(1)*y(2);f2=0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2);dy=f1;f2;,1.2 初值问题求解,例3：,当X较大时，两种方法计算结果有较大不同，为什么？,单个微分方程有零点问题？,function L43ODEs%在7.0 版本上调试通过%由华南理工大学方利国编写，2012年2月29日%欢迎读者调用，如有问题请告知 clear allclcy0=1;x1,y1=ode45(f,0:0.05:10,y0);%0到10，每0.05间隔一个计算点x2,y2=ode23(f,0:0.05:10,y0);%0到10，每0.05间隔一个计算点plot(x1,y1,r-)xlabel(x)ylabel(y)hold ondisp(Results by using ode45()：)disp(x y(1)disp(x1 y1)disp(Results by using ode23()：)disp(x y(2)disp(x2 y2)plot(x2,y2,b-)%-function dy=f(x,y)%定义微分方程dy=y2*cos(x);,计算值,x y(1)Columns 1 through 13 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 1.0000 1.0526 1.1109 1.1757 1.2479 1.3287 1.4195 1.5218 1.6378 1.7698 1.9210 2.0951 2.2970 Columns 14 through 26 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 2.5329 2.8107 3.1411 3.5381 4.0206 4.6152 5.3598 6.3078 7.5435 9.1928 11.4614 14.7283 19.5991 Columns 27 through 29 1.3000 1.3500 1.4000 27.4283 41.2030 68.6630,高阶微分方程求解,求解思路：将高阶微分方程通过变量转换，转变成一级微分方程组进行求解。,例4：,高阶微分方程求解程序及解,Columns 1 through 13 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 0 0.1002 0.2015 0.3052 0.4129 0.5263 0.6476 0.7789 0.9230 1.0827 1.2614 1.4627 1.6908 1.0000 1.0053 1.0227 1.0544 1.1026 1.1698 1.2588 1.3723 1.5136 1.6861 1.8935 2.1398 2.4296 Columns 14 through 21 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 1.9502 2.2461 2.5841 2.9705 3.4123 3.9171 4.4936 5.1513 2.7677 3.1595 3.6110 4.1286 4.7196 5.3918 6.1541 7.0159,刚性方程求解,有些微分组的系数变化很大，这时用ODE45就很难收敛求解，这时可用专门解决此类微分方程的ODE23S来求解，需要注意的是在解的图像绘制时，也需要考虑数值的波动幅度很大，需要引入对数坐标。,例5：dy1/dx=-0.03*y1+1e4*y2*y3 dy2/dx=0.03*y1-2e4*y2*y3-5e7*y22 dy3/dx=5e7*y22 y1(0)=1,y2(0)=0,y3(0)=0,function gangxinDEs%刚性问题计算,在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知%dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y22%dy3=5e7*y22clear allclcfigurexspan=0 6*logspace(-6,6);y0=1 0 0;x1,y1=ode15s(f,xspan,y0);%用 ode45计算刚性方程，可能有问题u=y1(:,1);v=1e4*y1(:,2);w=y1(:,3);semilogx(x1,u,r-,linewidth,2)xlabel(x)ylabel(1e4*v)hold onsemilogx(x1,v,k:,linewidth,2)hold onsemilogx(x1,w,g-,linewidth,2)gridaxis(10(-10)1010-0.2 1.2)legend(u,v,w)disp(Results by using ode45()：),disp(x u v w)format longdisp(x1 y1)%-function dy=f(x,y)%定义微分方程f1=-0.03*y(1)+1e4*y(2)*y(3);f2=0.03*y(1)-2e4*y(2)*y(3)-5e7*y(2)2;f3=5e7*y(2)2;dy=f1;f2;f3;,1.3 边值问题求解,边值问题相对于初值问题而言，多了一个端点的约束，如果在高阶或微分方程组中端点约束过多，微分方程组可能无解，端点约束有一定限制。可以通过建立离散的方程组，再利用ODE45进行求解，但可以利用MATLAB的专用工具求解最好。下面介绍ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通过实际例子介绍这些内部函数的功能。,1.3 边值问题求解,solinit=bvpinit（x,yinit):产生在初始网格上的初始解，以便bvp4c调用，其中x为自变量网格，yinit为对应函数的初值。sol=bvp4c(odefun,BCfun,solinit,)Bcfun:为定义边界条件方程，Bcfun(ya,yb),其中ya、yb分别表示左右边界。其他符号意义同上。deval(sol,xint):计算任意点处的函数值。,1.3 边值问题求解,例6：,程序,function BVP4c1%求解两点边值问题的示例在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知 clear allclca=0;b=10;solinit=bvpinit(linspace(a,b,101),0 0);sol=bvp4c(ODEfun,BCfun,solinit);format shortx=0:0.1:10;y=deval(sol,x);y1=y(1,:)y2=y(2,:)plot(x,y1,r-)xlabel(x)ylabel(y)hold ongridplot(x,y2,k:)legend(y,dy),disp(Results by using bvp4c：)disp(x y dy)disp(x1 y1)%-function dy=ODEfun(x,y)f1=y(2);f2=0.05*(1+x2)*y(1)+2;dy=f1;f2;%-function bc=BCfun(ya,yb)bc=ya(1)-40;yb(1)-80;,Results by using bvp4c：x y dy Columns 1 through 13 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 39.9998 38.1722 36.3840 34.6347 32.9244 31.2531 29.6214 28.0298 26.4791 24.9702 23.5038 22.0810 20.7025-18.4739-18.0779-17.6872-17.2985-16.9088-16.5158-16.1176-15.7125-15.2996-14.8781-14.4476-14.0080-13.5597,计算结果,1.3 边值问题求解,例7：,主要程序,function dy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=f1;f2;%-function bc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-exp(4);,function BVP4c2%求解两点边值问题的示例在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知 clear allclca=0;b=2solinit=bvpinit(linspace(a,b,21),0 0);sol=bvp4c(ODEfun,BCfun,solinit);format shortx=0:0.1:2;y=deval(sol,x);y1=y(1,:);y2=y(2,:);y3=exp(2*x);plot(x,y1,r-,linewidth,2)xlabel(x)ylabel(y)hold ongridplot(x,y2,k:,linewidth,2)hold onplot(x,y3,b:,linewidth,2)legend(y,dy,real-y)disp(Results by using bvp4c：)disp(x y dy)disp(x;y1;y2)%-function dy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=f1;f2;%-function bc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-exp(4);,2、偏微分方程（组）求解,2.1 问题描述 2.2一维动态PDE方程求解2.3 二维稳态PDE方程求解(自学）,问题描述,当A，B，C为常数时，称为拟线性偏微分方程，当A，B，C满足不同条件时，分为三种不同的类型：,不同类型的方程，MATLAB求解时，采用不同或相同的内置函数或工具箱,PDE 求解数学模型通式,m=0，表示平板，m=1表示圆柱，m=2表示球形，f项表示通量项，,s项表示源项，c项为对角阵，元素必须大于等于0才可以求解。,调用方法,sol=pdepe(m,pdefun,iCfun BCfun,xspan,tspan,)Bcfun:为定义边界条件方程iCfun:为定义初始条件方程。其他符号意义同上，注意边界条件必须写成以上形式：,应用策略,Pdepe 内部函数具体应用时，需将实际的偏微分方程对照标准模型，确定c、f、s函数的具体形式及边界条件p、q的具体形式及m 值。,套管动态传热问题,原方程：,代入具体数据，并对长度归一化无量纲处理后，得到以下方程：,对照标准模型，T就是标准模型中的函数u，m=0,a=0,b=1,t0=0,t f=1,标准模型：,初始条件：零时刻所有位置温度为30，即u0=30,边界条件：已知在零位置处任意时间温度为30，在1位置处，偏导为0 pa=u-30,qa=0;pb=0,qb=1,程序清单,function l51clcclear allglobal ua ubalpha=1.0;%cm/sua=30;ub=30;m=0;a=0;b=1;t0=0;tf=1x=linspace(a,b,11);t=linspace(t0,tf,101);sol=pdepe(m,PDEfun,ICfun,BCfun,x,t);,u=sol(:,:,1)%surface plot of the solutionfigure;surf(x,t,u);title(Numerical solution computed with 11 mesh points.);xlabel(Disuance x);ylabel(Time t);ppp=u%solution profile at t=1figure;subplot(2,2,1);plot(x,u(1,:),o-);title(Solutions at t=0.);xlabel(Distance x);ylabel(u(t=0,);grid on;subplot(2,2,2);,程序清单,plot(x,u(50,:),o-);title(Solutions at t=0.5.);xlabel(Distance x);ylabel(u(t=0.5);grid on;subplot(2,2,3);plot(x,u(end,:),o);title(Solutions at t=1.);xlabel(Distance x);ylabel(u(t=1,);grid on;%-function c1,f,s=PDEfun(x,t,u,Du),c1=1f=0.001*Du;s=2*(150-t)-3*Du%-function u0=ICfun(x)u0=30;%-function pa,qa,pb,qb=BCfun(xl,ul,xr,ur,t)global ua ubpa=ul-ua;qa=0;pb=0qb=1,计算结果,分析讨论,