高斯公式积分matlab,三用MATLAB实现定积分计算.PPT

令用式计算我们不妨只考虑二高斯求积公式各种近似求积公式都可以表示为若对于都有而当时则称的代数精度为梯形公式代数精度为辛甫森公式的代数精度为下面介绍的是取消对区间等分的限制给定后同时确定节点和系数使代数精度尽可能高的所谓高斯公式我们先考虑节点数为而使用计算的积分近似值有代数精度为而构造代数精度为的形如的求积公式成立依次将代入即可得到确定的方程组如何选择节点和系数使计算的精度更高对于应该有将依次代入得提高精度可以通过增加节点数的代数精度可达到但增加了解高维线性方程组的难度实用价值不大另一方法是将区间

令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]化为t∈[-1, 1] function s=gaussinteg(f,a,b,m) s=0; h=(b-a)/m;x=a:h:b; for j=1:m z1(j)=(x(j)+x(j+1))/2+h/2/sqrt(3); z2(j)=(x(j)+x(j+1))/2-h/2/sqrt(3); s=s+feval(f,z1(j))+feval(f,z2(j)); end s=s*h/2; function y=jifenfun(x) y=exp(x).*sin(x); return 定义高斯积分函数M文件： s=gaussinteg('jifenfun', ... 0,2*pi,1000) s = -267.2458 在命令窗口运行 定义任何一个被积函数 根据 Gauss-lobatto是改进的高斯积分方法,采取自适应求积方法 三 用MATLAB实现定积分计算： ⑴ 矩形公式与梯形公式 z1 = 0.9602 z2 = 1.0388 z11 = 0.9602 z12 = 1.0388 z3 = 0.9995 h=pi/40; x=0:h:pi/2; y=sin(x); z1=sum(y(1:20))*h, z2=sum(y(2:21))*h, z=cumsum(y)*h, z11=z(20)*h, z12=(z(21)-z(1))*h, z3=trapz(y)*h, z= 0 0.0062 0.3802 0.0184 0.4438 0.0368 0.5107 0.0611 0.5807 0.0911 0.6533 0.1268 0.7280 0.1678 0.8043 0.2140 0.8819 0.2650 0.9602 0.3205 1.0388 ⑵ 用辛甫森公式求定积分 z4=quad('sin',0,pi/2) z5=quadl('sin',0,pi/2) z4 = 1.0000 z5 = 1.0000 ⑶ 高斯积分公式求定积分 s=gaussinteg(‘sin', 0, pi/2,1000) s = 1.0000 前面已经保存了文件gaussinteg.m在搜索路径中，即可 直接调用。这里被积函数为内部函数，无需另外定义。 一 求卫星轨道长度 我国第一颗人造地球卫星轨道为平面椭圆曲线，近地点距 地球表面439km，远地点距地面2384km，地球半径为6371km， 求卫星轨道的长度。 椭圆的参数方程为： 椭圆长度由以下椭圆积分表示： a=7782.5 km b=7721.5 km §2 数值积分应用问题举例 c a b S1 S2 解得：a=7782.5km, b=7721.5km back 我们用梯形公式和辛甫生公式来求卫星轨道的长度。 t=0:pi/10:pi/2; y1=weixing(t); format long e l1=4*trapz(t,y1) l2=4*quad