5《数值逼近》课程设计报告文档格式.docx
《5《数值逼近》课程设计报告文档格式.docx》由会员分享，可在线阅读，更多相关《5《数值逼近》课程设计报告文档格式.docx（31页珍藏版）》请在冰点文库上搜索。
math.h>
#include<
iomanip>
#include<
iostream>
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;
endl;
}
voidcal_nh（）
intc;
几点的复化梯形公式？
c;
n=c-1;
n的值为：
n<
h=（b-a）/n;
h的值为：
h<
voidcal_x（）
inti=0;
doubletemp=0;
for（i;
i<
n+1;
i++）
{
temp=i*h;
x[i]=a+temp;
}
/*cout<
x的值为：
for（i=0;
cout<
x["
]"
="
x[i]<
"
endl<
*/
voidget_f（）
chartemp[ARRAY_LEN];
请输入f（x）的表达式;
temp;
strcpy（f,temp）;
voidcal_y（）
temp=13*（x[i]-x[i]*x[i]）*exp（-1.5*x[i]）;
y[i]=temp;
y的值为：
y["
y[i]<
doubleresult（）
doubletemp=y[0],sum;
inti=1;
for（i=1;
if（i==n）
temp=temp+y[i];
/////当有判断条件时，要先进行判断，不满足时才进行原始计算
elsetemp=temp+2*y[i];
sum=h*temp/2;
returnsum;
voiddisplay（）
{doublem;
m=result（）;
积分区间为["
a<
b<
],被积函数为"
f<
的积分"
用"
n+1<
点复化梯形公式"
计算结果为:
setprecision（8）<
m<
};
intmain（）
Comtcom;
com.getab（）;
com.get_f（）;
com.cal_nh（）;
com.cal_x（）;
com.cal_y（）;
com.display（）;
return0;
//复化抛物线公式程序
classComsim
几点的复化抛物线公式？
n=2*c-1;
n;
h=（b-a）/（2*n）;
h;
2*n+1;
if（i%5==0）
cout<
doubletemp1=y[0],temp2=0,temp3=0,sum1=0,sum=0;
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;
点复化抛物线公式"
Comsimsim;
sim.getab（）;
sim.get_f（）;
sim.cal_nh（）;
sim.cal_x（）;
sim.cal_y（）;
sim.display（）;
}
五、计算结果与分析
分别运用复化梯形公式和复化抛物线公式计算该积分的结果如下图：
图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
X0
x10
x11
x12
x13
x14
x15
x16
x17-x36
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阶差商为：
一阶差商
二阶差商
三阶差商
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&
=1884.96）则用分段线性插值进行计算，如果两种情况均不是则输出输入错误
计算出柱高S
输出S
4.1用牛顿插值和分段线性插值的程序
//nt.cpp:
//
StdAfx.h"
#defineARRAY_LEN18
doubledq（doublea）
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};
ARRAY_LEN;
x[i]=temp;
printf（"
请输出所有x的值：
\n"
）;
printf（"
%f"
x[i]）;
if（i%3==0）
printf（"
请输出所有y的值：
for（j=0;
j<
j++）
%lf"
y[j]）;
if（j%3==0）
///////////////////求差商
//printf（"
各阶插商值为：
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;
mul=1;
for（i=0;
=j-1;
{
mul=mul*（a-x[i]）;
}
s=s+mul*y[j];
returns;
doubles（doublen）
doubles=0;
doublex[2]={890.12,1884.96};
doubley[2]={80.45,502.75};
2;
if（i%2==0）
y[i]）;
/////////////线形函数
s=y[i-1]*（（n-x[i]）/（x[i-1]-x[i]））+y[i]*（（n-x[i-1]）/（x[i]-x[i-1]））;
doubleM=0;
doubleN=0;
doubleS=0;
doublex=0;
请先判断x的个数，如果超出宏定义ARRAY_LEN的值，请先修改宏定义的值！
请输入所要求的x的值：
scanf（"
%lf"
&
N）;
if（N>
N<
=16*52.36）
M=dq（N）;
printf（"
x=%f时的柱高的近似值为：
%lf\n"
N,M）;
elseif（N>
=1884.96）
S=s（N）;
N,S）;
else
输入的值不能进行计算！
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;
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次函数的系数
一次线性拟合曲线'
x=[52.36*1752.36*2152.36*2952.36*35]
x对应柱高的近似值为：
m=polyval（b,e）%%用于对已经拟合后的多项式系数%%当给出某个点时求其函数值
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
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中可以得到拟合的多项式的表达式，并且可以可到结果图，让人更直观的看出结果。
六、
升级会员 