西电计算方法上机实习题

计算方法实验报告

班级： 学号： 姓名： 成绩：

1 舍入误差及稳定性

一、实验目的

（1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性

二、实验内容

1、用两种不同的顺序计算10000

2

1n n -=∑，分析其误差的变化

2、已知连分数()

1

01223//(.../)n n a f b b a b a a b =+

+++，利用下面的算法计算f ：

1

1

,i n n i i i a d b d b d ++==+

(1,2,...,0i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分

1

041

n

n x y dx x =+? (0,1,...,1

n = 4、设2

2

11

N

N j S j ==

-∑

，已知其精确值为1311221N N ??

-- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序

（3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数

三、实验步骤、程序设计、实验结果及分析

1、用两种不同的顺序计算10000

2

1n n -=∑，分析其误差的变化

（1）实验步骤：

分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算

#include

#include

void main()

{

double sum=0;

int n=1;

while(1)

{

sum=sum+(1/pow(n,2));

if(n%1000==0)printf("sun[%d]=%-30f",n,sum);

if(n>=10000)break;

n++;

}

printf("sum[%d]=%f\n",n,sum);

}

b.逆序计算

#include

#include

void main()

{

double sum=0;

int n=10000;

while(1)

{

sum=sum+(1/pow(n,2));

if(n%1000==0)

printf("sum[%d]=%-30f",n,sum);

if(n<=1)break;

n--;

}

printf("sum[%d]=%f\n",n,sum);

}

（3）实验结果及分析：

程序运行结果：

a.顺序计算

b.逆序计算

结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，而逆序计算误差最开始十分大，后来结果正确。 2、已知连分数()

1

01223//(.../)n n a f b b a b a a b =++++，计算f ：

（1）实验步骤：

利用 1

1

,i n n i i i a d b d b d ++==+

(1,2,...,0i n n =--，0f d =，计算f （2）程序设计 #include #include void main() { int i=0,n; float a[1024],b[1024],d[1024]; printf("please input n,n="); scanf("%d",&n); printf("\nplease input a[1] to a[n]:\n"); for(i=1;i<=n;i++) {

printf("a[%d]=",i); scanf("%f",&a[i]); } printf("\nplease input b[0] to b[n]:\n"); for(i=0;i<=n;i++) { printf("b[%d]=",i); scanf("%f",&b[i]); }

d[n]=b[n]; for(i=n-1;i>=0;i--) d[i]=b[i]+a[i+1]/d[i+1]; printf("\nf=%f\n",d[0]); }

（3）实验结果 程序运行结果：

3、给出一个有效的算法和一个无效的算法计算积分

1

041

n

n x y dx x =+? (0,1,...,1

n = （1）实验步骤

利用C 语言编写程序，分别使用数值稳定的和数值不稳定的计算公式所建立的递推公式进行计算。 （2）程序设计 #include #include main() { double y_0=(1/4.0)*log(5),y_1; double y_2=(1.0/55.0+1.0/11.0)/2,y_3; int n=1,m=10; printf("有效算法输出结果：\n"); printf("y[0]=%-20f",y_0); while(1) { y_1=1.0/(4*n)+y_0/(-4.0); printf("y[%d]=%-20f",n,y_1); if(n>=10) break; y_0=y_1; n++; if(n%3==0) printf("\n"); } printf("\n 无效算法的输出结果：\n"); printf("y[10]=%-20f",y_2);

while(1) { y_3=1.0/n-4.0*y_2; printf("y[%d]=%-20f",m-1,y_3); if(m<=1) break; y_2=y_3; m--; if(m%2==0) printf("\n"); } }

（3）实验结果及分析 程序运行结果：

结果分析：无效算法数值不稳定，误差造成的影响特别大 4、设22

11

N

N j S j ==

-∑

，已知其精确值为1311221N N ??

-- ?+??

（1）实验步骤

先编程按从大到小的顺序计算N S 的程序，再编程按从小到大的顺序计算N S 的程序，然后按两种顺序分别计算10001000030000,,S S S 。

（2）程序设计 #include main() { int N; double SN[30000]; SN[30000]=(3.0/2.0-1.0/30000.0-1/30001.0)/2.0; for(N=30000;N>=2;N--) SN[N-1]=SN[N]-1.0/(N*N-1); printf("从大到小顺序计算：

\nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); SN[2]=(3.0/2-1.0/2.0-1/3.0)/2.0; for(N=3;N<=30000;N++) SN[N]=SN[N-1]+1.0/(N*N-1); printf("从小到大顺序计算：

\nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); }

