数值分析法求正弦余弦积分函数.docx

文档介绍：该【数值分析法求正弦余弦积分函数 】是由【森林书屋】上传分享，文档一共【12】页，该文档可以免费在线阅读，需要了解更多关于【数值分析法求正弦余弦积分函数 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。天津职业技术师范大学
课 程 设 计 任 务 书
理 学院 数学1403 班 学生张群
课程设计课题:
用数值积分法计算正弦积分函数和余弦积分函数
一、课程设计工作日自2016
年7
月4
日至2016
年7
月5日
二、同组学生: 无
三、课程设计任务要求(包括课题来源、类型、目的和意义、基本要求、完成时
间、主要参考资料等):
课题来源:教师自拟
类型:理论研究
目的和意义:培养学生对数值分析中主要算法的应用能力,探索相关算法之间的内在联系。
基本要求:根据数值分析课程所学的知识,应用Matlab软件编写程序,完成对算法及其内在原理的实验研究。
完成时间:
参考资料:《数值分析》 李庆扬等 清华大学出版社
指导教师签字: 教研室主任签字:
—
一、叙述
用数值积分法计算正弦积分函数和余弦积分函数
提示:正弦分,余弦si(x)
xsintdt函数ci(x)
xcostdt
0
t
t
要求:(1)写函数,任意定的 x,可求出。
2)使用尽可能少的正余弦函数的用,算更精确的。(用多种方法,新方法)
二、分析
、数分基本原理:用数分析求解分的数方法有很多,
如的梯形法、矩形法、辛普森(Simpson)法、牛-科斯特(Newton-Cotes)法等都是常用的方法。它的基本思想都是将整个分区[a,b]分成n个子区[xi,xi+1],i=1,2,⋯,n,其中x1=a,
xn+1=b。求定分就分解求和。
2、本要求用数分法算正弦分函数和余弦函数分,那么
从写函数的入手,建立 function 的m文件,通函数的
用,任意跟定的 x,求出分函数的。
三、数分法求解
1、 梯形公式、矩形公式
首先,分中定理告我,在分区 [a,b]内存在一点 ,成
立 bf(x)dx=(b-a)f( ),就是,底 b-a而高f( )的矩
a
形面恰等于所求区梯形的面。 如果我用两端点“高度”f(a)
与f(b)的算平均作平均高度 f( )的近似,出的
b
求公式
a

f(x)dx≈b-a[f(a)+f(b)]便是我熟悉的梯形公式。
2
欢迎下载 2
—
将分区[a,b]n等分,每个小区度均
(b-a)
h=
,h称布
n
步。 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




1




0
0 1 2 3
由可知种方法精度太低,其他方法。
欢迎下载 3
—
2、quad函数、quan1函数
正弦:functiony=si(t)
a=1e-8; %函数在0点无界,去掉0点
y=quad('sin(x)./x',a,t)
y=quadl('sin(x)./x',a,t)
余弦:functiony=ci(t)
a=-1e1; %函数在0点无界,去掉0点
y=quad('cos(x)./x',a,t)
y=quadl('cos(x)./x',a,t)
图像:
x=1:100;
for i=1:100
y2(i)=si(x(i));
end
plot(x,y2, 'r' )
title( '辛普森')
辛普森









1

0 10 20 30 40 50 60 70 80 90 100
x=1:100;
for i=1:100
欢迎下载 4
—
y2(i)=ci(x(i));
end
plot(x,y2, 'b' )
title( '辛普森')
辛普森
1400
1200
1000
800
600
400
200
0
-200
-400
0 10 20 30 40 50 60 70 80 90 100
给定任意x值,均可计算出对应的正弦、余弦函数积分。但从结
果可以看出精度不是很高。
3、复合求积公式
由于牛顿-科特斯公式在 n≥8时不具有稳定性,故不可能通过提
高阶的方法来提高求积精度。为了提高精度通常可把积分区间分成若
干子区间(通常是等分),再在每个子区间上用低级求积公式。这种
方法为复合求积法。
复合梯形公式
欢迎下载 5
—
将区间a,b
划分为n等分,分点xk
akh,h
b
a,k0,1,,n,在每个子
n
区间xk,xk1
k0,1,
,n1,上采用梯形公式,则得
b
n1
xk1
h
n
1
If(x)dx
xk
f(x)dx
f(xk)f(xk1)
Rn(f)
a
k0
2k
0
记
h
Tn
2

