1、数 值 分 析 课 程 设 计 报 告求积公式的实际应用学 院 数 学 与 统 计 学 院专 业 信 息 与 计 算 科 学学 号姓 名指导教师成 绩教师评语：指 导 教 师 签 字 ： 2018 年 1 月 8 日1 绪论数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现。随着计算机和计算方法的飞速发展，几乎所有学科都走向定量化和精确化，从而产生了一系列计算性的学科分支，如计算物理、计算化学、计算生物学、计算地质学、计算气象学和计算材料学等，计算数学中的数值计算方法则是解决“计算”问题的桥梁和工具。我们知道，计算能力是
2、计算工具和计算方法的效率的乘积，提高计算方法的效率与提高计算机硬件的效率同样重要。科学计算已用到科学技术和社会生活的各个领域中。数值计算方法，是一种研究并解决数学问题的数值近似解方法， 是在计算机上使用的解数学问题的方法，简称计算方法。在科学研究和工程技术中都要用到各种计算方法。 例如，在航天航空、地质勘探、汽车制造、桥梁设计、天气预报和汉字字样设计中都有计算方法的踪影。计算方法既有数学类课程中理论上的抽象性和严谨性，又有实用性和亚丁实验性的技术特征，计算方法是一门理论性和实践性都很强的学科。在 70 年代，大多数学校仅在数学系的计算数学专业和计算机系开设计算方法这门课程。随着计算机技术的迅速
3、发展和普及，现在计算方法课程几乎已成为所有理工科学生的必修课程。计算方法的计算对象是微积分，线性代数，常微分方程中的数学问题。内容包括：插值和拟合、数值微分和数值积分、求解线性方程组的直接法和迭代法、计算矩阵特征值和特征向量和常微分方程数值解等问题。2 Gauss 求积公式2.1 基本原理求积公式(2.1)0()d()nbkafxAfx含有 2n+2 个待定参数， ， 当 为等距节点时得到的插值求积公其kx,1.A k代数精度至少为 n 次，如果适当选取 ，有可能使求积公式(2.1)具有(,)kxn2n+1 次代数精度，这类求积公式称为高斯求积公式。为具有一般性，研究带权积分 ，这里 为权函数
4、，求积公式为()dbaIfx()x(2.2)0nkfAf为不依赖于 的求积系数， 为求积节点，可适当选(0,1)kAn ()x(,1)kx取 及 ，使(2.2) 具有 2n+1 次代数精度。xk如果求积公式(2.2)具有 2n+1 次代数精度，则称其节点 为高斯点，(0,)kxn相应求积公式(2.2) 称为高斯求积公式。根据定义要使(2.2) 式具有 2n+1 次代数精度，只要对 ，( )，)mf,12令(2.2)式精确成立，即. (2.3)0()dnbmkaAxx0,12n当给定权函数 ，求出右端积分，则可由(2.3)式解得 及()x (,)kx(0,1kAn2.2 程序实现建立 gauss
5、l.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);end 2.3 实例分析例 计算积分 10logdx解 建立gaussf.m 文件以调用gaussl.m文件中的函数，再写入如下内容：function y=gaussf(x)y=sqrt(x)*log(x);再在命令行中输入： s=gaussl(0,1,20)得出如下结果：s =-0.44563 高斯-勒让德求积公式3.1
6、基本原理在高斯求积公式(2.1)中，若取权函数 ， 区间为 ，则得公式()1x-1，(3.1)1-0d).nkfAf由于勒让德多项式是区间 上的正交多项式，因此，勒让德多项式 的零点就， 1P()nx是求积公式(3.1) 的高斯点。形如(3.1) 式的高斯公式称为高斯-勒让德求积公式。3.2 程序实现建立 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
7、=4tol=1e-8;elseif nargin=2|nargin5error(The Number of Input Arguments Is Wrong!);end% 计算求积节点syms xp=sym2poly(diff(x2-1)(n+1),n+1)/(2n*factorial(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); % 求积系数en
8、d% 积分变量代换，将a,b变换到-1,1xk=(b-a)/2*tk+(b+a)/2;% 检验积分函数 fun 的有效性fun=fcnchk(fun,vectorize);% 计算变量代换之后积分函数的值fx=fun(xk)*(b-a)/2;% 计算积分值ql=sum(Ak.*fx);参数说明：fun：积分表达式，可以是函数句柄a,b：积分上下限n：积分阶数tol：积分精度，默认 1e-6ql：积分结果Ak:积分系数xk：求积节点，满足 ql=sum(Ak.*fun(xk)3.3 实例分析例 用 4 点的高斯-勒让德公式求解定积分 的近似值。20cosdx解：打开 guasslegendre.
9、m 文件，并在命令行中输入如下内容 syms x; fun=inline(cos(x)*x2); ql,Ak,xk=guasslegendre(fun,0,pi/2,4)得出结果：ql =0.4674Ak =0.56890.23690.47860.23690.4786xk =0.78540.07370.36251.49711.2083即 的 4 点的高斯-勒让德积分结果为 ql=0.4674。20cosdx4 复化 Simpson 求积公式4.1 基本原理复化 Simpson 公式是一种比较实用的积分方法，可以给出误差估计。首先将区间a,b N 等分，子区间的长度为(4.1)nbahN在每个子
10、区间上采用 Simpson 公式，在用 Simpson 公式时，还要将子区间再二等分，因此有 2N+1 个分点。即(4.2)0 0,1,2,.2Nkhxkxa经推导得到，(4.3)1a 221()d()()4()6NNdefbNN kkkSxfabfxfx称为 为复化 Simpson 值，称（4.3）式为复化 Simpson 公式。N4.2 程序实现编写复化 Simpson 求积函数（函数名：s_quad.m ）Function I=S_quad(x,y)% 复化求积公式% x 为被积函数自变量的等距节点；y 为被积函数在节点处的函数值。n=length(x);m=length(y); % 积
11、分自变量的节点数应与它的函数值个数相同；if n=merror (The length of 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;endI=h/6*sum(a.*y);然后调用 s_quad 函数，来实现复化 Simpson
12、 公式法。建立一个文件 SPS，内容如下：clearx=input(请输入积分上下限及点间的间隔（例如-1:0.1:1）：);y=input(请输入被积公式：y=);I=S_quad(x,y);disp(得出积分值 I=)disp(I);4.3 实例分析例 1 用复化 Simpson 公式求积分 ，在积分区间中点与点之间的间隔取为 21dxe0.1。解：运行程序，按照提示输入积分上下限、点间的间隔及被积公式，如下所示：请输入积分上下限及点间的间隔（例如-1:0.1:1）：-1:0.1:1请输入被积公式：y=exp(-x.2)得出积分值 I=1.4936真值为：1.4937例 2 计算积分 ，将
13、区间 8 等分。120xd4解：运行程序，按照提示输入积分上下限、等分后的区间长度及被积公式，如下 所示：请输入积分上下限及点间的间隔（例如-1:0.1:1）：0:0.125:1请输入被积公式：y=x./(4+x.2)得出积分值 I=0.1116真值为：0.1115724.4 结果分析复化 Simpson 计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。5 数值方法的实际应用在实际问题中，往往会遇到一些困难。有些函数找不到用初等函数表示的原函数，例如，对于积分(5.1)10sindx而言，不存在用初等函数表示的原函数。而有些函数虽然能找到原函数，但计算过于
14、复杂，例如，椭圆型积分(5.2)21xdabcx而有些情况下，只能知道某些点处的函数值，并没有函数的具体表达式。这些情况，使我们有必要研究积分的数值计算问题。下面我们就以梯形公式为例做以说明。所谓梯形求积公式就是用梯形面积来近似曲边梯形面积，利用梯形公式和连续增加a,b的区间数来逼近：(5.3)21a1()d()(jbkkkhfxfxf第 j 次循环在 个等距节点处对 采样。12j f5.1 实例分析卫星轨道是一个椭圆，椭圆周长计算公式是，220cS4a1-sind（ ）这里 a 是椭圆半长轴， c 是地球中心与轨道中心（椭圆中心）的距离，记 h 为近地点距离，H 为远地点距离，R=637km
15、 为地球半径，则,2/)(,2/)2(ahHchR我国第一颗人造卫星近地点距离 h=439km，远地点距离 H=2384，试求卫星轨道的周长。解: 第一步：先利用 Matlab 画出被积函数的图形。输入程序如下：clearH=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(梯形法则);xlabel(x);ylabel(y);输出结果：a =7.782500000000000e+003c =9.725000000000000e+002输出图形：0 0.
16、5 1 1.50.9920.9930.9940.9950.9960.9970.9980.99911.001 值值值值xy图 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:n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);if n1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1; fc=
17、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),end然后建立一个名为 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)输出结果
18、：ans =0.00955791054630输出图形：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.600.0020.0040.0060.0080.010.012 值值值值值值值值值xy图 5.2 数值积分效果图积分结果为：0.00955791054630+0.99=0,99955791054630第三步：计算最后结果：041928e4.87035pi/2)*.9563(.7*sinc-aS22d）（第四步：考虑误差。clearn=1;format longfprintf(n Extended Trapezoidal Rulen);fprintf(n n I Errorn);I2
19、=0.00955791054630;for k=1:8n=n*2;I1=trapezg(f_name3,0,pi/2,n);format longif k=1;fprintf(%3.0f %10.8f %10.8fn, n, I1, I1-I2);endpauseend计算 7 步输出结果：Extended Trapezoidal Rulen I Error4 0.00956 0.000008 0.00956 0.0000016 0.00956 0.0000032 0.00956 -0.0000064 0.00956 -0.00000128 0.00956 0.00000256 0.00956
20、 -0.00000初始状态图：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.600.0020.0040.0060.0080.010.012 值值值值值值值值值xy图 5.3 初始状态图计算一步结果图：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.600.0020.0040.0060.0080.010.012 值值值值值值值值值xy图 5.4 计算一步结果图计算四步结果图：0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.600.0020.0040.0060.0080.010.012 值值值值值值值值值xy图 5.5 计算四步结果图最后计算八步结果图：0 0
21、.2 0.4 0.6 0.8 1 1.2 1.4 1.600.0020.0040.0060.0080.010.012 值值值值值值值值值xy图 5.6 计算八步结果图5.2 结果分析数值微积分在科学计算上有很大应用，复化梯形公式极大地简化了人们实际生活中的运算复杂性。不仅算法简明，几何意义明确，而且迭代结果准确科学，具有更好的收敛性和广泛的实用性，且精度高收敛速度快。本题就是利用复化积分原理，确定轨道，并加以计算，得到高精度的结果。结论通过本次课程设计，学习到了如何应用数学软件进行大量的数值求解方法，对该课程的理解进一步加强，对课堂知识的模糊点有了一个清晰的认识。在本文叙述的几种方法中，高斯-
22、勒让德求积公式可以用较少的节点数得到高精度的计算结果，是现实生活中经常用到的数值积分方法，但是当积分区间较大时，积分精度并不理想。复化 Simpson 计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。复化梯形公式收敛且稳定，极大地简化了人们实际生活中的运算复杂性。不仅算法简明，几何意义明确，而且迭代结果准确科学，具有更好的收敛性和广泛的实用性，且精度高收敛速度快。参考文献1 李庆扬,王能超,易大义.数值分析M. 北京: 清华大学出版社 , 2008.97-137.2 邓东皋,伊小玲.数学分析简明教程M. 北京: 高等教育出版社,2006.100-136.3 薛定宇,陈阳泉.高等应用数学问题的 MATLAB求解M. 北京: 清华大学出版社,2013.35-42.
