数值分析课程设计报告积分方法的实际应用

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

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

3、的计算对象是微积分，线性代数，常微分方程中的数学问题。容包括：插值和拟合、数值微分和数值积分、求解线性方程组的直接法和迭代法、计算矩阵特征值和特征向量和常微分方程数值解等问题。2Gauss求积公式2.1根本原理求积公式(2.1)含有2n+2个待定参数， 当为等距节点时得到的插值求积公其代数精度至少为n次，如果适中选取，有可能使求积公式(2.1)具有2n+1次代数精度，这类求积公式称为高斯求积公式。为具有一般性，研究带权积分，这里为权函数，求积公式为(2.2)为不依赖于的求积系数，为求积节点，可适中选取及，使(2.2)具有2n+1次代数精度。如果求积公式(2.2)具有2n+1次代数精度，则称其节

4、点为高斯点，相应求积公式(2.2)称为高斯求积公式。根据定义要使(2.2)式具有2n+1次代数精度，只要对，()，令(2.2)式准确成立，即.(2.3)当给定权函数，求出右端积分，则可由(2.3)式解得及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);end 2.3 实例分析例计算积分解建立gaussf.m文件以调用gaussl.m文件中的函数，

5、再写入如下容：function y=gaussf(*)y=sqrt(*)*log(*);再在命令行中输入： s=gaussl(0,1,20)得出如下结果：s = -0.44563高斯-勒让德求积公式3.1根本原理在高斯求积公式(2.1)中，假设取权函数，区间为，则得公式(3.1)由于勒让德多项式是区间上的正交多项式，因此，勒让德多项式的零点就是求积公式(3.1)的高斯点。形如(3.1)式的高斯公式称为高斯-勒让德求积公式。3.2程序实现建立guasslegendre.m文件，写入如下容：function ql,Ak,*k=guasslegendre(fun,a,b,n,tol)if nargi

