数值微积分、微分方程

数值微积分,数值微(差)分数值积分数值定积分数值重积分,主要内容,1、数值差分n维向量x=(x1,x2,xn)的差分定义为n-1维向量x=(x2-x1,x3-x2,xn-xn-1)。,Dx=diff(x)如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分。,调用格式：,数值差分,DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。例如：,V=vander(1:6);V=V(:,end:-1:1)%生成范德蒙行列式的矩阵DV=diff(V)%按列计算V的一阶差分DV=diff(V,1,2)%按行计算V的一阶差分DV=diff(V,2)%按列计算V的二阶差分DV=diff(V,2,2)%按行计算V的二阶差分,求下列定积分,数值积分,求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,1、梯形法：trapz,trapz(x,y)x为分割点（节点）组成的向量，y为被积函数在节点上的函数值组成的向量。,Matlab近似计算定积分的相关函数,Matlab计算定积分函数介绍,做法思想,例：用梯形法计算下面定积分(取n=100),解：,a=0,b=1,n=100,yi=f(xi)=1/(1+xi2),x=0:1/100:1;y=1./(1+x.2);trapz(x,y),trapz函数,trapz(x,1./(1+x.2),trapz举例,quad(f,a,b,tol)f=f(x)为被积函数，a,b为积分区间，tol为计算精度,将自变量看成是向量,2、抛物线法：quad,抛物线法,解：,quad(1./(1+x.2),0,1),quad(1./(1+x.2),0,1,10e-10),quad(1./(1+x.2),0,1,10e-16),函数表达式一定要用单引号括起来！涉及的运算一定要用数组运算！,例：用quad计算定积分：,quad举例,建立函数法求定积分,例求定积分：(1)建立被积函数文件fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)调用数值积分函数quad求定积分。S,n=quad(fesin,0,3*pi),n为迭代次数。,洛巴托数值积分,调用格式：Q=quadl(FUN,a,b,TOL),其中参数的含义和quad函数相似,不是数字1,高精度Lobatto数值积分：quadl,例：用函数trapz、quad和quadl分别数值积分将计算结果精确值进行比较.,应用举例,变限积分求导举例,例：求,symsxtF=int(exp(t)*sin(2+sqrt(t3),x,x2);dF=diff(F),综合应用举例,例：求下列函数所围成的平面区域D的面积S，并作出图形：,x=-1/2:0.001:3/2;F1=sin(x);F2=cos(x);plot(x,F1,b-,x,F2,g-),axis(-1,pi/4+1,-1.3,1.3)xlabel(x),ylabel(y),fill(x;x,F1;F2,r)title(y=sinx,y=cosx和x=-0.5及x=1.5所围成的图形),symsxf1=cos(x)-sin(x);f2=-f1;S1=int(f1,x,-0.5,pi/4);S2=int(f2,x,pi/4,1.5);S=S1+S2,Sj=double(S),抛物线法计算二重积分：dblquad,dblquad(f,a,b,c,d,tol),tol为计算精度，若不指定，则缺省精度为10-6,f(x,y)可以由inline定义，或通过一个函数句柄传递,a,b是第一积分变量的积分区间，c,d是第二积分变量的积分区间,按字母顺序，大写字母排在小写字母的前面,二重积分的计算,f=inline(4*x*y+3*y2);I=dblquad(f,-1,1,0,2),f(x,y)中关于第一自变量的运算是数组运算，即把x看成是向量，y看成是标量。也可以全部采用数组运算,例1：计算二重积分,dblquad举例,f(x,y)中关于第一自变量的运算是数组运算，即把x看成是向量，y看成是标量。也可以全部采用数组运算,例2：计算二重积分,dblquad(inline(4*x*y+3*x2),-1,1,0,2),dblquad(inline(4*x*y+3*x.2),-1,1,0,2),X,dblquad举例,例：计算二重积分,dblquad(x,y)4*x*y+3*x.2,-1,1,0,2),指定x、y分别是第一和第二积分变量,dblquad(inline(4*x*y+3*x.2),-1,1,0,2),被积函数f(x,y)的另一种定义方法：匿名函数,dblquad举例,抛物线法：,dblquad(inline(x+y2),0,2,-1,1),符号积分法：,f=int(x+y2,y,-1,1);int(f,0,2),数值实验,例：用Matlab函数近似计算二重积分,其中,利用梯形法，先将a,b区间n等分，hx=(b-a)/n,xi=a+ihx,i=0,1,n,重积分的数值计算可通过单积分组合计算,一般区域上的二重积分,非矩形区域重积分dblquad2.m,functionS=dblquad2(fun,a,b,clo,dhi,n)ifnarginy=dsolve(Dy+2*x*y=x*exp(-x2),x)symsx;diff(y)+2*x*y-x*exp(-x2),dsolve的使用,几点说明,如果省略初值条件，则表示求通解；,如果省略自变量，则默认自变量为t,dsolve(Dy=2*x,x);%dy/dx=2xdsolve(Dy=2*x);%dy/dt=2x,若找不到解析解，则返回其积分形式。,微分方程中用D表示对自变量的导数，如：,Dyy；D2yy；D3yy,dsolve举例,例2：求微分方程在初值条件下的特解，并画出解函数的图形。,y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)ezplot(y);,dsolve举例,例3：求微分方程组在初值条件下的特解，并画出解函数的图形。,x,y=dsolve(Dx+5*x+y=exp(t),Dy-t-3*y=0,.x(0)=1,y(0)=0,t)ezplot(x,y,0,1.3);,注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。,dsolve举例,例：,x,y=dsolve(Dx+5*x=0,Dy-3*y=0,.x(0)=1,y(0)=1,t),r=dsolve(Dx+5*x=0,Dy-3*y=0,.x(0)=1,y(0)=1,t),这里返回的r是一个结构类型的数据,r.x%查看解函数x(t)r.y%查看解函数y(t),只有很少一部分微分方程（组）能求出解析解。大部分微分方程（组）只能利用数值方法求数值解。,dsolve的输出个数只能为一个或与方程个数相等。,Matlab函数数值求解,T,Y=solver(odefun,tspan,y0),其中y0为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。solver为Matlab的ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。,Matlab提供的ODE求解器,参数说明及举例,odefun为显式常微分方程，可以用命令inline定义，或在函数文件中定义，然后通过函数句柄调用。,fun=inline(-2*y2+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);,注：也可以在tspan中指定对求解区间的分割，如：,x,y=ode23(fun,0:0.1:0.5,1);%此时x=0:0.1:0.5,T,Y=solver(odefun,tspan,y0),数值求解举例,如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义该常微分方程组。,令，则原方程可化为,数值求解举例,先编写函数文件verderpol.m,functionxprime=verderpol(t,x)globalmu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);,再编写脚本文件vdpl.m，再运行该文件。,clear;globalmu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0);plot(t,x(:,1),r-);figure(2)plot(t,x(:,2),r-),再说定义函数文件,T,Y=ode45(odefun,tspan,y0)T,Y=ode23(odefun,tspan,y0),odefun？,ode45、ode23能解什么样的ODE？,再说定义函数文件,T,Y=ode45(odefun,tspan,y0)T,Y=ode23(odefun,tspan,y0),ode45、ode23等函数用于求解显式常微分方程当是向量函数时，所对应的上面的方程即为微分方程组,odefun,举例说明,fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);,例：求初值问题的数值解。,解法一：使用inline定义微分方程odefun,odefun为方程右端项f(t,y)可以用inline定义（只适合于单个方程的情形）通过函数文件定义，然后用函数句柄调用（适合所有情形）,注：自变量必须在前面，因变量在后面！,举例说明（单个方程）,functiondy=myfun1(x,y)dy=-2*y+2*x2+2*x;,解法二：通过函数文件定义微分方程odefun,1、先编写函数文件myfun1.m,clear;x,y=ode23(myfun1,0,0.5,1);,2、编写主文件main1.m,举例说明（方程组）,解：此时只能通过函数文件定义微分方程odefun,例：求,的数值解。,functiondy=myfun2(t,y)dy=zeros(3,1);%dymustbeacolumnvector!dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,1、先编写函数文件myfun2.m,clear;T,Y=ode45(myfun2,0,12,0,1,1);,2、编写主文件main1.m,dy=y(2)*y(3);-y(1)*y(3);.-0.51*y(1)*y(2);,思考,functiondy=myfun2(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,1、函数文件myfun2.m,能不能写成下面形式？,functiondy=myfun2(t,x,y,z)dy=zeros(3,1);dy(1)=y*z;dy(2)=-x*z;dy(3)=-0.51*x*y;,X,说明,odefun,变量属性必须一一对应！,functiondy=myfun2(t,y),如果是常微分方程组，y就是列向量！,dy必须是列向量，长度为方程组的个数，通常与y的长度相同,函数中的输入参数和输出参数是形参，名字可以任意取，但必须满足上述条件。即输入参数有两个，第一个表示自变量，第二个是由因变量组成的列向量，输出参数必须是列向量。,例,functiondy=myfun3(t,y)dy=zeros(2,1);dy(1)=y(2)+t;dy(2)=t-2;,例：解初值问题：，,functionyprime=myfun3(t,y)yprime=y(2)+t;t2;,functionyprime=myfun3(x,y)yprime=y(2)+x;x2;,clear;T,Y=ode45(myfun3,0,10,1,1);,2、主文件main3.m,1、函数文件myfun3.m,高阶常微分方程,高阶常微分方程,例：VanderPol初值问题,令，则原方程可化为,一阶常微分方程组,变量替换化为,参数怎么处理？,用全局变量传递,思考,对于高阶微分方程组的数值求解。应用如何处理？,导弹追踪问题设位于坐标原点的甲舰向位于x轴上的A(1,0)处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度沿平行于y轴的直线行驶，导弹的速度是5，求导弹运行的曲线。当乙舰行驶多远时，导弹将它击中。,综合举例,设任意时刻t，乙舰的坐标为，导弹的坐标为。导弹速度恒为，从原点射出，且速度是横向和纵向距离的一阶导数的矢量和，则有：