（3）实验结果及分析 程序运行结果：

结果分析：不同顺序计算所得结果是一样的。

四、总结

通过这次上机，学习了解了舍入误差在不同算法时对结果的影响不同，稳定的算法才能获得正确的结果。

2 方程求根

一、实验目的

（1）通过对二分法与牛顿迭代法做编程练习和上机运算，进一步体会二分法和牛顿法的不同。

（2）编写割线迭代法的程序，求非线性方程的解，并与牛顿迭代法作比较。

二、实验内容

1、用牛顿法求下列方程的根 （1）2

0x

x e -= （2）10x xe -= （3）lg 20x x +-=

2、编写割线法程序求解第一问的方程

三、实验步骤、程序设计、实验结果及分析

1、牛顿法

（1）实验步骤

通过定义牛顿法求方程的子函数，用main 函数调用子函数求根 （2）程序设计 #include #include

typedef float (*p)(float );

float ff1(float x) {

return x*x-exp(x);

}

float ff2(float x)

{

return x*exp(x)-1;

}

float ff3(float x)

{

return log(x)+x-2;

}

float answer(float(*p)(float))

{

int k=2;

float m=1,n=-1,x2,a,b,c;

if (p==ff3)n=2;

printf("x[0] = %.4f, x[1] = %.4f, ",m,n);

while (1)

{

if (fabs(m-n)<1e-4) break;

a=p(n)*(n-m);

b=p(n)-p(m);

c=a/b;

x2=n-c;

m = n;

n = x2;

printf("x[%d] = %.4f, ",k,x2);

k++;

if (k%3==0) printf("\n");

}

if (k%3!=0) printf("\n");

printf("iteration times: %d, roots: %.4f\n ",k-2,n);

return 0;

}

main()

{

printf("x*x-exp(x),\n");

answer(ff1);

printf("x*exp(x)-1,\n");

answer(ff2);

printf("lg(x)+x-2,\n");

answer(ff3);

return 0;

}

（3）实验结果及分析

2、割线法

（1）程序设计

#include

#include

float gexian(float,float);

float f(float);

main()

{

int i,j;

float x1=2.2;

float x2=2,x3;

scanf("%d",&i);

if(i==1) printf("%f",x1); else if(i==2) printf("%f",x2); else

{

for(j=3;j<=i;j++)

{

x3=gexian(x1,x2);

x1=x2;

x2=x3;

}

printf("%f",gexian(x1,x2)) ; }

}

float f(float x)

{

return (x*x-exp(x));

}

float gexian(float x1,float x2)

{

return (x2-(f(x2)/(f(x2)-f(x1)))*(x2-x1));

}

（3）实验结果及分析

四、总结

了解和学习了二分法和牛顿迭代法的思想以及程序设计的方法，比较了迭代法和牛顿法的特点：牛顿法收敛速度较快，但对初值选取要求较高；割线法计算量少。

3 线性方程组数值解法

一、实验目的

（1）熟悉求解线性方程组的有关理论和方法；

（2）会编制列主元消去法，LU分解法，雅可比及高斯-赛德尔迭代法的程序；

（3）通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。

二、实验内容

1、用列主元消去法解方程组

2、用LU分解法解方程组

三、实验步骤、程序设计、实验结果及分析

1、用列主元消去法解方程组

（1）程序设计

#include

#include

void ColPivot(float*,int,float[]);

void ColPivot(float*c,int n,float x[])

{

int i,j,t,k;

float p;

for(i=0;i<=n-2;i++)

{

k=i;

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

if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))

k=j;

if(k!=i)

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

{

p=*(c+i*(n+1)+j);

*(c+i*(n+1)+j)=*(c+k*(n+1)+j);

*(c+k*(n+1)+j)=p;

}

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

{

p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));

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

*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));

}

}

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

{

for(j=n-1;j>=i+1;j--)

(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));

x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));

}

}

void main()

{

int i;

float x[4];

float c[4][5]={1,1,0,3,4,2,1,-1,1,1,3,-1,-1,3,-3,-1,2,3,-1,4};

ColPivot(c[0],4,x);

for(i=0;i<=3;i++)

printf("x[%d]=%f\n",i,x[i]);

}

（2）实验结果及分析

(1)题

(2)题

2、用LU分解法解方程组

（1）程序设计

#include

void main()