6、n=1 a=-1;b=1;n=7;tol=1e-8;elseif nargin=3 n=7;tol=1e-8;elseif nargin=4 tol=1e-8;elseif nargin=2|nargin5 error(The Number of Input Arguments Is Wrong!);end% 计算求积节点syms *p=sym2poly(diff(*2-1)(n+1),n+1)/(2n*factorial(n);tk=roots(p); % 求积节点% 计算求积系数Ak=zeros(n+1,1);for i=1:n+1 *kt=tk; *kt(i)=; pn=poly(*kt

7、); fp=(*)polyval(pn,*)/polyval(pn,tk(i); Ak(i)=quadl(fp,-1,1,tol); % 求积系数end% 积分变量代换，将a,b变换到-1,1*k=(b-a)/2*tk+(b+a)/2;% 检验积分函数fun的有效性fun=fchk(fun,vectorize);% 计算变量代换之后积分函数的值f*=fun(*k)*(b-a)/2;% 计算积分值ql=sum(Ak.*f*);参数说明：fun：积分表达式，可以是函数句柄a,b：积分上下限n：积分阶数tol：积分精度，默认1e-6ql：积分结果Ak:积分系数*k：求积节点，满足ql=sum(Ak.

8、*fun(*k)3.3实例分析例用4点的高斯-勒让德公式求解定积分的近似值。解：翻开guasslegendre.m文件，并在命令行中输入如下容 syms *; fun=inline(cos(*)*2); ql,Ak,*k=guasslegendre(fun,0,pi/2,4)得出结果：ql = 0.4674Ak = 0.5689 0.2369 0.4786 0.2369 0.4786*k = 0.7854 0.0737 0.3625 1.4971 1.2083即的4点的高斯-勒让德积分结果为ql=0.4674。4复化Simpson求积公式4.1根本原理 复化Simpson公式是一种比拟实用的积

9、分方法，可以给出误差估计。首先将区间a,b N等分，子区间的长度为(4.1)在每个子区间上采用Simpson公式，在用Simpson公式时，还要将子区间再二等分，因此有2N+1个分点。即(4.2)经推导得到，(4.3)称为为复化Simpson值，称4.3式为复化Simpson公式。4.2程序实现 编写复化Simpson求积函数函数名：s_quad.mFunction I=S_quad(*,y)% 复化求积公式% *为被积函数自变量的等距节点；y为被积函数在节点处的函数值。n=length(*);m=length(y); % 积分自变量的节点数应与它的函数值个数一样；if n=m error (

10、The length of * and Y must be equal); return;endif rem(n-1,2)=0 % 如果n-1不能被2整除，则调用复化公式 error (节点数不满足要求); return;endN=(n-1)/2;h=(*(n)-*(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公式法。建立一个文件SPS，容如下：clear*=input(请输入积分上

11、下限及点间的间隔例如-1:0.1:1：);y=input(请输入被积公式：y=);I=S_quad(*,y);disp(得出积分值I=)disp(I);4.3实例分析例1 用复化Simpson公式求积分，在积分区间中点与点之间的间隔取为 0.1。解：运行程序，按照提示输入积分上下限、点间的间隔及被积公式，如下所示：请输入积分上下限及点间的间隔例如-1:0.1:1：-1:0.1:1请输入被积公式：y=e*p(-*.2)得出积分值I= 1.4936真值为：1.4937例2 计算积分，将区间8等分。解：运行程序，按照提示输入积分上下限、等分后的区间长度及被积公式，如下 所示：请输入积分上下限及点间的

12、间隔例如-1:0.1:1：0:0.125:1请输入被积公式：y=*./(4+*.2)得出积分值I= 0.1116真值为：0.1115724.4结果分析复化Simpson计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。5数值方法的实际应用在实际问题中，往往会遇到一些困难。有些函数找不到用初等函数表示的原函数，例如，对于积分 (5.1)而言，不存在用初等函数表示的原函数。而有些函数虽然能找到原函数，但计算过于复杂，例如，椭圆型积分 (5.2) 而有些情况下，只能知道*些点处的函数值，并没有函数的具体表达式。这些情况，使我们有必要研究积分的数值计算问题。下面我们

13、就以梯形公式为例做以说明。 所谓梯形求积公式就是用梯形面积来近似曲边梯形面积，利用梯形公式和连续增加a,b的区间数来逼近： (5.3)第j次循环在个等距节点处对采样。5.1实例分析 卫星轨道是一个椭圆，椭圆周长计算公式是，这里a是椭圆半长轴，c是地球中心与轨道中心椭圆中心的距离，记h为近地点距离，H为远地点距离，R=637km为地球半径，则我国第一颗人造卫星近地点距离h=439km，远地点距离H=2384，试求卫星轨道的周长。解: 第一步：先利用Matlab画出被积函数的图形。 输入程序如下：clearH=2384;h=439;R=6371;a=(2*R+H+h)/2c=(H-h)/2*=0:

14、0.1:pi/2;y=sqrt(1-(c/a)2*(sin(*).2);plot(*,y,-)title(梯形法则);*label(*);ylabel(y);输出结果：a = 7.782500000000000e+003c = 9.725000000000000e+002输出图形：图5.1 被积函数的图形第二步：应用数值积分梯形公式。 首先建立一个名为trapezg.m的M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%输出用15位数字表示n=n;h=(b-a)/n; *=a+(0:n)*h;f=feval(f_name3,*);I=h

15、/2*(f(1)+f(n+1);if n1 I=I+h*sum(f(2:n);endh1=(b-a)/100;*c=a+(0:100)*h1; fc=feval(f_name3,*c);plot(*c,fc,r);hold on*label(*);ylabel(y);plot(*,f)title(数值积分梯形效果图);plot(*,zeros(size(*),.)for i=1:n;plot(*(i),*(i),0,f(i),end 然后建立一个名为f_name3.m的M文件定义函数，Matlab命令如下：function y=f_name3(*)y=sqrt(1-(9.725000e+002

16、/7.782500e+003)2*(sin(*).2)-0.99; 输入命令程序： trapezg(f_name3,0,pi/2,30) 输出结果：ans = 0.00955791054630输出图形：图5.2 数值积分效果图积分结果为：0.00955791054630+0.99=0,99955791054630第三步：计算最后结果：第四步：考虑误差。clearn=1;format longfprintf(n E*tended Trapezoidal Rulen);fprintf(n n I Errorn);I2=0.00955791054630;for k=1:8n=n*2;I1=trape

17、zg(f_name3,0,pi/2,n);format longif k=1;fprintf(%3.0f %10.8f %10.8fn, n, I1, I1-I2);endpauseend计算7步输出结果：E*tended Trapezoidal Rulen I Error4 0.00956 0.00000 8 0.00956 0.0000016 0.00956 0.0000032 0.00956 -0.0000064 0.00956 -0.00000128 0.00956 0.00000256 0.00956 -0.00000初始状态图：图5.3 初始状态图计算一步结果图：图5.4计算一步结果图计算四步结果图：图5.5 计算四步结果图最后计算八步结果图：图5.6 计算八步结果图5.2结果分析 数值微积分在科学计算上有很大应用，复化梯形公式极简化了人们实际生活中的运算复杂性。不仅算法简明，几何意义明确，而且迭代结果准确科学，具有更好的收敛性和广泛的实用性，且精度高收敛速度快。此题就是利用复化积分原理，确定轨道，并加以计算，得到高精度的结果。结论通过本次课程设计，学习到了如何应用数学软件进展大量的数值求解方法，对该课程的理解进一步加强，对课堂知识的模糊点有了一个清晰的认识。在本文表达的几种方法中，高斯-勒让德求积公式可以用较少的节点数得到高精度的计算结果，是现实生活中经常