数值分析教学规划积分方法的实际应用

1、_数值分析课程设计报告求积公式的实际应用学院数学与统计学院专业信息与计算科学学号姓名指导教师成绩教师评语：指导教师签字：_2018年1月8日_1 绪论数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现。随着计算机和计算方法的飞速发展，几乎所有学科都走向定量化和精确化，从而产生了一系列计算性的学科分支，如计算物理、计算化学、计算生物学、计算地质学、计算气象学和计算材料学等，计算数学中的数值计算方法则是解决“计算”问题的桥梁和工具。我们知道，计算能力是计算工具和计算方法的效率的乘积，提高计算方法的效率与提高计算机硬件的效率同

2、样重要。科学计算已用到科学技术和社会生活的各个领域中。数值计算方法， 是一种研究并解决数学问题的数值近似解方法，是在计算机上使用的解数学问题的方法， 简称计算方法。在科学研究和工程技术中都要用到各种计算方法。例如，在航天航空、地质勘探、汽车制造、桥梁设计、天气预报和汉字字样设计中都有计算方法的踪影。计算方法既有数学类课程中理论上的抽象性和严谨性，又有实用性和亚丁实验性的技术特征，计算方法是一门理论性和实践性都很强的学科。在70 年代，大多数学校仅在数学系的计算数学专业和计算机系开设计算方法这门课程。随着计算机技术的迅速发展和普及，现在计算方法课程几乎已成为所有理工科学生的必修课程。计算方法的计

3、算对象是微积分，线性代数，常微分方程中的数学问题。内容包括：插值和拟合、数值微分和数值积分、求解线性方程组的直接法和迭代法、计算矩阵特征值和特征向量和常微分方程数值解等问题。2Gauss 求积公式2.1基本原理求积公式bf (x)d xa含有 2n+2 个待定参数， xk ， Ak (k0,1,代数精度至少为n 次，如果适当选取 xk (k_nAk f ( xk )(2.1)k 0, n). 当 xk 为等距节点时得到的插值求积公其0,1, n) ，有可能使求积公式 (2.1) 具有 2n+1次代数精度，这类求积公式称为高斯求积公式。Ib( x) 为权函数，求积公式为为具有一般性，研究带权积分

4、f ( x) (x)d x ，这里abnf (x)( x)d xAk f ( xk )(2.2)ak 0Ak (k 0,1, , n) 为不依赖于 f ( x) 的求积系数， xk ( k0,1, n) 为求积节点，可适当选取xk 及 Ak ，使 (2.2) 具有 2n+1 次代数精度。如果求积公式 (2.2) 具有 2n+1次代数精度，则称其节点xk ( k 0,1, , n) 为高斯点，相应求积公式 (2.2) 称为高斯求积公式。根据定义要使(2.2)式具有2n+1次代数精度， 只要对 fxxm ，( m0,1,2n 1 )，( )令(2.2) 式精确成立，即nAk xkmbxm ( x)