{

float x[4];

int i;

float a[4][5]={48,-24,0,-12,4,

-24,24,12,12,4,

0,6,20,2,-2,

-6,6,2,16,-2

};

void DirectLU(float*,int,float[]);

DirectLU(a[0],4,x);

for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]);

}

void DirectLU(float*u,int n,float x[])

{

int i,r,k;

for(r=0;r<=n-1;r++)

{

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

for(k=0;k<=r-1;k++)

*(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));

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

{

for(k=0;k<=r-1;k++)

*(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r));

*(u+i*(n+1)+r)/=*(u+r*(n+1)+r);

}

}

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

{

for(r=n-1;r>=i+1;r--)

*(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r];

x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i));

}

}

四、总结

掌握了用列主元消去法和LU分解法求解方程组程序编写的技巧。

4 插值法

一、实验目的

（1）熟悉拉格朗日插值法多项式和牛顿插值多项式，注意其不同点；（2）掌握三次样条插值解决一些实际问题。

二、实验内容

1、按所给数据做二次插值，并求给定点的函数值

2、按所给数据做五次插值，并求给定点的函数值

3、牛顿前插公式计算函数值

三、实验步骤、程序设计、实验结果及分析

1、二次插值

（1）程序设计

#include

float Lagrange(float x[],float y[],float xx,int n) //n为（n+1）次插值；{

int i,j;

float *a,yy=0;

a=new float[n];

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

{

a[i]=y[i];

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

if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);

yy+=a[i];

}

delete a;

return yy;

}

void main()

{

float x[5]={-3.0,-1.0,1.0,2.0,3.0};

float y[5]={1.0,1.5,2.0,2.0,1.0};

float xx1=-2,xx2=0,xx3=2.75,yy1,yy2,yy3;

yy1=Lagrange(x,y,xx1,3);

yy2=Lagrange(x,y,xx2,3);

yy3=Lagrange(x,y,xx3,3);

printf("x1=%-20f,y1=%f\n",xx1,yy1);

printf("x2=%-20f,y2=%f\n",xx2,yy2);

printf("x3=%-20f,y3=%f\n",xx3,yy3);

}

（2）实验结果

2、五次插值

（1）程序设计

#include

float Lagrange(float x[],float y[],float xx,int n) //n为（n+1）次插值；{

int i,j;

float *a,yy=0;

a=new float[n];

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

{

a[i]=y[i];

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

if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);

yy+=a[i];

}

delete a;

return yy;

}

int main()

{

float x[6]={0.30,0.42,0.50,0.58,0.66,0.72};

float y[6]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};

float xx1=0.46,xx2=0.55,xx3=0.60,yy1,yy2,yy3;

yy1=Lagrange(x,y,xx1,6);

yy2=Lagrange(x,y,xx2,6);

yy3=Lagrange(x,y,xx3,6);

printf("x1=%-20f,y1=%f\n",xx1,yy1);

printf("x2=%-20f,y2=%f\n",xx2,yy2);

printf("x3=%-20f,y3=%f\n",xx3,yy3);

}

（2）实验结果

3、牛顿前插公式计算函数值

（1）程序设计

#include

#define N 3

void Difference(float y[],float f[4][4],int n)

{

int k,i;

f[0][0]=y[0];f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];

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

for(i=0;i<=(N-k);i++)

f[i][k]=f[i+1][k-1]-f[i][k-1];

return;

}

int main()

{

int i,k=1;

float a,b=1,m=21.4,t=1.4,f[4][4]={0};

float x[5]={20,21,22,23,24};

float y[5]={1.30103,1.32222,1.34242,1.36173,1.38021};

Difference(y,f,N);

a=f[0][0];

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

{

k=k*i;

b=b*(t-i+1);

a=a+b*f[0][i]/k;

}

printf("x(k)\n");

for (i=0;i<=4;i++)

printf( "%-20f",x[i]);

printf("\ny(k)\n");

for (i=0;i<=4;i++)

printf("%-20f",y[i]);

for(k=1;k<=3;k++)

{

printf("\nF(%d)\n ",k);

for(i=0;i<=(3-k);i++)

{

printf("%-20f",f[i][k]);

}

}

printf ("\n");

printf("f(%f)=%-20f",m,a);

printf ("\n");

}

（2）实验结果

四、总结

学习了插值法，学会了利用插值法编程求多项式的解，可以求解很多问题，让求解多项式解变得非常简单。

5 曲线拟合

一、实验目的

（1）了解最小二乘法的基本原理，通过计算机解决实际问题； （2）了解超定方程组的最小二乘解法。

二、实验内容

