《数值逼近》课程设计报告

《数值逼近》课程设计报告
数值逼近课程设计报告
一、目的意义
(1)进一步熟悉掌握复化梯形和复化抛物线公式
(2)学会比较复化梯形公式和复化抛物线公式如何达到所要求的精度
(3) 提高编程能力
(4)通过数值方法求出很难求得原函数的积分和解析表达是没有明确的给出积分的近似
值
二、内容要求
积分计算问题：分别用复化梯形和复化Simpsom 求积公式计算积分
dx e x x x 5.14
02)(13-⎰-，并比较计算量（精度为10-8）。

三、问题解决的方法与算法
方法：解决该积分问题时，运用了数值积分近似解法的方法，运用复化梯形和复化
Simpsom 求积公式进行计算
3.1 复化梯形积分
3.1.1复化梯形积分公式表达式
()()()1122n n i i h T f a f b f x -=⎡⎤=++⎢⎥⎣⎦
∑ 3.1.2复化梯形积分误差表达式
[]()()()2,,12n b a h R f f a b ηη-''=-∈
3.2复化抛物线积分
3.2.1复化抛物线积分公式表达式
()()()()()11/21142n n b
k k a k k f x dx f a f x f x f b --==⎡⎤≈+++⎢⎥⎣⎦
∑∑⎰
3.2.2复化抛物线积分误差表达式
[]
()()()()()()()54444,,28802880n b a h b a R f f f a b n ηηη--=-=-∈
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<stdio.h>
//#include "stdafx.h"
#include<math.h>
#include <iomanip>
# include <iostream>
using namespace std;
const int ARRAY_LEN (10000);
class Comt
{protected:
double a,b,h,n,x[ARRAY_LEN],y[ARRAY_LEN];
char f[ARRAY_LEN];
public:
void getab()
{
cout<<"请输入该积分的上下限(即区间):";
cin>>a>>b;
cout<<endl;
}
void cal_nh()
{
int c;
cout<<"几点的复化梯形公式？";
cin>>c;
cout<<endl;
n=c-1;
cout<<"n的值为："<<n<<endl;
h=(b-a)/n;
cout<<"h的值为："<<h<<endl;
}
void cal_x()
{
int i=0;
double temp=0;
for(i;i<n+1;i++)
{
temp=i*h;
x[i]=a+temp;
}
/*cout<<"x的值为："<<endl;
for(i=0;i<n+1;i++)
{
cout<<"x["<<i<<"]"<<"="<<x[i]<<" ";
}
cout<<endl<<endl;*/
}
void get_f()
{
char temp[ARRAY_LEN];
cout<<"请输入f(x)的表达式;";
cin>>temp;
strcpy(f,temp);
cout<<endl;
}
void cal_y()
{
int i=0;
double temp=0;
for(i=0;i<n+1;i++)
{
temp=13*(x[i]-x[i]*x[i])*exp(-1.5*x[i]);
y[i]=temp;
}
/*cout<<"y的值为："<<endl;
for(i=0;i<n+1;i++)
{
cout<<"y["<<i<<"]"<<"="<<y[i]<<" ";
}
cout<<endl<<endl;*/
}
double result()
{
double temp=y[0],sum;
int i=1;
for(i=1;i<n+1;i++)
{
if(i==n)
temp=temp+y[i];/////当有判断条件时，要先进行判断，不满足时才进行原始计算
else temp=temp+2*y[i];
}
sum=h*temp/2;
return sum;
}
void display()
{ double m;
m=result();
cout<<"积分区间为["<<a<<" "<<b<<"],被积函数为"<<f<<"的积分"<<endl;
cout<<"用"<<n+1<<"点复化梯形公式"<<endl<<"计算结果为:"<<setprecision(8)<<endl<<m<<endl;
}
};
int main()
{
Comt com;
com.getab();
com.get_f();
com.cal_nh();
com.cal_x();
com.cal_y();
com.display();
return 0;
}
//复化抛物线公式程序
#include <stdio.h>
#include<iomanip>
#include<math.h>
# include <iostream>
using namespace std;
const int ARRAY_LEN (10000);
class Comsim
{protected:
double a,b,h,n,x[ARRAY_LEN],y[ARRAY_LEN]; char f[ARRAY_LEN];
public:
void getab()
{
cout<<"请输入该积分的上下限(即区间):";
cin>>a>>b;
cout<<endl;
}
void cal_nh()
{
int c;
cout<<"几点的复化抛物线公式？";
cin>>c;
n=2*c-1;
cout<<"n的值为："<<n;
h=(b-a)/(2*n);
cout<<"h的值为："<<h;
}
void cal_x()
{
int i=0;
double temp=0;
for(i;i<2*n+1;i++)
{
temp=i*h;
x[i]=a+temp;
}
/*cout<<"x的值为："<<endl;
for(i=0;i<2*n+1;i++)
{
cout<<"x["<<i<<"]"<<"="<<x[i]<<" ";
if(i%5==0)
cout<<endl;
}
cout<<endl<<endl;*/
}
void get_f()
{
char temp[ARRAY_LEN];
cout<<"请输入f(x)的表达式;";
cin>>temp;
strcpy(f,temp);
}
void cal_y()
{
int i=0;
double temp=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的值为："<<endl;
for(i=0;i<2*n+1;i++)
{
cout<<"y["<<i<<"]"<<"="<<setprecision(8)<<y[i]<<" ";
if(i%5==0)
cout<<endl;
}
cout<<endl<<endl;*/
}
double result()
{
double temp1=y[0],temp2=0,temp3=0,sum1=0,sum=0;
int i=1;
for(i=1;i<2*n+1;i++)
{
if(i==2*n)
{
temp1=temp1+y[i];/////当有判断条件时，要先进行判断，不满足时才进行原始计算
}
else if(i%2==0)
{
temp2=temp2+2*y[i];
}
else
temp3=temp3+4*y[i];
}
sum1=temp1+temp2+temp3;
sum=h*sum1/3;
return sum;
}
void display()
{ double m;
m=result();
cout<<"积分区间为["<<a<<" "<<b<<"],被积函数为"<<f<<"的积分"<<endl;
cout<<"用"<<n+1<<"点复化抛物线公式"<<endl<<"计算结果为:"<<setprecision(8)<<m<<endl;
}
};
int main()
{
Comsim sim;
sim.getab();
sim.get_f();
sim.cal_nh();
sim.cal_x();
sim.cal_y();
sim.display();
return 0;
}
五、计算结果与分析
分别运用复化梯形公式和复化抛物线公式计算该积分的结果如下图：
图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等分为x i (i=0,1, …，36)，每一圆弧段长为h=52.36mm,对应于每一分点的柱高为y i (i=0,1, …，36)。

为方便，将圆柱展开为平面，柱面的的顶端成为图2所示的平面曲线，并已知该曲线上的37个点的坐标（表1）。

B
y i
y i
x0
O x i x17x36x x i
图1 凸轮模型图2 展开曲线
x i =ih, x 0 =0, x 36=1884.96mm, h =52.36mm 。

BC 是直线段，AB 是曲线段，为了数控加工，需要计算出圆周上任一点
处的柱高,试构造算法、设计程序、编程计算。

三、问题解决的方法与算法
方法：Newton 差值和分段线性插值 3.1 Newton 插值
3.1.1差商公式
设y=f(x)在节点处的函数值已知，则()f x 关于点j i x
x ,的一阶差商为：
[]
()()j
i j i j i x x x f x f x x f --=
,
则()x f 关于点k j i x
x x ,,的二阶差商为：
[][][]
k
i k
j j i k
j
i
x x x x f x x f x x x f --=
,,,,
则()x f 关于点
k
j i x x x ,,, 的k 阶差商为
：
[]
[][]
k
k k k j i x x x x x f x x x f x x x f --=
-021110,,,,,,,,,
3.1.2 插值公式
[][][]
00010101201101
()()(),()(),,()
()
(),,n n n N x f x x x f x x x x x x f x x x x x x x x x f x x x -=+-+--++---
3.1.3 插值余项 ()()()()[]
n n n x x x x f x x x x x x x R ,,,,1010 ---=
3.2分段线性差值
3.2.1 插值函数
设y=f(x)在节点
b x x x a n =〈〈〈=...10处的函数值为()n i x f y i i ,...,2,1,0,==为了提
高近似程度，可以考虑用分段线性插值来逼近原函数，这时的插值函数为分段函数：
()()()()⎪⎪⎩
⎪⎪
⎨
⎧∈∈∈=-]12,121,01,[,.....
..........][,][,n n n x x x x s x x x x s x x x x s x s
在区间],[1i i x x -上的线性函数为：
()n
i x x x x y x x x x y x s i i i i i i i i i ,...,2,11
1
11
=--+--=----
3.2.2 插值公式
()n
i x x x x y x x x x y x s i i i i i i i i i ,...,2,11
1
11
=--+--=----
3.2.3 插值余项
()()()()i i i i i i i x x x x x x f x R 〈〈--=
--ξξ11,!2''
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<stdio.h>
#define ARRAY_LEN 18
double dq(double a)
{
int i=0;
int j=0;
double mul=1;
double h=52.36;
double x[ARRAY_LEN]={0};
double
y[ARRAY_LEN]={502.75,520.96,525,523.6,514.3,492,451,394.6,326.5,256.7,188.6,132.1,9 2.2,68.9,59.6,58.2,62.24};
double temp=0;
for(i=0;i<ARRAY_LEN;i++)
{
temp=i*h;
x[i]=temp;
}
printf("请输出所有x的值：\n");
for(i=0;i<ARRAY_LEN;i++)
{
printf("%f ",x[i]);
if(i%3==0)
printf("\n");
}
printf("\n");
printf("请输出所有y的值：\n");
for(j=0;j<ARRAY_LEN;j++)
{
printf("%lf ",y[j]);
if(j%3==0)
printf("\n");
}
printf("\n");
///////////////////求差商
//printf("各阶插商值为：\n");
for(i=0;i<ARRAY_LEN;i++)
{
for(j=ARRAY_LEN-1;j>i;j--)
{
y[j]=(y[j]-y[j-1])/(x[j]-x[j-1-i]);
}
}
///////////////////////求和
double s=y[0];
for(j=1;j<ARRAY_LEN;j++)
{
mul=1;
for(i=0;i<=j-1;i++)
{
mul=mul*(a-x[i]);
}
s=s+mul*y[j];
}
return s;
}
double s(double n)
{
int i=0;
int j=0;
double s=0;
double x[2]={890.12,1884.96};
double y[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]));
}
return s;
}
int main()
{
double M=0;
double N=0;
double S=0;
double x=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);
}
else if (N>=52.36*17&&N<=1884.96)
{
S=s(N);
printf("\n");
printf("x=%f时的柱高的近似值为：%lf\n",N,S);
}
else
printf("输入的值不能进行计算！\n");
return 0;
}
4.2用多项式拟合的matlab程序
4.2.1前曲线段拟合五次多项式的值
clear all;
clc;
x=0:52.36:16*52.36;
y=[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];
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('五次多项式拟合曲线');
hold on;
fprintf('x=53时柱高的近似值为：\n')
e=[53 104 208 409]
m=polyval(b,e)%%用于对已经拟合后的多项式系数，
%%当给出某个点时求其函数值;计算插值多项式在pi/8处的值plot(e,m,'b*')
4.2.2线性拟合后半段
x=[52.36*17 1884.96];
y=[80.45 502.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('一次线性拟合曲线');
hold on;
x=[52.36*17 52.36*21 52.36*29 52.36*35]
fprintf('x对应柱高的近似值为：\n')
m=polyval(b,e)%%用于对已经拟合后的多项式系数%%当给出某个点时求其函数值plot(e,m,'b*')
五、计算结果与分析
5.1用牛顿插值和分段线性插值的运行结果
图5.1.1
[]
017
,
x x x
∈
时程序的运行结果
图5.1.2
[]
1736
,
x x x
∈
时程序的运行结果
5.2用多项式拟合的matlab程序运行结果
5.2.1前曲线段拟合五次多项式的结果
下面为自己给定的x值经过五次多项式拟合后计算的结果，图5.2.1为拟合的曲线图。

运行结果为：
拟合出的五次多项式P5(x)为：
f =-1.678e-011 x^5 + 3.6684e-008 x^4 - 2.4645e-005 x^3 + 0.0041844 x^2
- 0.042005 x + 506.7499
x=53 104 208 409时柱高的近似值为：
m =
512.8908 524.0042 519.3960 337.8354
图5.2.1五次多项式拟合图
5.2.2线性拟合后半段结果
下面为自己给定的x值经过线性拟合后计算的结果，图5.2.2为拟合的结果图。

f =
0.42449 x - 297.3974
x =
1.0e+003 *
0.8901 1.0996 1.5184 1.8326
x对应柱高的近似值为：
m =
80.4500 169.3553 347.1658 480.5237
图5.2.2一次线性拟合结果图
5.3结果分析
在求解该题的过程中，我分别运用了C语言和matlab两个软件进行编程实现，在
C语言中，在
[]
017
,
x x x
∈
时，运用了牛顿插值多项式来进行逼近原曲线，为了验证程序
的准确性，在这里取x=53来进行计算，计算结果为512.241797在[520.96 525]之间，
且接近520.96，即证明了程序的准确性。

同理，可以取任意
[]
1736
,
x x x
∈
来验证，在这里
取x=1880进行验证，运行结果见图5.1.2。

在matlab中，进行了不同区间的多项式拟合，拟合的结果见图5.2.1和图5.2.2，并且拟合的多项式表达式见上面的运行结果。

比较C语言和matlab程序可知，两者的结果较为接近，但是在matlab中可以得到拟合的多项式的表达式，并且可以可到结果图，让人更直观的看出结果。

六、参考文献
[1] 谭浩强. C语言程序设计[M]. 北京：清华大学出版社，2005.
一、目的意义
(1)提高matlab编程能力
(2) matlab中插值函数以及其他画图函数的运用
(3)提高分析程序的能力
二、内容要求
丘陵高程测定问题
在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如下。

试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。

三、问题解决的方法与算法
方法：用插值法计算出各个节点的值，用matlab中的mesh函数绘出地貌图，即可以直观的观察。

算法：
1.输入原始数据
2.用插值函数计算出每个节点的值
3.画出地貌图
4.判断寻找最大值点标记位置并保存
5.将最大值点及其对应坐标用stem在原图中绘出
四、计算程序
程序1：
function [s,x0,y0]=high1(N)
x=[100 200 300 400];
y=[100 200 300 400]';
z=[636 697 624 478;698 712 630 478;680 674 598 412;662 626 552 334];
%mesh(x,y,z) %绘原始数据图
xx=linspace(100,400,N);%即为从100到400分成N等分
yy=linspace(100,400,N)'; %加密纵坐标数据到N个
zh=interp2(x,y,z,xx,yy,'cubic') %插值三次多项式插值:'cubic'
mesh(xx,yy,zh) %加密后的地貌图
s=0;
for i=1:N^2
if zh(i)>s
s=zh(i); %循环找最大值
n=mod(i,N);%i对N取余
m=(i-n)/N;
end
end
x0=100+300/N*m %最大值点对应的X值
y0=100+300/N*n
hold on;
stem3(x0,y0,s,'fill')
程序2：
clc;
clear all;
x=[100 100 100 100 200 200 200 200 300 300 300 300 400 400 400 400];
y=[100 200 300 400 100 200 300 400 100 200 300 400 100 200 300 400];
z=[636 697 624 478 698 712 630 478 680 674 598 412 662 626 552 334];
xi=100:5:400;
yi=100:5:400;
[X,Y]=meshgrid(xi,yi);
H=griddata(x,y,z,X,Y,'cubic')%二元非等距插值
surf(X,Y,H);
view(-120,23);
hold on;
maxh=vpa(max(max(H)),6)%控制变量计算结果的显示位数
[r,c]=find(H>=single(maxh))% 检索R中最大元素所在的位置（行标r和列标c）
stem3(X(r,c),Y(r,c),maxh,'fill')
五、计算结果与分析
下面为两个不同matlab程序的运行结果，程序1是用三次多项式插值编写函数实现的，程序2是用二元非等距插值编写M文件实现的。

在程序1中，由于存在等分区间的问题，即给定不同的区间个数，会产生不同的结果，但两个结果的值分别为718.1243，716.526结果较为接近，即证明两个程序的正确性。

两个程序的运行结果见图5.1和图5.2。

程序1运行结果：
>> high(100)
x0 =
166
y0 =
196
ans =
718.1243
图5.1 程序1丘陵地貌图以及最高值位置
程序2运行结果：
maxh =
716.526
r =
16
c =
20
图5.2程序2丘陵地貌图以及最高值位置图
六、参考文献
[1] 谭浩强. C语言程序设计[M]. 北京：清华大学出版社，2005.。