n1
hfa
n1
fxb,
fxk
fxk1
2fxk
k0
2
k1
称为复合梯形公式。
复合梯形公式的余项
n
1
h
3
RnfITn
f''
k
,k
xkxk1
k
0
12
由于f(x)
C2a,b,
且
0
minf''
k
1n
1
f''
k
max,
kn1
nk
0
0
kn1
所以
a,b使
f''
1n
1
f''
k
nk
0
于是复合梯形公式的余项为
Rn
f
b
ah2f''
12
事实上只要设
fx
Ca,b
,则可得收敛性,只要把Tn改写成为
1
b
an1
b
a
n
Tn
n
fxk
n
fxk
2
k0
k1
程序如下:
正弦:
function T_n=fhtxs(a,b,n)
h=(b-a)/n;
for k=0:n
欢迎下载 6
—
x(k+1)=a+k*h;
if x(k+1)==0
x(k+1)=10^(-10);
end
end
T_1=h/2*(SS(x(1))+SS(x(n+1)));
for i=2:n
F(i)=h*SS(x(i));
end
T_2=sum(F);
T_n=T_1+T_2;
余弦:
function T_n=fhtxc(a,b,n)
h=(b-a)/n;
for k=0:n
x(k+1)=a+k*h;
if x(k+1)==0
x(k+1)=10^(-10);
end
end
T_1=h/2*(CC(x(1))+CC(x(n+1)));
for i=2:n
F(i)=h*CC(x(i));
end
T_2=sum(F);
T_n=T_1+T_2;
图像:
正弦
余弦
复化梯形
10
x10
2
10

8

6

4

2
1
0

10
20
30
40
50
-2
20
40
60
80
100
0
0
欢迎下载 7
—
复合新普斯求积公式
将区间[a,b]划分为n
等分,在每个子区间xk,xk1
上采用辛普森公式,若记
xk12
1
xkh,则得
2
b
n1
I
f(x)dx
f(x)dx
a
k0
hn1
6k
[f(xk)4f(xk12)f(xk)]Rn(f).
0
称为复合辛普森求积公式。
程序如下:
正弦
function S_n=fhxpss(a,b,n)
h=(b-a)/n;
for k=0:n
x(k+1)=a+k*h;
x_k(k+1)=x(k+1)+1/2*h;
if (x(k+1)==0)||(x_k(k+1)==0)
x(k+1)=10^(-10);
x_k(k+1)=10^(-10);
end
end
S_1=h/6*(SS(x(1))+SS(x(n+1)));
for i=2:n
F_1(i)=h/3*SS(x(i));
end
for j=1:n
F_2(j)=2*h/3*SS(x_k(j));
end
S_2=sum(F_1)+sum(F_2);
S_n=S_1+S_2;
余弦:
function S_n=fhxpsc(a,b,n)
h=(b-a)/n;
欢迎下载 8
—
for k=0:n
x(k+1)=a+k*h;
x_k(k+1)=x(k+1)+1/2*h;
if (x(k+1)==0)||(x_k(k+1)==0)
x(k+1)=10^(-10);
x_k(k+1)=10^(-10);
end
end
S_1=h/6*(CC(x(1))+CC(x(n+1)));
for i=2:n
F_1(i)=h/3*CC(x(i));
end
for j=1:n
F_2(j)=2*h/3*CC(x_k(j));
end
S_2=sum(F_1)+sum(F_2);
S_n=S_1+S_2;
图像与复合梯形所得图像基本相同,深入分析两只复合函数的优劣,
对于积分函数 si(x) xsintdt 假设x=1,则将区间[0,1]划分为8等份,
t
应用复合梯形求得
T8=
而如果将[0,1]分为4等份,应用复合辛普森有
S4=
通过参考数值分析(李庆阳)的结论,发现无论是复合梯形公式还是复合辛普森
公式,最终结果都会随着 h值的减小而更加精确。对复合梯形公式和复合辛普森
公式计算出的结果进行比较,发现复合梯形法的结果 T8只有两位有效数字,而
复合辛普森的结果却有六位有效数字, 所以复合辛普森公式计算出的结果更加的
精确。
4、插值型的求积公式
clc,clear
x0=0::5;
y0=[ Inf - - -
- - -
欢迎下载 9
—
]; %所求积分函数的数值
pp=csape(x0,y0); %默认的边界条件, Lagrange 边界条件
format long g
chazhi= %显示每个区间上三次多项式的系数
s=quadl(@(t)ppval(pp,t),0,5) %求积分
format %恢复短小数的显示格式
x=0::5;
y=cos(x)/x;
y1=spline(x0,y0,x);
z=0*x;
hold on
plot(x,z,x,y, 'k--' ,x,y1, 'r' )
plot(x0,y0, '*' )
hold off
clear
x0=0::5;
y0=[ Inf - - - -
- -
]; %所求积分函数的数值
pp=csape(x0,y0); %默认的边界条件, Lagrange边界条件
formatlongg
chazhi= %显示每个区间上三次多项式的系数
s=quadl(@(t)ppval(pp,t),0,5) %求积分
format %恢复短小数的显示格式
x=0::5;
y=cos(x)/x;
y1=spline(x0,y0,x);
z=0*x;
holdon
plot(x,z,x,y,'k--',x,y1,'r')
plot(x0,y0,'*')
holdoff
如图所示:
欢迎下载 10