《2019第三章数值计算ppt课件》由会员分享，可在线阅读，更多相关《2019第三章数值计算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 计算定积分函数介绍,前面的做法,例：用梯形法计算下面定
4、积分 ( 取 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 为计算精度,将自变量看成是向量,抛物线法：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),函数表达式
5、一定要用 单引号 括起来！ 涉及的运算一定要用 数组运算！,例：用 quad 计算定积分：,quad 举例,抛物线法计算二重积分： 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 看成是向
6、量，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,例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
7、(x.(-2); trapz(x,y),梯形法：,抛物线法：, quad(exp(x.(-2),1,2,10e-10),符号积分法：, syms x int(exp(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、od
8、e23t、ode15s、 ode23s、ode23tb,dsolve 求解析解,dsolve 的使用,y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v),其中 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);
9、dy/dx = 2x dsolve(Dy=2*x); dy/dt = 2x,若找不到解析解，则返回其积分形式。,微分方程中用 D 表示对 自变量 的导数，如：,Dy y； D2y y； D3y y,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-x-3*y=0, . x(0)=1, y(0)=0, t) ezp
10、lot(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),只有很少一部分微分方程（组）能求出解析解。大部分微分方程（组）只能利用数值
11、方法求数值解。,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求解器，对于不
12、同的ODE，可以调用不同的求解器。,Matlab提供的ODE求解器,参数说明,odefun 为显式常微分方程，可以用命令 inline 定义，或在函数文件中定义，然后通过函数句柄调用。,fun=inline(-2*y+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),数值求解举例,如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时需用函数文件来定义
13、该常微分方程组。,令 ，则原方程可化为,数值求解举例,先编写函数文件 verderpol.m,function xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);,再编写脚本文件 vdpl.m，在命令窗口直接运行该文件。,clear; global mu;mu=7; y0=1;0; t,x=ode45(verderpol,0,40,y0); plot(t,x(:,1),r-);,问题一 ： 任务分配问题：某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可用台时数分别为800和900，三种工件
14、的数量分别为400、600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？,三、线性规划问题,解 设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：,问题二： 某厂每日8小时的产量不低于1800件。为了进行质量控制，计划聘请两种不同水平的检验员。一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时。检验员每错检一次
15、，工厂要损失2元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？,解 设需要一级和二级检验员的人数分别为x1、x2人, 则应付检验员的工资为：,因检验员错检而造成的损失为：,故目标函数为：,约束条件为：,线性规划模型：,用MATLAB优化工具箱解线性规划,命令：x=linprog（c，A，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）,注意
16、：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.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;
