天津职业技术师范大学文档格式.docx

天津职业技术师范大学文档格式.docx

《天津职业技术师范大学文档格式.docx》由会员分享，可在线阅读，更多相关《天津职业技术师范大学文档格式.docx（13页珍藏版）》请在冰豆网上搜索。

指导教师签字：

教研室主任签字：

一、问题叙述

用数值积分法计算正弦积分函数和余弦积分函数

提示：

正弦积分，余弦函数

要求：

（1）编写函数，对任意给定的x，可求出值。

（2）使用尽可能少的正余弦函数的调用，计算更精确的值。

（用多种方法，创新方法）

2、问题分析

1、数值积分基本原理：

用数值分析求解积分的数值方法有很多，如简单的梯形法、矩形法、辛普森（Simpson）法、牛顿-科斯特（Newton-Cotes）法等都是常用的方法。

它们的基本思想都是将整个积分区间[a，b]分成n个子区间[xi,xi+1]，i=1，2，…，n，其中x1=a，xn+1=b。

这样求定积分问题就分解为求和问题。

2、本题要求用数值积分法计算正弦积分函数和余弦函数积分，那么应该从编写函数的入手，建立function的m文件，通过对函数的调用，对任意跟定的x值，求出积分函数的值。

3、数值积分法求解问题

1、梯形公式、矩形公式

首先，积分中值定理告诉我们，在积分区间[a，b]内存在一点，成立dx=（b-a）f（），就是说，底为b-a而高为f（）的矩形面积恰等于所求区边梯形的面积。

如果我们用两端点“高度”f（a）与f（b）的算术平均值作为平均高度f（）的近似值，这样导出的求积公式dx≈[f（a）+f（b）]便是我们熟悉的梯形公式。

将积分区间[a,b]n等分，每个小区间宽度均为h=，h称为积布步长。

记a=x0＜x1＜…＜xk…＜xo=b，在小区间上用小矩形面积近似小曲边梯形的面积,若分别取左端点和右端点的函数值为小矩形的高,则分别得到两个曲边梯形面积的近似计算公式。

具体程序如下：

clear

x=linspace（0,pi）;

dx=x

（2）;

y=sin（x）;

s1=sum（y）*dx

s2=trapz（y）*dx

sc1=cumsum（y）*dx;

sc2=cumtrapz（y）*dx;

plot（x,-cos（x）+1,x,sc1,'

.'

x,sc2,'

o'

）

holdon

由图可知这种方法精度太低，应选择其他方法。

2、quad函数、quan1函数

正弦：

functiony=si（t）

a=1e-8;

%函数在0点无界，去掉0点

y=quad（'

sin（x）./x'

a,t）

y=quadl（'

余弦：

functiony=ci（t）

a=-1e1;

cos（x）./x'

图像:

x=1:

100;

fori=1:

100

y2（i）=si（x（i））;

end

plot（x,y2,'

r'

）

title（'

辛普森'

y2（i）=ci（x（i））;

b'

给定任意x值，均可计算出对应的正弦、余弦函数积分。

但从结果可以看出精度不是很高。

3、复合求积公式

由于牛顿-科特斯公式在n≥8时不具有稳定性，故不可能通过提高阶的方法来提高求积精度。

为了提高精度通常可把积分区间分成若干子区间（通常是等分），再在每个子区间上用低级求积公式。

这种方法为复合求积法。

3.3.1复合梯形公式

将区间划分为n等分，分点在每个子区间上采用梯形公式，则得

记

，

称为复合梯形公式。

复合梯形公式的余项

由于且

所以使

于是复合梯形公式的余项为

事实上只要设，则可得收敛性，只要把改写成为

程序如下：

functionT_n=fhtxs（a,b,n）

h=（b-a）/n;

fork=0:

n

x（k+1）=a+k*h;

ifx（k+1）==0

x（k+1）=10^（-10）;

end

T_1=h/2*（SS（x

（1））+SS（x（n+1）））;

fori=2:

F（i）=h*SS（x（i））;

T_2=sum（F）;

T_n=T_1+T_2;

functionT_n=fhtxc（a,b,n）

T_1=h/2*（CC（x