5、 d xm 0,1, 2n1.(2.3)k 0a当 给 定 权 函 数 ( x)，求出右端积分，则可由(2.3)式 解 得 xk (k0,1,n) 及Ak (k0,1, , n)2.2程序实现建立 gaussl.m文件，写入如下内容：function s=gaussl(a,b,n)_h=(b-a)/n;s=0.0;for m=0:(1*n/2-1)s=s+h*(gaussf(a+h*(1-1/sqrt(3)+2*m)+gaussf(a+h*(1+1/sqrt(3)+2*m); end2.3实例分析1例计算积分x log xd x0解建立 gaussf.m 文件以调用 gaussl.m 文件中的

6、函数，再写入如下内容：function y=gaussf(x)y=sqrt(x)*log(x);再在命令行中输入：>> s=gaussl(0,1,20)得出如下结果：s =-0.44563 高斯 - 勒让德求积公式3.1基本原理在高斯求积公式 (2.1) 中，若取权函数( x) 1， 区间为 -1，1，则得公式1n(3.1)f ( x) d xAk f ( xk ).-1k 0_由于勒让德多项式是区间-1，1 上的正交多项式，因此，勒让德多项式Pn 1( x) 的零点就是求积公式 (3.1) 的高斯点。形如 (3.1) 式的高斯公式称为高斯 - 勒让德求积公式。3.2程序实现建立

7、guasslegendre.m文件，写入如下内容：function ql,Ak,xk=guasslegendre(fun,a,b,n,tol)if nargin=1a=-1;b=1;n=7;tol=1e-8;elseif nargin=3n=7;tol=1e-8;elseif nargin=4tol=1e-8;elseif nargin=2|nargin>5error('The Number of Input Arguments Is Wrong!');end% 计算求积节点syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factori

8、al(n);tk=roots(p); %求积节点% 计算求积系数Ak=zeros(n+1,1);for i=1:n+1xkt=tk;_xkt(i)=;pn=poly(xkt);fp=(x)polyval(pn,x)/polyval(pn,tk(i);Ak(i)=quadl(fp,-1,1,tol); %求积系数end% 积分变量代换，将 a,b 变换到 -1,1 xk=(b-a)/2*tk+(b+a)/2;% 检验积分函数 fun 的有效性fun=fcnchk(fun,'vectorize');% 计算变量代换之后积分函数的值fx=fun(xk)*(b-a)/2;% 计算积分值

9、ql=sum(Ak.*fx);参数说明：fun ：积分表达式，可以是函数句柄a,b ：积分上下限n ：积分阶数tol ：积分精度，默认1e-6ql ：积分结果Ak: 积分系数xk ：求积节点，满足ql=sum(Ak.*fun(xk)_3.3实例分析例用 4 点的高斯 - 勒让德公式求解定积分2 x2 cos x dx 的近似值。0解：打开 guasslegendre.m文件，并在命令行中输入如下内容>> syms x;>> fun=inline(cos(x)*x2);>> ql,Ak,xk=guasslegendre(fun,0,pi/2,4)得出结果：ql

10、 =0.4674Ak =0.56890.23690.47860.23690.4786xk =_0.78540.07370.36251.49711.2083即2 x2 cos x dx 的 4 点的高斯 - 勒让德积分结果为ql=0.4674 。04复化 Simpson求积公式4.1基本原理复化 Simpson公式是一种比较实用的积分方法，可以给出误差估计。首先将区间a,b N 等分，子区间的长度为b a(4.1)hnN在每个子区间上采用 Simpson公式，在用 Simpson公式时，还要将子区间再二等分，因此有 2N+1 个分点。即xkx0k hN , k 0,1,2 N , x0 a.(4

11、.2)2经推导得到，defhN f (a) f (b) 2N 1NSNab f ( x)d xf (x2k )4f ( x2k 1)(4.3)6k 1k1称为 SN 为复化 Simpson值，称（ 4.3 ）式为复化 Simpson公式。4.2程序实现_编写复化 Simpson求积函数（函数名： s_quad.m ）FunctionI=S_quad(x,y)% 复化求积公式% x 为被积函数自变量的等距节点； y 为被积函数在节点处的函数值。n=length(x);m=length(y);% 积分自变量的节点数应与它的函数值个数相同；if n=merror ('The length o

12、f X and Y must be equal');return;endif rem(n-1,2)=0% 如果 n-1 不能被 2 整除，则调用复化公式error (' 节点数不满足要求 ');return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros (1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;end_I=h/6*sum(a.*y);然后调用 s_quad 函数，来实现复化Simpson公式法。建立一个文件SPS，内容如下：clearx=inpu

13、t('请输入积分上下限及点间的间隔（例如-1:0.1:1 ）：');y=input('请输入被积公式： y=');I=S_quad(x,y);disp(' 得出积分值 I=')disp(I);4.3实例分析例 1用复化 Simpson公式求积分1e x2d x ，在积分区间中点与点之间的间隔取为10.1 。解：运行程序，按照提示输入积分上下限、点间的间隔及被积公式，如下所示：请输入积分上下限及点间的间隔（例如-1:0.1:1 ）： -1:0.1:1请输入被积公式： y=exp(-x.2)得出积分值 I=1.4936真值为： 1.4937_1 x例

14、 2 计算积分 0 4 x2 d x ，将区间 8 等分。解：运行程序，按照提示输入积分上下限、等分后的区间长度及被积公式，如下所示：请输入积分上下限及点间的间隔（例如-1:0.1:1 ）：0:0.125:1请输入被积公式： y=x./(4+x.2)得出积分值 I=0.1116真值为： 0.1115724.4结果分析复化 Simpson计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。5 数值方法的实际应用在实际问题中，往往会遇到一些困难。有些函数找不到用初等函数表示的原函数，例如，对于积分1 sin xd x(5.1)0 x而言，不存在用初等函数表示的原函

15、数。而有些函数虽然能找到原函数，但计算过于复_杂，例如，椭圆型积分x2(5.2)ax2 bx c d xx1而有些情况下， 只能知道某些点处的函数值， 并没有函数的具体表达式。 这些情况，使我们有必要研究积分的数值计算问题。下面我们就以梯形公式为例做以说明。所谓梯形求积公式就是用梯形面积来近似曲边梯形面积，利用梯形公式和连续增加a,b 的区间数来逼近：bh 2j(5.3)f ( x)d x f ( xk 1 ) f (xk )a2 k 1第 j 次循环在 2j1 个等距节点处对 f x采样。5.1实例分析卫星轨道是一个椭圆，椭圆周长计算公式是S 4a 21-（c 22d，） sin0 a这里

16、a 是椭圆半长轴， c 是地球中心与轨道中心（椭圆中心）的距离，记h 为近地点距离， H 为远地点距离， R=637km为地球半径，则a( 2RHh) / 2,c(Hh) / 2,我国第一颗人造卫星近地点距离h=439km，远地点距离 H=2384 ，试求卫星轨道的周长。解: 第一步：先利用 Matlab 画出被积函数的图形。输入程序如下：clear_H=2384;h=439;R=6371;a=(2*R+H+h)/2c=(H-h)/2x=0:0.1:pi/2;y=sqrt(1-(c/a)2*(sin(x).2);plot(x,y,'-')title(' 梯形法则 

17、9;);xlabel('x');ylabel('y');输出结果：a =7.782500000000000e+003c =9.725000000000000e+002_输出图形：梯形法则1.00110.9990.9980.997y0.9960.9950.9940.9930.99200.511.5x图 5.1 被积函数的图形第二步：应用数值积分梯形公式。首先建立一个名为trapezg.m的 M 文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%输出用 15 位数字表示n=n;h=(b-a)/n;x=a+(0:

18、n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);_if n>1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1;fc=feval(f_name3,xc);plot(xc,fc,'r');hold onxlabel('x');ylabel('y');plot(x,f)title(' 数值积分梯形效果图 ');plot(x,zeros(size(x),'.')for i=1:n;plot(x(i),x(i),0,f(i),en

19、d然后建立一个名为f_name3.m的 M 文件定义函数， Matlab命令如下：function y=f_name3(x)y=sqrt(1-(9.725000e+002/7.782500e+003)2*(sin(x).2)-0.99;输入命令程序：>> trapezg('f_name3',0,pi/2,30)输出结果：ans =_0.00955791054630输出图形：数值积分梯形效果图0.0120.010.008y 0.0060.0040.002000.20.40.60.811.21.41.6x图 5.2 数值积分效果图积分结果为：0.009557910546

20、30+0.99=0,99955791054630第三步：计算最后结果：2（c 22S 4a）sin d1-0 a4* 7782.5* (0.009557910546300.99* pi/2)4.870743851190028e 004第四步：考虑误差。_clearn=1;format longfprintf('n Extended Trapezoidalfprintf('nnIRulen');Errorn');I2=0.00955791054630;for k=1:8n=n*2;I1=trapezg('f_name3',0,pi/2,n);for

21、mat longif k=1;fprintf('%3.0f%10.8f%10.8fn',n,I1,I1-I2);endpauseend计算 7 步输出结果：ExtendedTrapezoidalRulenIError40.009560.0000080.009560.00000160.009560.00000320.00956-0.00000_640.00956-0.000001280.009560.000002560.00956-0.00000初始状态图：数值积分梯形效果图0.0120.010.008y 0.0060.0040.002000.20.40.60.811.21.41.6x图 5.3 初始状态图计算一步结果图：_数值积分梯形效果图0.0120.010.008y 0.0060.0040.002000.20.40.60.811.21.41.6x图 5.4 计算一步结果图计算四步结果图：数值积分梯形效果图0.0120.010.008y 0.0060.0040.00