《《数值计算》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《数值计算》PPT课件.ppt（49页珍藏版）》请在三一办公上搜索。
1、第三章 数值运算基础,定积分的近似计算微分方程的求解线性规划问题的求解,一、定积分的近似计算,定积分计算的基本公式是牛顿莱布尼兹公式。但当被积函数的原函数不知道时，如何计算？这时就需要利用近似计算。特别是在许多实际应用中，被积函数甚至没有解析表达式，而是一条实验记录曲线，或一组离散的采样值，此时只能用近似方法计算定积分。,主要研究定积分的三种近似计算算法：矩形法、梯形法和抛物线法。同时介绍 Matlab 计算定积分的相关函数。,矩形法,定积分的定义：,矩形法,n 充分大，x 充分小,定积分的近似：,左点法,右点法,中点法,点 可以任意选取，常见的取法有：左端点，右端点 和中点。,步长,节点,右
2、点法：,中点法：,左点法：,左点法、右点法和中点法,解：,矩形法举例,a=0,b=1,n=100,例：用不同的矩形法计算下面的定积分(取 n=100)，并比较这三种方法的相对误差。,左点法：,右点法：,中点法：,(i=0,1,2,.,100),理论值：,左点法相对误差：,误差分析,矩形法举例,右点法相对误差：,中点法相对误差：,不同的方法有不同的计算精度,有没有更好的近似计算定积分的方法?,定积分几何意义,曲边小梯形的面积可以由直边小梯形的面积来近似,整个曲边梯形的面积：,梯形法,如果我们 n 等分区间 a,b，即令：,则,=,梯形公式,梯形法,梯形公式与中点公式有什么区别?,解：,=,例：用
3、梯形法计算下面定积分(取 n=100)，并计算相对误差,梯形法举例,a=0,b=1,n=100,f(x)=1/(1+x2),相对误差：,2n 等分区间 a,b，得,该直线用抛物线代替，计算精度是否会更好？,计算每个节点上的函数值：,抛物线法,在区间 x0,x2 上，用过以下三点,的抛物线来近似原函数 f(x)。,梯形法：trapz,trapz(x,y)x 为分割点（节点）组成的向量，y 为被积函数在节点上的函数值组成的向量。,Matlab 近似计算定积分的相关函数,Matlab 计算定积分函数介绍,前面的做法,例：用梯形法计算下面定积分(取 n=100),解：,a=0,b=1,n=100,yi
4、=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 为计算精度,将自变量看成是向量,抛物线法：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 举例,抛
5、物线法计算二重积分：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 看成是标量。也可以全部采用数组运算,例2：计算二重积分,dblquad(inline(4*x*y+3*x2),-1,
6、1,0,2),dblquad(inline(4*x*y+3*x.2),-1,1,0,2),X,例1：计算二重积分,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 举例,x=1:0.001:2;y=exp(x.(-2);trapz(x,y),梯形法：,抛物线法：,quad(exp(x.(-2),1,2,10e-10),符号积分法：,syms x int(ex
7、p(x(-2),x,1,2),例 1：用 Matlab 函数近似计算积分,数值实验,抛物线法：,dblquad(inline(x+y2),0,2,-1,1),符号积分法：,f=int(x+y2,y,-1,1);int(f,0,2),数值实验,例 2：用 Matlab 函数近似计算二重积分,Matlab 微分方程,用 Maltab自带函数 解微分方程,求解析解：dsolve,求数值解：ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tb,dsolve 求解析解,dsolve 的使用,y=dsolve(eq1,eq2,.,cond1,cond2,.,v),
8、其中 y 为输出，eq1、eq2、.为微分方程，cond1、cond2、.为初值条件，v 为自变量。,例 1：求微分方程 的通解，并验证。,y=dsolve(Dy+2*x*y=x*exp(-x2),x),syms x;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 表示对 自变量 的导数，如：,Dy y；D2y y；D3y y,dsolve
9、举例,例 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-x-3*y=0,.x(0)=1,y(0)=0,t)ezplot(x,y,0,1.3);,注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。,dsolve 举例,例：,x,y
10、=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在数值求解时自动对求解区间进行分割，
11、T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。solver 为Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的 ODE 问题，因此MATLAB 提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。,Matlab提供的ODE求解器,参数说明,odefun 为显式常微分方程，可以用命令 inline 定义，或在函数文件中定义，然后通过函数句柄调用。,fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode2
12、3(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,function xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1-x(1)2)*x(2)-x(1);,再编写脚本文件 vdpl.
13、m，在命令窗口直接运行该文件。,clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0);plot(t,x(:,1),r-);,问题一：任务分配问题：某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400、600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？,三、线性规划问题,解 设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量
14、分别为x4、x5、x6。可建立以下线性规划模型：,问题二：某厂每日8小时的产量不低于1800件。为了进行质量控制，计划聘请两种不同水平的检验员。一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时。检验员每错检一次，工厂要损失2元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？,解 设需要一级和二级检验员的人数分别为x1、x2人,则应付检验员的工资为：,因检验员错检而造成的损失为：,故目标函数为：,约束条件为：,线性规划模型：,用MATLAB优化工具箱解线性规划,命令：x=linprog（c，A，
15、b）,2、模型：min z=cX,命令：x=linprog（c，A，b，Aeq,beq）,注意：若没有不等式：存在，则令A=，b=.,命令：1 x=linprog（c，A，b，Aeq,beq,VLB，VUB）2 x=linprog（c，A，b，Aeq,beq,VLB，VUB,X0）,注意：1 若没有等式约束:,则令Aeq=,beq=.2其中X0表示初始点,4、命令：x,fval=linprog()返回最优解及处的目标函数值fval.,解 编写M文件xxgh1.m如下：c=-0.4-0.28-0.32-0.72-0.64-0.6;A=0.01 0.01 0.01 0.03 0.03 0.03;0
16、.02 0 0 0.05 0 0;0 0.02 0 0 0.05 0;0 0 0.03 0 0 0.08;b=850;700;100;900;Aeq=;beq=;vlb=0;0;0;0;0;0;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),解:编写M文件xxgh2.m如下：c=6 3 4;A=0 1 0;b=50;Aeq=1 1 1;beq=120;vlb=30;0;20;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),例3 问题一的解答,编写M文件xxgh3.m如下:f=13 9 10 11 12 8;A=0.4
17、1.1 1 0 0 0 0 0 0 0.5 1.2 1.3;b=800;900;Aeq=1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1;beq=400;600;500;vlb=zeros(6,1);vub=;x,fval=linprog(f,A,b,Aeq,beq,vlb,vub),结果:x=0.0000 600.0000 0.0000 400.0000 0.0000 500.0000fval=1.3800e+004 即在甲机床上加工600个工件2,在乙机床上加工400个工件1、500个工件3，可在满足条件的情况下使总加工费最小为13800。,例2 问题二的解答,改写为
18、：,编写M文件xxgh4.m如下：c=40;36;A=-5-3;b=-45;Aeq=;beq=;vlb=zeros(2,1);vub=9;15;%调用linprog函数：x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),结果为：x=9.0000 0.0000fval=360即只需聘用9个一级检验员。,注：本问题应还有一个约束条件：x1、x2取整数。故它是一个整数线性规划问题。这里把它当成一个线性规划来解，求得其最优解刚好是整数：x1=9，x2=0，故它就是该整数规划的最优解。若用线性规划解法求得的最优解不是整数，将其取整后不一定是相应整数规划的最优解，这样的整数规划应用专门的方法求解。,实验练习,某厂生产甲乙两种口味的饮料,每百箱甲饮料需用原料6千克,工人10名,可获利10万元;每百箱乙饮料需用原料5千克,工人20名,可获利9万元.今工厂共有原料60千克,工人150名,又由于其他条件所限甲饮料产量不超过8百箱.问如何安排生产计划,即两种饮料各生产多少使获利最大.进一步讨论:1)若投资0.8万元可增加原料1千克,问应否作这项投资.2)若每百箱甲饮料获利可增加1万元,问应否改变生产计划.,
