数值计算课程设计

四川理工学院

《数值计算方法》课程设计

题目利用带权高斯求积公式求定积分专业信息与计算科学

班级2014级1班

姓名裴秋波覃勇潮王萌孙波

一、摘要 (3)

二、应用计算方法的基本原理 (4)

1.高斯-切比雪夫求积公式.......... 错误！未定义书签。

1.1算法描述 (4)

2.高斯-拉盖尔求积公式 (4)

2.1算法描述 (4)

三、例题的计算结果及误差分析........ 错误！未定义书签。

四、总结及心得体会 (6)

五、参考文献 (7)

六、附录程序 (7)

一、摘要

数值积分一直都是数值计算领域一个重要的分支，很多工程问题都要运用数值积分的相关知识，而高斯求积公式又是一种重要的数值积分的方法，对比与牛顿-科特斯公式，高斯公式不但是高精确度的，而且是数值稳定的，收敛的。高斯-切比雪夫求积公式又是在高斯公式的基础上发展而来的，是高斯公式的补充与完善。本文就是重点介绍基于MATLAB软件下高斯求积公式的求解，通过MATLAB数学软件我们不仅能够比较精确的解决高斯公式的求积问题，还可以比较形象易懂地反映出随着高斯点的增加误差的变化。

关键字：数值积分高斯公式稳定精确度 MATLAB

二、应用计算方法的基本原理

1.高斯-切比雪夫求积公式

1.1算法描述

在区间[-1,1]上取权函数()2

11x

x -=ρ 的正交多项式是Chebyshev 正交多项

式。n+1次Chebyshev 多项式

()()[]x n x T n arccos 1cos 1

+=+

的零点为n k n k x k ,,1,0,2

21

2cos

=++=π 以此为Gauss 点,利用Chebyshev 多项式的性质可得相应的求积系数为

()n k n dx x x

l A k

k ,,1,0,1

111

1

2 =+=-=?

-π

其中()x l k 是关于Gauss 点的Lagrange 插值基函数.从而有Gauss-

Chebyshev 求积公式如下

()(),1

110

1

1

2

∑?

=-+≈

-n

k k

x f n dx x f x

π

对于n=1, 二次Chebyshev 多项式为12-x x ，二点Gauss-Chebyshev 求积公式为

()???

?????? ??+??? ?

?

-≈

-?

-21212111

1

2

f f dx x f x π 对于n=2,三次Chebyshev 多项式为x x 343-，三点Gauss-Chebyshev 求积公式为

()()????

???????? ??++???

? ??-≈

-?

-230233111

1

2

f f f dx x f x π

2.高斯-拉盖尔求积公式

2.1算法描述

将插值型求积公式中的区间[]b a ,换成区

间[]∞,0，权函数取为()x e x -=ρ，取节点()n k x k ,,1,0 =为n+1次拉盖尔多项式

()()x

n n n x n e x dx d e x L -++++=1111

的零点，称这样的高斯型求积公式为高斯-拉盖尔 求积公式，其表示式为

()()

n

x

k

k k e

f

x dx A f

x +∞

-=≈∑?

其中

2

21[(1)!](0,1,,)

[()]k k n

k n A k n x L x ++=

='

截断误差为

三、例题的计算结果及误差分析

1、高斯-切比雪夫求积公式

1.1例 用三点和四点Gauss —Chebyshev 求积分计算公式?--=1

1

2

1I x

e x dx,并

估计误差。

