计物05-matlab05-微积分与微分方程数值解

2020/5/30,1,第五讲数值微积分方法,以微分（差分）、积分的数值法以及微分方程数值解为例介绍数值方法,微分（差分）积分的数值法常微分方程数值解特殊微分方程数值方法偏微分方程数值法介绍,2020/5/30,2,5.1差分,数值微分非常困难。积分描述了一个函数的整体或宏观性质，而微分则描述一个函数在一点处的斜率，这是函数的微观性质。因此积分对函数的形状在小范围内的改变不敏感。而微分却很敏感。一个函数小的变化，容易产生相邻点的斜率的大的改变。由于微分这个固有的困难，所以尽可能避免数值微分，特别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这种数据，然后对所得到的多项式进行微分。,2020/5/30,3,*diff*diff(x)对向量,x(2:n)-x(1:(n-1);diff(X)对矩阵列的差分X(2:n,:)-X(1:(n-1),:);diff(X,n)n阶差分;diff(X,n,DIM)沿DIM维n阶差分;,给定一些描述某函数的数据，MATLAB提供了一个计算其非常粗略的微分的函数。这个函数命名为diff，它计算数组中元素间的差分。y=f(x)的微分可近似为：,2020/5/30,4,例1,x=0.1.2.3.4.5.6.7.8.91y=-.4471.9783.286.167.087.347.669.569.489.3011.2;%datady=diff(y)./diff(x);%computedifferencesandusearraydivisionxd=x(1:length(x)-1);%createnewxaxissincedyisshorterthanyplot(xd,dy);title(ApproximateDerivativeUsingDIFF)ylabel(dy/dx),xlabel(x),2020/5/30,5,*gradient*近似梯度fx,fy=gradient(F)各差分间隔为1;fx,fy=gradient(F,H)H数量，即各方向差分间隔;fx,fy=gradient(F,HX,HY)H矢量指定差分间隔;fx,fy,fz=gradient(F,HX,HY,HZ)三维梯度；fx,fy,fz,=gradient(F,)n维梯度;,例2v=-2:0.2:2;x,y=meshgrid(v);z=x.*exp(-x.2-y.2);px,py=gradient(z,.2,.2);contour(v,v,z),holdon,quiver(v,v,px,py),holdoff,2020/5/30,6,函数1quad、quadl、quad8功能数值定积分，自适应Simpleson积分法。格式q=quad(fun,a,b)%近似地从a到b计算函数fun的数值积分，误差为10-6。若给fun输入向量x，应返回向量y，即fun是一单值函数。q=quad(fun,a,b,tol)%用指定的绝对误差tol代替缺省误差。tol越大，函数计算的次数越少，速度越快，但结果精度变小。q=quad(fun,a,b,tol,trace,p1,p2,)%将可选参数p1,p2,等传递给函数fun(x,p1,p2,)，再作数值积分。若tol=或trace=，则用缺省值进行计算。q,n=quad(fun,a,b,)%同时返回函数计算的次数n=quadl(fun,a,b,)%用高精度进行计算，效率可能比quad更好。=quad8(fun,a,b,)%该命令是将废弃的命令，用quadl代替。,5.2数值积分,5.2.1一元函数的数值积分,2020/5/30,7,例3,fun=inline(3*x.2./(x.3-2*x.2+3);Q1=quad(fun,0,2)Q2=quadl(fun,0,2)结果为：Q1=3.7224Q2=3.7224,2020/5/30,8,函数2trapz功能梯形法数值积分格式T=trapz(Y)%用等距梯形法近似计算Y的积分。若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。T=trapz(X,Y)%用梯形法计算Y在X点上的积分。若X为一列向量，Y为矩阵，且size(Y,1)=length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。T=trapz(,dim)%沿着dim指定的方向对Y进行积分。若参量中包含X，则应有length(X)=size(Y,dim)。,2020/5/30,9,例4,X=-1:.1:1;Y=1./(1+25*X.2);T=trapz(X,Y)计算结果为：T=0.5492,2020/5/30,10,函数3cumsum功能元素累积和函数。格式cumsum(X,DIM)%参数DIM指明求和从第DIM维进行。,例5,X=0,1,2;3,4,5cumsum(X,1)ans=012357,X=0,1,2;3,4,5cumsum(X,2)ans=0133712,2020/5/30,11,此外，还有一些求积算法和函数,Guass求积算法%参考1，p238Romberg求积算法%参考1，p240Mote-Carlo求积方法%参考1，p242,2020/5/30,12,函数1dblquad功能矩形区域上的二重积分的数值计算格式q=dblquad(fun,xmin,xmax,ymin,ymax)%调用函数quad在区域xmin,xmax,ymin,ymax上计算二元函数z=f(x,y)的二重积分。输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。q=dblquad(fun,xmin,xmax,ymin,ymax,tol)%用指定的精度tol代替缺省精度10-6，再进行计算。q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)%用指定的算法method代替缺省算法quad。method的取值有quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,)%将可选参数p1,p2,.等传递给函数fun(x,y,p1,p2,)。若tol=，method=，则使用缺省精度和算法quad。,5.2.1二元函数的数值积分,2020/5/30,13,例6,fun=inline(y./sin(x)+x.*exp(y);Q=dblquad(fun,1,3,5,7)计算结果为：Q=3.8319e+003,2020/5/30,14,函数2quad2dggen功能任意区域上二元函数的数值积分格式q=quad2dggen(fun,xlower,xupper,ymin,ymax)%在由xlower,xupper,ymin,ymax指定的区域上计算二元函数z=f(x,y)的二重积分。q=dblquad(fun,xlower,xupper,ymin,ymax,tol)%用指定的精度tol代替缺省精度10-6，再进行计算。q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)%用指定的算法method代替缺省算法。method的取值有缺省算法或用户指定的、与缺省命令有相同调用次序的函数句柄。q=dblquad(fun,xlower,xupper,ymin,ymax,tol,method,p1,p2,)%将可选参数p1,p2,.等传递给函数fun(x,y,p1,p2,)。若tol=，method=，则使用缺省精度和算法。,2020/5/30,15,例7,计算单位圆域上的积分：先把二重积分转化为二次积分的形式：,fun=inline(exp(-x.2/2).*sin(x.2+y),x,y);xlower=inline(-sqrt(1-y.2),y);xupper=inline(sqrt(1-y.2),y);Q=quad2dggen(fun,xlower,xupper,-1,1,1e-4)计算结果为：Q=0.5368603818,2020/5/30,16,微分方程问题算法概述四阶定步长Runge-Kutta算法及MATLAB实现一阶微分方程组的数值解微分方程转换,5.3常微分方程问题的数值解法,2020/5/30,17,5.3.1常微分方程问题算法概述（了解）,2020/5/30,18,微分方程求解的误差与步长问题(了解）,2020/5/30,19,2020/5/30,20,2020/5/30,21,5.3.2四阶定步长Runge-Kutta算法及MATLAB实现(了解),2020/5/30,22,2020/5/30,23,5.3.3一阶微分方程组的数值解,1.四阶五级Runge-Kutta-Felhberg算法（了解）,2020/5/30,24,2.基于MATLAB的微分方程求解函数,2020/5/30,25,2020/5/30,26,2020/5/30,27,【例8】,2020/5/30,28,2020/5/30,29,【例9】,2020/5/30,30,2020/5/30,31,3.MATLAB下带有附加参数的微分方程求解,【例10】,2020/5/30,32,2020/5/30,33,2020/5/30,34,2020/5/30,35,5.3.4微分方程转换（了解）,1.单个高阶常微分方程处理方法,2020/5/30,36,2020/5/30,37,【例11】,2020/5/30,38,2020/5/30,39,2.高阶常微分方程组的变换方法,2020/5/30,40,【例12】,2020/5/30,41,2020/5/30,42,2020/5/30,43,2020/5/30,44,2020/5/30,45,【例13】,2020/5/30,46,2020/5/30,47,【例14】,2020/5/30,48,2020/5/30,49,5.4特殊微分方程的数值解（了解）,刚性微分方程的求解隐式微分方程求解微分代数方程的求解延迟微分方程求解,2020/5/30,50,5.4.1刚性微分方程的求解,【例15】,2020/5/30,51,【例16】,2020/5/30,52,2020/5/30,53,2020/5/30,54,【例17】,2020/5/30,55,2020/5/30,56,2020/5/30,57,5.4.2隐式微分方程求解,【例18】,2020/5/30,58,2020/5/30,59,【例19】,2020/5/30,60,2020/5/30,61,2020/5/30,62,【例20】,2020/5/30,63,2020/5/30,64,5.4.3微分代数方程的求解,2020/5/30,65,【例21】,2020/5/30,66,2020/5/30,67,2020/5/30,68,2020/5/30,69,2020/5/30,70,2020/5/30,71,【例22】,2020/5/30,72,2020/5/30,73,5.4.4延迟微分方程求解,2020/5/30,74,【例23】,2020/5/30,75,2020/5/30,76,【例24】中性延迟微分方程,2020/5/30,77,5.5偏微分方程求解入门(了解）,偏微分方程组求解二阶偏微分方程的数学描述偏微分方程的求解界面应用举例,2020/5/30,78,5.5.1偏微分方程组求解,2020/5/30,79,边界条件的函数描述：,2020/5/30,80,【例7-25】,2020/5/30,81,2020/5/30,82,2020/5/30,83,2020/5/30,84,5.5.2二阶偏微分方程的数学描述,1.椭圆型偏微分方程,2020/5/30,85,2020/5/30,86,2.抛物线型偏微分方程,2020/5/30,87,3.双曲型偏微分方程,2020/5/30,88,4.特征值型偏微分方程,2020/5/30,89,5.5.3偏微分方程的求解界面应用举例,启动偏微分方程求解界面在MATLAB下键入pdetool该界面分为四个部分菜单系统工具栏集合编辑求解区域,1.偏微分方程求解程序概述,2020/5/30,90,2.偏微分方程求解区域绘制,3.偏微分方程边界条件描述,2020/5/30,91,4.偏微分方程求解举例,【例26】,5.时变解的动画显示,2020/5/30,92,6.函数参数的偏微分方程求解,【例27】,2020/5/30,93,本章内容简介,2020/5/30,94,本章介绍了基于MATLAB的差分函数diff()，梯度函数gradient()的用法；介绍了数值积分的几个函数；,2020/5/30,95,引入了数值解的概念，并以最简单的一阶微分方程的Euler算法为例，介绍了一般数值解法的思路并介绍了变步长求解的概念，还介绍了MATLAB下的微分方程数值求解函数ode45()，通过例子演示了该函数的使用方法。微分方程初值函数能直接求解的方程是一阶显式微分方程组，若给出的方程不是这类函数，则需要通过本书介绍的方法选择一组状态变量，将原方程变换成一阶显式微分方程组，以便用给定的求解函数直接求解。若某微分方程模型求解速度极慢，则有可能为刚性方程，需要调用ode15s()等函数来求解，此外，其他