1、分别用抛物线2y a bx cx =++和指数曲线bx y ae =拟合所给数据，并比较这两个拟合函数的优劣。

2、按所给实验数据，用形如2y a bx =+的抛物线进行最小二乘拟合。 三、程序设计、结果分析

1、分别用抛物线2y a bx cx =++和指数曲线bx y ae =拟合所给数据

a.抛物线

（1）程序设计： #include #include int main() { int i; float a[3]; float x[15]={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8}; float

y[15]={33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,301.40,310.40,318.15,325.15}; void Approx(float[],float[],int,int,float[]); Approx(x,y,15,2,a); for(i=0;i<=2;i++) printf("a[%d]=%f\n",i,a[i]); }

void Approx(float x[],float y[],int m,int n,float a[]) { int i,j,t; float *c=new float[(n+1)*(n+2)]; float power(int,float); void ColPivot(float *,int,float[]); for(i=0;i<=n;i++) { for(j=0;j<=n;j++) {

*(c+i*(n+2)+j)=0;

for(t=0;t<=m-1;t++)

*(c+i*(n+2)+j)+=power(i+j,x[t]);

}

*(c+i*(n+2)+n+1)=0;

for(j=0;j<=m-1;j++)

*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);

}

ColPivot(c,n+1,a);

delete c;

}

void ColPivot(float *c,int n,float x[])

{

int i,j,t,k;

float p;

for(i=0;i<=n-2;i++)

{

k=i;

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

if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))

k=j;

if(k!=i)

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

{

p=*(c+i*(n+1)+j);

*(c+i*(n+1)+j)=*(c+k*(n+1)+j);

*(c+k*(n+1)+j)=p;

}

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

{

p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));

for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));

}

}

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

{

for(j=n-1;j>=i+1;j--)

(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));

x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));

}

}

float power(int i,float v)

{

float a=1;

while(i--)a*=v;

return a;

}

（2）实验结果

2、最小二乘拟合

（1）程序设计

#include

#include

void main()

{

int i,n;

float a[2];

float x[15]={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8},z[15];

float

y[15]={33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,301.40,310. 40,318.15,325.15};

for(n=0;n<=14;n++) //增加了数组z；

{z[n]=log(y[n]/x[n]);}

void Approx(float[],float[],int,int,float[]);

Approx(x,z,15,1,a); //变成一次拟合；

//for(i=0;i<=1;i++)

//printf("a[%d]=%f\n",i,a[i]);

printf("a=exp(a[0])=%f\n",exp(a[0]));

printf("b=-a[1]=%f\n",-a[1]);

}

void Approx(float x[],float y[],int m,int n,float a[])

{

int i,j,t;

float *c=new float[(n+1)*(n+2)];

float power(int,float);

void ColPivot(float *,int,float[]);

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

{

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

{

*(c+i*(n+2)+j)=0;

for(t=0;t<=m-1;t++)

*(c+i*(n+2)+j)+=power(i+j,x[t]);

}

*(c+i*(n+2)+n+1)=0;

for(j=0;j<=m-1;j++)

*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);

}

ColPivot(c,n+1,a);

delete c;

}

void ColPivot(float *c,int n,float x[])

{

int i,j,t,k;

float p;

for(i=0;i<=n-2;i++)

{

k=i;

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

if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))

k=j;

if(k!=i)

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

{

p=*(c+i*(n+1)+j);

*(c+i*(n+1)+j)=*(c+k*(n+1)+j);

*(c+k*(n+1)+j)=p;

}

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

{

p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));

for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));

}

}

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

{

for(j=n-1;j>=i+1;j--)

(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));

x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));

}

}

float power(int i,float v)

{

float a=1;

while(i--)a*=v;

return a;

}

（2）实验结果

四、总结

通过曲线拟合，最小二乘法的基本原理的学习，我学会了利用计算机解决现实实际问题中的曲线拟合。

6 数值积分

一、实验目的

（1）通过实际计算体会各种方法的精确度；

（2）会编写用龙贝格算法求定积分的程序。

二、实验内容

编写复化柯特斯求积分公式，并计算例题1例题2，观察n为多少时有6位有效数字。

三、程序设计、实验结果及分析

（1）程序设计

#include

#include

float Cotes(float(*f)(float),float a,float b,int n)

{

int k;

float c,c1=0,c2,c3,c4;

float h=(b-a)/n;

c2=(*f)(a+h/4);

c3=(*f)(a+h/2);

c4=(*f)(a+3*h/4);