解：这里x e x f =)(，()e x =x ，由Gauss —Chebyshev 求积公式（5.4.9）可得

?

∑-=+≈

-=

1

1

2

1

1I n

k x x k e n dx x e π

当n=2时，(),2,1,0612cos k =??

?

??+=k k x π求得

2365cos ,063cos ,2

36cos 210-

=====???

??=πππx x x 代入上式得

()97732.3042062.137744.2332

30

2

3

≈++=???

?

??++≈

π

πe

e e I 估计误差可用余项表达式（2.8），因()()()x x e x

f e x f ==6,，故

()()()45652107065.3!

62!

62|

||-?≤?≤

?=e f f R π

ηπ

，

当n=3时， ().3,2,1,0812cos =??

? ??+=k k x k π，求得

2

(22)

[(1)!][](),(0,)

(22)!

n n R f f

n ζζ++=∈∞+

9238795

.0,3826834.03826834

.0.,9238795.04210-=-===x x x x ()∑=+++=

≈

3

3969760

.06820288.04662138.15190441.24

4

k x k e I π

π

误差 ()()()6787310655.1!

82|!

82|||-?≤?≤

?=e f f R π

ηπ

2、高斯-拉盖尔求积公式

2.1例 用高斯型积分公式计算积分：20

1

1x

dx e ∞

+?

四、总结及心得体会

在做课程设计前,我以为不会难做,就像以前做物理实验一样,做完实验,然后两下子就将实验报告做完.直到做完测试实验时,我才知道其实并不容易做,但学到的知识与难度成正比,使我受益匪浅。

本次试验的目的是基于MATLAB 软件下高斯求积公式的求解，通过MATLAB 数学软件我们能够比较精确的解决高斯公式的求积问题，还可以比较形象易懂地反映出随着高斯点的增加误差的变化。

首先，，做实验之前一定要把书本上的只是吃透，这是做实验的基础，否则

做实验时的难度会很大，做实验时,一定要亲力亲为,务必要将每个细节弄清楚,弄明白,实验后,还要复习,思考,这样,你的印象才深刻,记得才牢固,否则,过后不久你就会忘得一干二净,这还不如不做.这样能使我们认识到这门课程在数学应用是那么的广泛.

通过这次试验，不仅让我们更加清楚高斯-切比雪夫求积公式的原理以及应用，更是加强了我们的MATLAB软件的应用，这让我们受益良多，此次的设计是1个新的体验、新的挑战，我在此过程中不但应用了所学的知识，而且还不断的学习新的知识、工具，以完成设计的需要，

此次设计，我们以下的能力得到提升：

1 我们对高斯-切比雪夫求积公式的理解与应用

2 对MATLAB软件的应用得到加强

3 由于是由一个小组合作完成，这让我们意识到团队合作的重要性，加强了我们小组之间的互动与合作能力。

在这次课程设计中，我们用到了以前所学的专业知识，虽然过去这些知识用于解题中，未有独立应用过它们，但在学习的过程中带着问题去学返现效率很高，设计过程，好比我们的成长过程，常有一些不如意，难免会遇到各种各样的问题，这也激发了我们今后努力学习的兴趣，通过这次设计，我们懂得了学习的重要性，了解到理论知识与实践结合的重要意义，学会了坚持丶耐心和努力，这将为今后的学习和工作打下基础。

由于这是1个新的尝试，还有我们有限的知识和仓促的时间所限，所以设计过程中难免有缺点和不足的地方，望各位老师能给予我批评指正。

五、参考文献

[1]李庆扬,王能超,易大义. 数值分析(第5版). 北京：清华大学出版社，2005

[2]刘琼荪编.数学实验[M].北京:高等教育出版社,20l4年7月

[3]王建卫曲中水凌演编著.MATLAB 7.X程序设计.中国水利水电出版社,2007年9月

[4]谢东秀左军编著.数值计算方法与实验：国防工业出版社.2014年8月

六、附录程序

高斯-切比雪夫求积公式：

新建M文件guasschebyshev.m

function [gc,Ak,xk]=guasschebyshev(fun,n)

% 高斯-勒让德数值积分

%

% 参数说明

% fun：积分表达式，可以是函数句柄、inline函数、匿名函数、字符串表达式，但是必须可以接受矢量输入

% a,b：积分上下限只能[-1 1]，故忽略

% n：积分阶数，可以任意正整数，但是不建议设置过大，大不一定能得到更好的精度

% gc：积分结果

% Ak：求积系数

% xk：求积节点，满足gc=sum(Ak.*fun(xk))

%

% 举例说明

% fun=@(x)exp(x).*cos(x); % 必须可以接受矢量输入

% quadl(fun,-1,1) % 调用MATLAB内部积分函数检验

% [gc,Ak,xk]=guasschebyshev(fun,7)

%

% 高斯积分说明

% 1.高斯积分是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。一般插值型数值积分精度为至少n阶，且具有高阶不稳定性

% 2.高斯求积节点就是对应n阶正交多项式的根，构建首项系数为1的正交多项式参见https://www.360docs.net/doc/c17802852.html,/thread-4798-1-1.html

% 3.高斯求积系数，可以由Lagrange多项式插值系数进行积分得到

% 4.由高斯求积节点为根构成的多项式，与任何不超过n阶的多项式正交

%

% 勒让德正交多项式说明

% 1.区间[-1,1]上关于权rho(x)=1/sqrt(1-x^2)的正交多项式系，满足

%

|- pi/2 (i=j≠0)

% ∫(Ti(x)*Tj(x)/sqrt(1-x^2),-1,1)= | pi (i=j=0) %

|- 0 (i≠j)

% 2.切比雪夫正交多项式的通式为：T0=1, Tn=cos(n*arcos(x)) (n=1,2,...) % 3.Tn(x)=0在[-1,1]有n个零点xk=cos((2k-1)/(2n)*pi) (k=1,2,...n)，

% 有n+1个最值点pk=cos(k*pi/n) (k=0,2,...,n)

% 3.关于高斯-切比雪夫的数值积分的求积系数Ak=pi/(n+1)，求积节点xk=cos((2k+1)/(2n+2)*pi) (k=0,2,...,n)

%

% by dynamic of Matlab技术论坛

% see also https://www.360docs.net/doc/c17802852.html,

% contact me matlabsky@https://www.360docs.net/doc/c17802852.html, % 2010-01-15 23:16:33

%

% 计算求积系数

k=0:n;

Ak=ones(1,n+1)*pi/(n+1);

% 计算求积节点

xk=cos((2*k+1)/(2*n+2)*pi);

% 计算权函数

rho=@(x)sqrt(1-x.^2);

% 积分函数有效性检验

fun=fcnchk(fun,'vectorize');

% 带权的积分函数值

fx=fun(xk).*rho(xk);

% 计算积分值

gc=sum(Ak.*fx);

在MATLAB命令窗口输入

n=10；

k=0:n;

Ak=ones(1,n+1)*pi/(n+1)

% 计算求积节点

xk=cos((2*k+1)/(2*n+2)*pi)

% 计算权函数

rho=@(x)sqrt(1-x.^2)；

% 积分函数有效性检验

fun=fcnchk(fun,'vectorize');

% 带权的积分函数值

fx=fun(xk).*rho(xk)

% 计算积分值

gc=sum(Ak.*fx)；

得到结果

高斯-拉盖尔求积公式：

#include

#include

double Lag(double x)

{

double z;

z=1/(1+exp(2*x));

return z;

}

void main()

{

double x[7],A[7],I=0;

int i,n;

printf("请输入点数n:");

scanf("%d",&n);

switch(n)

{

case 2: x[1]=0.5857864376,x[2]=3.4142135624;

A[1]=0.853*******,A[2]=0.1464466094;break;

case 3: x[1]=0.4157745567,x[2]=2.2942803602,x[3]=6.2899450829;

A[1]=0.7110930099,A[2]=0.2785177335,A[3]=0.010*******;break;

case 4: x[1]=0.3225476896,x[2]=1.7457611011,x[3]=4.5366202969,x[4]=9.39507091 23;

A[1]=0.6031541043,A[2]=0.3574186924,A[3]=0.0388879085,A[4]=0.00053929 47;break;

case 5: x[1]=0.2635603197,x[2]=1.4034030591,x[3]=3.5964257710,x[4]=7.08581000 58,x[5]=12.6408008442;

A[1]=0.5217556105,A[2]=0.3986668110,A[3]=0.0759424497,A[4]=0.00361175 87,A[5]=0.0000233700;break;

case 6: x[1]=0.2228466041,x[2]=1.1889321016,x[3]=2.9927363260,x[4]=5.77514356 91,x[5]=9.8374674183,x[6]=15.9828739806;

A[1]=0.4589646793,A[2]=0.4170008307,A[3]=0.1133733820,A[4]=0.01039919 75,A[5]=0.0002610172,A[6]=0.0000008985;break;

default :printf("errer\n");

}

for(i=1;i<=n;i++)

I=I+A[i]*Lag(x[i]);

printf("原积分I=%f\n",I);

}

任务分工：

裴秋波：25%（书写格式规范，拉盖尔求积公式程序代码）

覃勇潮：25%（原理描述，切比雪夫求积公式程序代码）

王萌：25%（例题及误差分析）

孙波：25%（总结）