数值分析法求正弦余弦积分函数.docx
《数值分析法求正弦余弦积分函数.docx》由会员分享,可在线阅读,更多相关《数值分析法求正弦余弦积分函数.docx(16页珍藏版)》请在冰豆网上搜索。
数值分析法求正弦余弦积分函数
天津职业技术师范大学
课程设计 任 务 书
理 学院 数学1403 班 学生张群
课程设计课题:
用数值积分法计算正弦积分函数和余弦积分函数
一、课程设计工作日自 2016 年 7 月4 日至 2016年 7月 5日
二、同组学生:
无
三、课程设计任务要求(包括课题来源、类型、目的和意义、基本要求、完成时间、主要参考资料等):
课题来源:
教师自拟
类型:
理论研究
目的和意义:
培养学生对数值分析中主要算法的应用能力,探索相关算法之间的内在联系。
基本要求:
根据数值分析课程所学的知识,应用Matlab软件编写程序,完成对算法及其内在原理的实验研究。
完成时间:
参考资料:
《数值分析》李庆扬等清华大学出版社
指导教师签字:
教研室主任签字:
一、问题叙述
用数值积分法计算正弦积分函数和余弦积分函数
提示:
正弦积分,余弦
函数
要求:
(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函数
正弦:
function y=si(t)
a=1e—8; %函数在0点无界,去掉0点
y=quad('sin(x)./x’,a,t)
y=quadl('sin(x)./x',a,t)
余弦:
function y=ci(t)
a=-1e1; %函数在0点无界,去掉0点
y=quad('cos(x)./x',a,t)
y=quadl('cos(x)。
/x',a,t)
图像:
x=1:
100;
fori=1:
100
y2(i)=si(x(i));
end
plot(x,y2,'r’)
title(’辛普森')
x=1:
100;
for i=1:
100
y2(i)=ci(x(i));
end
plot(x,y2,’b’)
title('辛普森')
给定任意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
end
T_1=h/2*(SS(x
(1))+SS(x(n+1)));
fori=2:
n
F(i)=h*SS(x(i));
end
T_2=sum(F);
T_n=T_1+T_2;
余弦:
functionT_n=fhtxc(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
end
T_1=h/2*(CC(x
(1))+CC(x(n+1)));
fori=2:
n
F(i)=h*CC(x(i));
end
T_2=sum(F);
T_n=T_1+T_2;
图像:
正弦 余弦
3.3.2复合新普斯求积公式
将区间
划分为n等分,在每个子区间
上采用辛普森公式,若记
则得
●
称为复合辛普森求积公式。
程序如下:
正弦
functionS_n=fhxpss(a,b,n)
h=(b-a)/n;
fork=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)));
fori=2:
n
F_1(i)=h/3*SS(x(i));
end
forj=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;
余弦:
functionS_n=fhxpsc(a,b,n)
h=(b-a)/n;
fork=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)));
fori=2:
n
F_1(i)=h/3*CC(x(i));
end
forj=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;
图像与复合梯形所得图像基本相同,深入分析两只复合函数的优劣,
对于积分函数
假设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=[ Inf 1.7552 0。
5403 0。
0472-0.2081—0。
3205 -0。
3300 -0。
2676 -0。
1634 -0.0468 0。
0567
];%所求积分函数的数值
pp=csape(x0,y0) ;%默认的边界条件,Lagrange边界条件
format longg
chazhi=pp。
coefs %显示每个区间上三次多项式的系数
s=quadl(@(t)ppval(pp,t),0,5) %求积分
format %恢复短小数的显示格式
x=0:
0。
1:
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,’*’)
holdoff
clear
x0=0:
0.5:
5;
y0=[Inf 1.7552 0.5403 0.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:
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,'*')
holdoff
如图所示:
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;tol=1e-8;
elseifnargin==4
tol=1e-8;
elseifnargin==2||nargin>5
error('TheNumberofInputArgumentsIsWrong!
');
end
%计算求积节点
symsx
p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n));
tk=roots(p);%求积节点
%计算求积系数
Ak=zeros(n+1,1);
for i=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);%求积系数
end
% 积分变量代换,将[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程序。
梯形
②辛普森公式
③复化辛普森、复化梯形④高斯勒让德
⑤插值
7。
5
编写程序,对所找的方法逐一处理编程,思考跟好的方法.
运行编程结果,进行检查改进
7。
6
编写报告总结,对程序进行系统性总结,完成课程设计
7。
7
修改论文,修改程序,检查修正出现的错误