5《数值逼近》课程设计报告.docx
《5《数值逼近》课程设计报告.docx》由会员分享，可在线阅读，更多相关《5《数值逼近》课程设计报告.docx（31页珍藏版）》请在冰豆网上搜索。
5《数值逼近》课程设计报告
数值逼近
课程设计报告
一、目的意义
（1）进一步熟悉掌握复化梯形和复化抛物线公式
（2）学会比较复化梯形公式和复化抛物线公式如何达到所要求的精度
（3）提高编程能力
（4）通过数值方法求出很难求得原函数的积分和解析表达是没有明确的给出积分的近似值
二、内容要求
积分计算问题：
分别用复化梯形和复化Simpsom求积公式计算积分
，并比较计算量（精度为10-8）。
三、问题解决的方法与算法
方法：
解决该积分问题时，运用了数值积分近似解法的方法，运用复化梯形和复化Simpsom求积公式进行计算
3.1复化梯形积分
3.1.1复化梯形积分公式表达式
3.1.2复化梯形积分误差表达式
3.2复化抛物线积分
3.2.1复化抛物线积分公式表达式
3.2.2复化抛物线积分误差表达式
3.3算法
3.3.1复化梯形积分算法
第一步：
根据精度计算n的值，输入两端点的值，计算步长h
第二步：
根据步长计算出各个节点x[i]的值，i=0,1,2,…,n
第三步：
根据x[i]计算出各个节点对应y[i]的值，i=0,1,2,…,n
第四步：
对各个节点的值进行求和
第五步：
输出最终的积分的值
3.3.2复化抛物线积分算法
第一步：
根据精度计算n的值，端点a,b的值，计算步长h
第二步：
根据步长计算出各个节点x[i]的值，i=0,1,2,…,n
第三步：
根据x[i]计算出各个节点对应y[i]的值，i=0,1,2,…,n
第四步：
对各个节点的值进行求和,分情况，对左右端点先求和，对剩下的端点，奇数的求和后乘以4倍，偶数的求和后乘以2倍，最终将各个值进行加和
第五步：
对加和的值乘以步长除以3
第六步：
输出最终的积分的值
四、计算程序
//复化梯形公式.cpp:
定义控制台应用程序的入口点。
//n+1点的复化梯形公式
#include
//#include"stdafx.h"
#include
#include
#include
usingnamespacestd;
constintARRAY_LEN（10000）;
classComt
{protected:
doublea,b,h,n,x[ARRAY_LEN],y[ARRAY_LEN];
charf[ARRAY_LEN];
public:
voidgetab（）
{
cout<<"请输入该积分的上下限（即区间）:
";
cin>>a>>b;
cout<}
voidcal_nh（）
{
intc;
cout<<"几点的复化梯形公式？
";
cin>>c;
cout<n=c-1;
cout<<"n的值为：
"<h=（b-a）/n;
cout<<"h的值为：
"<}
voidcal_x（）
{
inti=0;
doubletemp=0;
for（i;i{
temp=i*h;
x[i]=a+temp;
}
/*cout<<"x的值为：
"<for（i=0;i{
cout<<"x["<}
cout<}
voidget_f（）
{
chartemp[ARRAY_LEN];
cout<<"请输入f（x）的表达式;";
cin>>temp;
strcpy（f,temp）;
cout<}
voidcal_y（）
{
inti=0;
doubletemp=0;
for（i=0;i{
temp=13*（x[i]-x[i]*x[i]）*exp（-1.5*x[i]）;
y[i]=temp;
}
/*cout<<"y的值为：
"<for（i=0;i{
cout<<"y["<}
cout<}
doubleresult（）
{
doubletemp=y[0],sum;
inti=1;
for（i=1;i{
if（i==n）
temp=temp+y[i];/////当有判断条件时，要先进行判断，不满足时才进行原始计算
elsetemp=temp+2*y[i];
}
sum=h*temp/2;
returnsum;
}
voiddisplay（）
{doublem;
m=result（）;
cout<<"积分区间为["<cout<<"用"<"<}
};
intmain（）
{
Comtcom;
com.getab（）;
com.get_f（）;
com.cal_nh（）;
com.cal_x（）;
com.cal_y（）;
com.display（）;
return0;
}
//复化抛物线公式程序
#include
#include
#include
#include
usingnamespacestd;
constintARRAY_LEN（10000）;
classComsim
{protected:
doublea,b,h,n,x[ARRAY_LEN],y[ARRAY_LEN];
charf[ARRAY_LEN];
public:
voidgetab（）
{
cout<<"请输入该积分的上下限（即区间）:
";
cin>>a>>b;
cout<}
voidcal_nh（）
{
intc;
cout<<"几点的复化抛物线公式？
";
cin>>c;
n=2*c-1;
cout<<"n的值为：
"<h=（b-a）/（2*n）;
cout<<"h的值为：
"<}
voidcal_x（）
{
inti=0;
doubletemp=0;
for（i;i<2*n+1;i++）
{
temp=i*h;
x[i]=a+temp;
}
/*cout<<"x的值为：
"<for（i=0;i<2*n+1;i++）
{
cout<<"x["<if（i%5==0）
cout<}
cout<}
voidget_f（）
{
chartemp[ARRAY_LEN];
cout<<"请输入f（x）的表达式;";
cin>>temp;
strcpy（f,temp）;
}
voidcal_y（）
{
inti=0;
doubletemp=0;
for（i=0;i<2*n+1;i++）
{
temp=13*（x[i]-x[i]*x[i]）*exp（-1.5*x[i]）;
y[i]=temp;
}
/*cout<<"y的值为：
"<for（i=0;i<2*n+1;i++）
{
cout<<"y["<if（i%5==0）
cout<}
cout<}
doubleresult（）
{
doubletemp1=y[0],temp2=0,temp3=0,sum1=0,sum=0;
inti=1;
for（i=1;i<2*n+1;i++）
{
if（i==2*n）
{
temp1=temp1+y[i];/////当有判断条件时，要先进行判断，不满足时才进行原始计算
}
elseif（i%2==0）
{
temp2=temp2+2*y[i];
}
else
temp3=temp3+4*y[i];
}
sum1=temp1+temp2+temp3;
sum=h*sum1/3;
returnsum;
}
voiddisplay（）
{doublem;
m=result（）;
cout<<"积分区间为["<cout<<"用"<"<}
};
intmain（）
{
Comsimsim;
sim.getab（）;
sim.get_f（）;
sim.cal_nh（）;
sim.cal_x（）;
sim.cal_y（）;
sim.display（）;
return0;
}
五、计算结果与分析
分别运用复化梯形公式和复化抛物线公式计算该积分的结果如下图：
图1.5.1为运用复化梯形公式的计算结果截图，
图1.5.2为运用复化抛物线公式的计算结果截图。
图1.5.1复化梯形公式计算结果图
图1.5.2复化抛物线公式计算结果图
由上面计算结果可以看出运用复化抛物线公式的计算量比较大，但是利用复化梯形公式要将积分区间等分的个数要大于复化抛物线公式。
并且由该程序对于已知结果的积分进行运算，可知复化simpson积分公式计算精度要高一些。
六、参考文献
[1]谭浩强.C语言程序设计[M].北京：
清华大学出版社，2005.
一、目的意义
（1）机械设计目的是要设计出一种能达到预定功能要求、具有性能好、成本低、价值最优，能满足市场需求的机械产品。
（2）掌握牛顿插值和分段线性插值多项式的算法实现。
（3）学会分析误差和精度，熟悉运用变成语句
（4）提高用数值算法解决实际问题的能力
二、内容要求
机械设计问题：
万能拉拨机中有一个圆柱形凸轮（见图1），其底圆半径R=30cm,凸轮的上端面不在同一平面上，要根据从动杆位移变化的需要进行设计制造。
将底圆周长36等分为xi（i=0,1,…，36），每一圆弧段长为h=52.36mm,对应于每一分点的柱高为yi（i=0,1,…，36）。
为方便，将圆柱展开为平面，柱面的的顶端成为图2所示的平面曲线，并已知该曲线上的37个点的坐标（表1）。
yC
A
B
yiyi
x0
Oxix17x36x
xi
图1凸轮模型图2展开曲线
表1：
测量数据表
xi
x0
x1
x2
x3
x4
x5
x6
x7
x8
yi
502.75
520.96
525
523.6
514.3
492
451
394.6
326.5
xi
X0
x10
x11
x12
x13
x14
x15
x16
x17-x36
yi
256.7
188.6
132.1
92.2
68.9
59.6
58.2
62.24
80.45-502.75
xi=ih,x0=0,x36=1884.96mm,h=52.36mm。
是直线段，AB是曲线段，为了数控加工，需要计算出圆周上任一点处的柱高,试构造算法、设计程序、编程计算。
三、问题解决的方法与算法
方法：
Newton差值和分段线性插值
3.1Newton插值
3.1.1差商公式
设y=f（x）在节点处的函数值已知，则
关于点
的一阶差商为：
则
关于点
的二阶差商为：
则
关于点
的k阶差商为：
一阶差商
二阶差商
三阶差商
X0
f（X0）
X1
f（X1）
f[X0,X1]
X2
f（X2）
f[X1,X2]
f[X0,X1,X2]
X3
f（X3）
f[X2,X3]
f[X1,X2,X3]
f[X0,X1,X2X3]
3.1.2插值公式
3.1.3插值余项
3.2分段线性差值
3.2.1插值函数
设y=f（x）在节点
处的函数值为
为了提高近似程度，可以考虑用分段线性插值来逼近原函数，这时的插值函数为分段函数：
在区间
上的线性函数为：
3.2.2插值公式
3.2.3插值余项
3.3差值算法
第一步：
输入要进行计算的x的值
第二步：
输入各组y[i]及x[i]的值，i=0,1,2,…,n
第三步：
判断条件：
if（x>=0&&x<=16*52.36）则计算各阶插商，用牛顿插值进行计算if（x>=52.36*17&&x<=1884.96）则用分段线性插值进行计算，如果两种情况均不是则输出输入错误
第四步：
计算出柱高S
第五步：
输出S
四、计算程序
4.1用牛顿插值和分段线性插值的程序
//nt.cpp:
定义控制台应用程序的入口点。
//
//#include"StdAfx.h"
#include
#defineARRAY_LEN18
doubledq（doublea）
{
inti=0;
intj=0;
doublemul=1;
doubleh=52.36;
doublex[ARRAY_LEN]={0};
doubley[ARRAY_LEN]={502.75,520.96,525,523.6,514.3,492,451,394.6,326.5,256.7,188.6,132.1,92.2,68.9,59.6,58.2,62.24};
doubletemp=0;
for（i=0;i{
temp=i*h;
x[i]=temp;
}
printf（"请输出所有x的值：
\n"）;
for（i=0;i{
printf（"%f",x[i]）;
if（i%3==0）
printf（"\n"）;
}
printf（"\n"）;
printf（"请输出所有y的值：
\n"）;
for（j=0;j{
printf（"%lf",y[j]）;
if（j%3==0）
printf（"\n"）;
}
printf（"\n"）;
///////////////////求差商
//printf（"各阶插商值为：
\n"）;
for（i=0;i{
for（j=ARRAY_LEN-1;j>i;j--）
{
y[j]=（y[j]-y[j-1]）/（x[j]-x[j-1-i]）;
}
}
///////////////////////求和
doubles=y[0];
for（j=1;j{
mul=1;
for（i=0;i<=j-1;i++）
{
mul=mul*（a-x[i]）;
}
s=s+mul*y[j];
}
returns;
}
doubles（doublen）
{
inti=0;
intj=0;
doubles=0;
doublex[2]={890.12,1884.96};
doubley[2]={80.45,502.75};
printf（"请输出所有x的值：
\n"）;
for（i=0;i<2;i++）
{
printf（"%f",x[i]）;
if（i%2==0）
printf（"\n"）;
}
printf（"\n"）;
printf（"请输出所有y的值：
\n"）;
for（j=0;j<2;j++）
{
printf（"%f",y[i]）;
if（i%2==0）
printf（"\n"）;
}
/////////////线形函数
for（i=1;i<2;i++）
{
s=y[i-1]*（（n-x[i]）/（x[i-1]-x[i]））+y[i]*（（n-x[i-1]）/（x[i]-x[i-1]））;
}
returns;
}
intmain（）
{
doubleM=0;
doubleN=0;
doubleS=0;
doublex=0;
printf（"请先判断x的个数，如果超出宏定义ARRAY_LEN的值，请先修改宏定义的值！
\n"）;
printf（"请输入所要求的x的值：
\n"）;
scanf（"%lf",&N）;
if（N>=0&&N<=16*52.36）
{
M=dq（N）;
printf（"\n"）;
printf（"x=%f时的柱高的近似值为：
%lf\n",N,M）;
}
elseif（N>=52.36*17&&N<=1884.96）
{
S=s（N）;
printf（"\n"）;
printf（"x=%f时的柱高的近似值为：
%lf\n",N,S）;
}
else
printf（"输入的值不能进行计算！
\n"）;
return0;
}
4.2用多项式拟合的matlab程序
4.2.1前曲线段拟合五次多项式的值
clearall;
clc;
x=0:
52.36:
16*52.36;
y=[502.75520.96525523.6514.3492451394.6326.5256.7188.6132.192.268.959.658.262.24];
b=polyfit（x,y,5）;%拟合出的五次函数的系数
fprintf（'运行结果为：
\n\n'）%%输出语句
disp（'拟合出的五次多项式P4（x）为：
'）%%输出语句
f=poly2str（b,'x'）%%将拟合后的多项式系数（双精度数组）转换为字符形式的函数
poly2sym（b）;%%将该向量转换为多项式
xx=linspace（min（x）,max（x））;%绘图用到的点的横坐标
yy=polyval（b,xx）;%拟合曲线的纵坐标
%subplot（2,2,2）;
plot（x,y,'m.',xx,yy,'-c'）;%绘图，原始数据+拟合曲线
xlabel（'圆周上任一点x'）;
ylabel（'对应的柱高'）;
legend（'原始数据','拟合曲线'）;%图示
title（'五次多项式拟合曲线'）;
holdon;
fprintf（'x=53时柱高的近似值为：
\n'）
e=[53104208409]
m=polyval（b,e）%%用于对已经拟合后的多项式系数，
%%当给出某个点时求其函数值;计算插值多项式在pi/8处的值
plot（e,m,'b*'）
4.2.2线性拟合后半段
x=[52.36*171884.96];
y=[80.45502.75];
b=polyfit（x,y,1）;%拟合出的1次函数的系数
f=poly2str（b,'x'）%%将拟合后的多项式系数（双精度数组）转换为字符形式的函数
poly2sym（b）;%%将该向量转换为多项式
xx=linspace（min（x）,max（x））;%绘图用到的点的横坐标
yy=polyval（b,xx）;%拟合曲线的纵坐标
plot（x,y,'m.',xx,yy,'-c'）;%绘图，原始数据+拟合曲线
xlabel（'圆周上任一点x'）;
ylabel（'对应的柱高'）;
legend（'原始数据','拟合曲线'）;%图示
title（'一次线性拟合曲线'）;
holdon;
x=[52.36*1752.36*2152.36*2952.36*35]
fprintf（'x对应柱高的近似值为：
\n'）
m=polyval（b,e）%%用于对已经拟合后的多项式系数%%当给出某个点时求其函数值
plot（e,m,'b*'）
五、计算结果与分析
5.1用牛顿插值和分段线性插值的运行结果
图5.1.1
时程序的运行结果
图5.1.2
时程序的运行结果
5.2用多项式拟合的matlab程序运行结果
5.2.1前曲线段拟合五次多项式的结果
下面为自己给定的x值经过五次多项式拟合后计算的结果，图5.2.1为拟合的曲线图。
运行结果为：
拟合出的五次多项式P5（x）为：
f=-1.678e-011x^5+3.6684e-008x^4-2.4645e-005x^3+0.0041844x^2
-0.042005x+506.7499
x=53104208409时柱高的近似值为：
m=
512.8908524.0042519.3960337.8354
图5.2.1五次多项式拟合图
5.2.2线性拟合后半段结果
下面为自己给定的x值经过线性拟合后计算的结果，图5.2.2为拟合的结果图。
f=
0.42449x-297.3974
x=
1.0e+003*
0.89011.09961.51841.8326
x对应柱高的近似值为：
m=
80.4500169.3553347.1658480.5237
图5.2.2一次线性拟合结果图
5.3结果分析
在求解该题的过程中，我分别运用了C语言和matlab两个软件进行编程实现，在C语言中，在
时，运用了牛顿插值多项式来进行逼近原曲线，为了验证程序的准确性，在这里取x=53来进行计算，计算结果为512.241797在[520.96525]之间，且接近520.96，即证明了程序的准确性。
同理，可以取任意
来验证，在这里取x=1880进行验证，运行结果见图5.1.2。
在matlab中，进行了不同区间的多项式拟合，拟合的结果见图5.2.1和图5.2.2，并且拟合的多项式表达式见上面的运行结果。
比较C语言和matlab程序可知，两者的结果较为接近，但是在matlab中可以得到拟合的多项式的表达式，并且可以可到结果图，让人更直观的看出结果。
六、
升级会员 