08实验3插值与数值积分2.pptVIP

* * * * * * * * * * * * * * * * * * * * * * 本次课的主要内容 数值积分的四种方法 1、矩形公式 sum 或 cumsum 2、梯形公式 trapz 3、辛普森公式 quad 4、高斯公式 编程2：gaussinteg Gauss-Lobatto方法 quadl 编程1：simp 3.3 数值积分 数值积分:用数值的方法近似地求一个定积分 (教材P54,参考书P222) 数值分析中主要介绍三种等距节点的求积公式(牛顿-科茨公式) : 1、矩形公式（k=0） a b h h=(b-a)/n fk=f(xk) x1 x2 x3 xi xi+1 xn-1 f1 f2 f3 fi fi+1 f0 fn fn-1 x y o 2、梯形公式（k=1） 3、辛甫森公式（k=2） 已知n+1对节点数据(xi , yi)(i=0,1,…n)，求积分. trapz(y) 按梯形公式计算定积分（单位步长）。 trapz(x,y) x , y同长度，输出 y 对 x 的按梯形公式计算的积分 （变步长）。 quad(fun,a,b) 用辛甫森（2阶）公式计算函数fun在区间 [ a, b]的积分，自动选择步长。 quad(fun,a,b,tol) 与上同，但指定了相对误差 tol。 §1 MATLAB作数值积分 一 用于数值积分的几种命令： sum(x) 输入数组x，输出为x的和，用于矩形公式求积分。 通常还要乘以等分小区间的长度(b-a)/n cumsum(x) 返回与x一样长的向量。此向量的第n个元素为x的 前n个元素之和.用于矩形公式求积分，同样要乘 (b-a)/n（参考书P223） quadl(‘fun’,a,b,tol) 用自适应Gauss-Lobatto公式计算，精度更高。 上机熟悉cumsum,sum,trapz,quad(x,y,a,b) function s=simp(x,y) if mod(length(x),2)==0 error(数据点必须为奇数) end n=length(x);m=(n-1)/2; h=(x(n)-x(1))/2/m; s1=0; s2=0; for i=1:m s1=s1+y(2*i); end for j=1:m-1 s2=s2+y(2*j+1); end s=(y(1)+y(n)+4*s1+2*s2)*h/3 根据辛甫森求积公式 当被积函数不是解析表示时, 比如离散数据表表示的函数 通常就用这个函数按辛甫森 公式计算积分。 将[a,b]区间等分成2m个小区间： 每2个区间构成一个组,即： 每一组有3个节点，如第一组： 由3个节点构造一个拉格朗日2次多项式： 返回 令f(x)= xk ，用(11)式计算 我们不妨只考虑 二 高斯（Gauss)求积公式 各种近似求积公式都可以表示为 (11) 若对于k=0,1,...,m 都 有In = I ，而当k =m+1时，In ≠ I ，则称In 的代数精度为m 。梯形公式代数精度为 1，辛甫森公式的代数精度为 3。 下面介绍的是取消对区间等分的限制，n 给定后同时确定节点xi和系数Ai，使代数精度尽可能高的所谓高斯公式。 我们先考虑节点数为2 而使用(11)计算的积分近似值有代数精度为3。 而构造代数精度为3的形如 G2=A1 f(x1)+ A2 f(x2) (12) 的求积公式。 成立，依次将f(x)= 1, x, x 2, x 3代入，即可得到确定A1,A2 ,x1 ,x2 的方程组。 如何选择节点xi 和系数Ai ，使(11)计算的精度更高？ 对于 f (x)=1, x, x 2, x 3，应该有 将 f (x) =1, x , x 2 , x 3 依次代入，得： 提高精度可以通过增加节点数n ，(11)的代数精度可达到2n-1, 但增加了解高维线性方程组的难度，实用价值不大；另一方法是将区间分小，在小区间上用G2。 n=2的高斯公式为： 将区间(a, b)作m等分，记h=(b-a)/m，xk=a+kh, k=0,1,...,m, 作变换 其中 常用的高斯公式就是： 下面我们来编写M文件，应用高斯公式计算定积分。 将x∈ [xk-1, xk]化为