（1））+CC（x（n+1）））;

F（i）=h*CC（x（i））;

图像：

正弦余弦

3.3.2复合新普斯求积公式

将区间划分为n等分，在每个子区间上采用辛普森公式，若记则得

●

称为复合辛普森求积公式。

正弦

functionS_n=fhxpss（a,b,n）

x_k（k+1）=x（k+1）+1/2*h;

if（x（k+1）==0）||（x_k（k+1）==0）

x_k（k+1）=10^（-10）;

S_1=h/6*（SS（x

（1））+SS（x（n+1）））;

F_1（i）=h/3*SS（x（i））;

forj=1:

F_2（j）=2*h/3*SS（x_k（j））;

S_2=sum（F_1）+sum（F_2）;

S_n=S_1+S_2;

functionS_n=fhxpsc（a,b,n）

S_1=h/6*（CC（x

（1））+CC（x（n+1）））;

F_1（i）=h/3*CC（x（i））;

F_2（j）=2*h/3*CC（x_k（j））;

图像与复合梯形所得图像基本相同，深入分析两只复合函数的优劣，

对于积分函数假设x=1，则将区间[0,1]划分为8等份，应用复合梯形求得

T8=0.9456909

而如果将[0，1]分为4等份，应用复合辛普森有

S4=0.9460832

通过参考数值分析（李庆阳）的结论，发现无论是复合梯形公式还是复合辛普森公式，最终结果都会随着h值的减小而更加精确。

对复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合梯形法的结果T8只有两位有效数字，而复合辛普森的结果却有六位有效数字，所以复合辛普森公式计算出的结果更加的精确。

4、插值型的求积公式

clc,clear

x0=0:

0.5:

5;

y0=[Inf1.75520.54030.0472-0.2081-0.3205-0.3300-0.2676-0.1634-0.04680.0567

];

%所求积分函数的数值

pp=csape（x0,y0）;

%默认的边界条件，Lagrange边界条件

formatlongg

chazhi=pp.coefs%显示每个区间上三次多项式的系数

s=quadl（@（t）ppval（pp,t）,0,5）%求积分

format%恢复短小数的显示格式

x=0:

0.1:

y=cos（x）/x;

y1=spline（x0,y0,x）;

z=0*x;

holdon

plot（x,z,x,y,'

k--'

x,y1,'

plot（x0,y0,'

*'

holdoff

clear

%所求积分函数的数值

%默认的边界条件，Lagrange边界条件

chazhi=pp.coefs%显示每个区间上三次多项式的系数

s=quadl（@（t）ppval（pp,t）,0,5）%求积分

format%恢复短小数的显示格式

如图所示：

5、高斯求积公式

function[ql,Ak,xk]=gsqj（fun,a,b,n,tol）

ifnargin==1

a=-1;

b=1;

n=7;

tol=1e-8;

elseifnargin==3

n=7;

elseifnargin==4

tol=1e-8;

elseifnargin==2||nargin>

5

error（'

TheNumberofInputArgumentsIsWrong!

'

）;

%计算求积节点

symsx

p=sym2poly（diff（（x^2-1）^（n+1）,n+1））/（2^n*factorial（n））;

tk=roots（p）;

%求积节点

%计算求积系数

Ak=zeros（n+1,1）;

n+1

xkt=tk;

xkt（i）=[];

pn=poly（xkt）;

fp=@（x）polyval（pn,x）/polyval（pn,tk（i））;

Ak（i）=quadl（fp,-1,1,tol）;

%求积系数

%积分变量代换，将[a,b]变换到[-1,1]

xk=（b-a）/2*tk+（b+a）/2;

%检验积分函数fun有效性

fun=fcnchk（fun,'

vectorize'

%计算变量代换之后积分函数的值

SS=fun（xk）*（b-a）/2;

%计算积分值

ql=sum（Ak.*SS）;

计划表

7月3号

熟悉问题，准备工作，借阅相应的资料，搞清楚题目的用意

题目要求多种方法计算，并尽量减少函数的调用。

7.4

归纳总结多种数值求积的方法，找到各种方法对应的matlab程序。

梯形

辛普森公式

复化辛普森、复化梯形④高斯勒让