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

课程设计报告

课程名称数值逼近

专业信息与计算科学

班级计算091

姓名

学号

指导教师秦新强、胡钢

日期2011-07-01

理学院应用数学系

一、目的意义

(1) 进一步掌握复化梯形积分公式，提高运用能力。

(2) 进一步掌握复化抛物线积分公式，提高运用能力。 (3) 进一步掌握c 语言，提高编程能力。

二、内容要求

积分计算问题：分别用复化梯形和复化Simpson 求积公式计算积分

dx e x x x 5.14

2)(13-?

-，并比较计算量（精度为10-8）。

三、问题解决的方法与算法

3.1解决方法：

利用复化梯形和复化Simpsom 积分公式。

（1）根据两个方程式的误差估计式，求出各自满足精度的等分份数n 。

（2）带入所编写的c 语言程序，检验并计算出结果。 3.2算法设计：

（1）输入所求的等分数n 。

（2）调用复化积分公式，算出n+1点和n 点的积分估计数值。 （3）用n+1点和n 点的数值检验，比较数值是否达到精度要求。 （4）若达不到精度要求，输出提示，结束程序。

（5）若达到精度要求，输出积分近似结果，程序结束。

四、计算程序

4.1复化梯形积分程序源代码： #include #include

#define f(x)(13*((x)-(x)*(x))*exp(-1.5*(x))) #define epsilon 0.00000001

double Tixing(double aa,double bb,int n){ int i;

double fz;

double h=(bb-aa)/n; fz=(f(aa)+f(bb))/2; for(i=1;i

数值积分及其应用

报告

1

{fz+=f(aa+h*i);}

fz*=h;

return fz;

}

int main(){

double a=0.0;

double b=4.0;

int n;

printf("请输入等分值：");

scanf("%d",&n);

printf("\n");

double T1;

T1=Tixing(a,b,n);

double T2;

T2=Tixing(a,b,(n+1));

if(fabs(T2-T1)>0.00000001){

printf("达不到精度要求，请检查等分值是否求解正确！\n");

}

else printf("用复化梯形积分方法所计算的积分的数值为：%.8lf\n",T1);

return 0;

}

4.2复化Simpson积分程序源代码:

#include

#include

#define f(x)(13*((x)-(x)*(x))*exp(-1.5*(x)))

#define epsilon 0.00000001

double Simp(double aa,double bb,int n){

int i;

double S_n;

double h=(bb-aa)/(2*n);

double fz,fm,fk;

fz=f(aa)+f(bb);

fm=0;

fk=0;

for(i=1;i

{fm+=f(aa+(2*i-1)*h);}

for(i=1;i

{fk+=f(aa+(2*i)*h);}

S_n=h*(fz+4*fm+2*fk)/3;

return S_n;

}

void main()

{ int n;

double a=0.0;

double b=4.0;

printf("请输入等分值：");

scanf("%d",&n);

double T1;

T1=Simp(a,b,n);

double T2;

T2=Simp(a,b,(n+1));

if(fabs(T2-T1)>0.00000001){

printf("达不到精度要求，请检查等分值是否求解正确！\n");

}

else printf("用复化Simpson积分方法所计算的积分的数值为：%.8lf\n",T1);

}

五、计算结果与分析

5.1

复化梯形积分计算结果与分析：

利用复化梯形积分的误差估计式

()()

()

3

2

22

2

|[]|||,max||

12a x b

b a

Rn f M M f x

n≤≤

-

≤=

，求出n

的值，带入公式计算，是满足精度要求的。

5.2复化Simpson积分计算结果与分析：

这里同样利用复化Simpson的误差估计式，求出n的值，带入公式计算，求出数值。

比较两种方法，可见满足同样精度的积分算式，复化Simpson比复化梯形要求的等分数小得多。

六、参考文献

[1] 谭浩强. C语言程序设计[M]. 北京:清华大学出版社,2007.

[2] 秦新强. 数值逼近[M]. 西安:西安理工大学印刷厂,2011.

一、目的意义

(1)利用有限个离散点构造n 次Lagrange 插值计算式来解决实际问题。 (2)掌握利用分段线性插值计算式来解决实际问题。

(3)进一步提高c 语言编程能力，提高实际问题分析能力。

二、内容要求

机械设计问题：万能拉拨机中有一个圆柱形凸轮(见图1)，其底圆半径R=30cm,凸轮的上端面不在同一平面上，要根据从动杆位移变化的需要进行设计制造。

将底圆周长36等分为x i (i =0,1, …， 36)，每一圆弧段长为h=52.36mm,对应于每 一分点的柱高为y i (i =0,1, …， 36)。为方便，将圆柱展开为平面，柱面的的顶端成为图2所示的平面曲线，并已知该曲线上的37个点的坐标（表1）。

y C

A

B y i y i

x 0

O x i x 17 x 36 x

x i

图1 凸轮模型 图2 展开曲线

表1：测量数据表

x i x 0 x 1 x 2 x 3 x 4 x 5 x 6

x 7 x 8 y i 502.75 520.96 525 523.6 514.3 492 451 394.6 326.5 x i X 0 x 10 x 11 x 12 x 13 x 14 x 15 x 16 x 17- x 36 y i

256.7 188.6 132.1 92.2 68.9 59.6

58.2

62.24 80.45-502.75

x i =jh, x 0 =0, x 36=1884.96mm, h =52.36mm 。

BC 是直线段，AB 是曲线段，为了数控加工，需要计算出圆周上任一点处的柱

高,试构造算法、设计程序、编程计算。

三、问题解决的方法与算法

机械设计问题

报告

2

3.1解决方法：

（1）利用Excle 算出 0x

到 36x

各点的值；

x 0 x 1 x 2 x 3 x 4 x 5 x 6 x 7 0 52.36 104.72 157.08 209.44 261.8 314.16 366.52 x 8 x 9 x 10 x 11 x 12 x 13 x 14 x 15 418.88 471.24 523.6 575.96 628.32 680.68 733.04 785.4 x 16 x 17 x 18 x 19 x 20 x 21 x 22 x 23

837.76 890.12 942.48 994.84 1047.2 1099.56 1151.92 1204.28 x 24 x 25 x 26 x 27 x 28 x 29 x 30 x 31 1256.64 1309 1361.36 1413.72 1466.08 1518.44 1570.8 1623.16 x 32 x 33 x 34 x 35 x 36

1675.52 1727.88 1780.24 1832.6 1884.96

（2）对于前17个点和后20个点分别构造Lagrange 插值函数和分段线性插值函数。 3.2算法设计：

（1）输入要计算的坐标点值。

（2）用if 语句进行控制，对坐标值进行判断裁定，选择相应的插值计算式进行。 （3）输出计算的圆柱体高度，结束程序。 四、计算程序

#include #include void main(){ double a; printf("********请输入要求的点的坐标值（范围为0到1884.96mm ）："); scanf("%lf",&a);

if(a>=0&&a<=837.76){ int i,j,n;

double tmp=1,lagrange=0,e; double fm=1,fz=1; double l[20]; n=16; double x[17]={0,52.36,104.72,157.08,209.44,261.8,314.16,366.52,418.88,575.96, 628.32,680.68,733.04,785.4,837.76}, y[17]={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}; e=a; for(j=0;j<=n;j++){ fm=1;fz=1;

for(i=0;i<=n;i++){

表2-3-1 各坐标点对应的数值

if(i!=j)

{fm*=e-x[i];fz*=x[j]-x[i];}

l[j]=fm/fz;

}

}

for(j=0;j<=n;j++)

lagrange=lagrange+l[j]*y[j];

printf("\n得到圆柱体的高度为：%lf（mm）\n",lagrange);

}

if(a>837.76){

double h=a;

double q[2]={890.12,1884.96},w[2]={80.45,502.75};

double fq=(h-q[1])/(q[0]-q[1]);

double fw=(h-q[0])/(q[1]-q[0]);

double s=w[0]*fq+w[1]*fw;

printf("\n得到圆柱体的高度为：%lf（mm）\n",s);

}

}

五、计算结果与分析

在程序中任意输入范围内的一点，得到（截图）：

对于计算的结果，前17个点用了Lagrange插值法，计算出的高度肯定存在一定的误差，但是可以满足数控加工的精度要求。后面的点运用点斜式编程求出计算式，在认为数据没有误差的前提下，求出的圆柱体的高度精确成立。

c语言程序运用if结构，来选择输入的坐标值点满足哪个计算式，从而选择正确的计算式，计算出相应的近似结果。

六、参考文献

[1] 谭浩强. C语言程序设计[M]. 北京:清华大学出版社,2007.

[2] 秦新强. 数值逼近[M].西安:西安理工大学印刷厂,2011.

一、目的意义

（1）掌握三次样条的推导过程，并熟练运用。

（2）熟练运用三弯矩算法来构造三次样条插值函数。

（3）了解运用Matlab作三次周期样条函数轮廓图。

（4）熟练运用c++编程来解决问题，学会运用追赶法来解方程式，并会正确设计程序。

二、内容要求

三、问题解决的方法和算法

3.1解决方法：

（1）根据三弯矩算法，利用周期函数的条件，求解出弯矩向量

,1,2,...,18

i

M i

。

（2）带入三次样条函数计算式，求解出每个小区间上的插值关系式。

（3）利用Matlab，画出三次样条周期函数的近似图。

3.2算法设计：

（1）运用追赶法，求出弯矩向量。

（2）求出每个小区间上的三次样条函数插值，并输出。

（3）检验，输入一个坐标值，输出两个函数值，程序结束。

四、计算程序

#include

#include

#include

using namespace std;

报告3综合应用问题

const int MAX = 50;

const double x[19]={100,134,164,180,198,195,186,160,136,100,66,35,15,0,5,17,32,63,100}; const double y[19]={503,525,514.3,451.0,326.5,188.6,92.2,59.6,62.2,102.7,147.1,191.6, 236.0,280.5,324.9,369.4,413.8,458.3,503};

float h[MAX];

float c[MAX], a[MAX], g[MAX];

float f(int x1, int x2, int x3){

float a=(y[x3]-y[x2])/(x[x3]-x[x2]);

float b=(y[x2]-y[x1])/(x[x2]-x[x1]);

return (a-b)/(x[x3]-x[x1]);

}

void M(int n){ //用追赶法求解出弯矩向量M

float B[18];

B[0]=c[0]/2;

for(int i=1;i

B[i]=c[i]/(2-a[i]*B[i-1]);

for(i=1;i<=n;i++)

g[i]=(g[i]-a[i]*g[i-1])/(2-a[i]*B[i-1]);

for(i=n-1;i>=0;i--)

g[i]=g[i]-B[i]*g[i+1];

}

void printout(int n);

void main()

{

int n=18, i;

for(i=0;i

h[i]= x[i+1]-x[i];

for(i=1;i

g[i]=6 * f(i-1,i,i+1);

g[0]=(h[18]/6*g[17]+(503-458.3)/(100-63)+h[1]/6*g[1]-11/17)/(-h[1]/3-h[18]/6);

g[n]=g[0];

for(i=1;i

a[i]=h[i]/(h[i]+h[i+1]);

c[i]=1-a[i];

}

a[n]=h[n]/(h[1]+h[n]);

c[n]=h[1]/(h[n]+h[1]);

M(n);

cout<<"\n输出三次样条插值函数：\n";

printout( n);

double u;

cout<<"请输入横坐标x的值：";

cin>>u;

double y1;

double y2;

int j=0;

for(i=0;i

if(u>=x[i]&&u<=x[i+1])

y1=g[i]/(6*h[i])*(x[i+1]-u)*(x[i+1]-u)*(x[i+1]-u)+g[i+1]/(6*h[i])*(u-x[i])*(u-x[i])*(u-x[i ] )+(y[i]-g[i]*h[i]*h[i]/6)/h[i]*(x[i+1]-u)+(y[i+1]-g[i+1]*h[i]*h[i]/6)/h[i]*(u-x[i]);

}

for(i=0;i

if(u<=x[i]&&u>=x[i+1])

y2=fabs(g[i]/(6*h[i])*(u-x[i+1])*(u-x[i+1])*(u-x[i+1])+g[i+1]/(6*h[i])*(x[i]-u)*(x[i]-u)* (x[i]-u)+(y[i] - g[i]*h[i]*h[i]/6)/h[i]*(u- x[i+1])+(y[i+1]-g[i+1]*h[i]*h[i]/6)/h[i]*(x[i]-u));

}

cout<<"x对应的函数值为:"<

cout<<"y1="<

}

void printout(int n){

cout<

for(int i=0;i

cout<<"["<

float t=g[i]/(6*h[i]);

if(t>0)

cout<

else

cout<<-t<<"*(x-"<

t=g[i+1]/(6*h[i]);

if(t>0)

cout<<"+"<

else

cout<<"-"<<-t<<"*(x-"<

t=(y[i]-g[i]*h[i]*h[i]/6)/h[i];

if(t>0)

cout<<"+"<

else

cout<<"-"<<-t<<"*("<

t=(y[i+1]-g[i+1]*h[i]*h[i]/6)/h[i];

if(t>0)

cout<<"+"<

else

cout<<"-"<<-t<<"*(x-"<

cout<

} }

五、计算结果与分析

利用Matlab 作出周期样条的轮廓图：

20

40

60

80

100

120

140

160

180

200

50100150200250300350400450500

550

求出的周期三次样条插值函数为（截图）：

图3-5-1 周期样条函数轮廓图

从程序运算结果截图中可以看出，作者对三次样条函数进行了检验，随便输入一个范围内的坐标值，得到两个函数值，与图3-5-1的函数值大体一致。可见求出的三次样条插值函数计算式近似正确。

六、参考文献

[1] 姚全珠. C++面向对象程序设计[M]. 北京:电子工业出版社,2010.

[2] 秦新强. 数值逼近[M]. 西安:西安理工大学印刷厂,2011.