计算方法实验报告.docx

计算方法实验报告.docx

《计算方法实验报告.docx》由会员分享，可在线阅读，更多相关《计算方法实验报告.docx（20页珍藏版）》请在冰豆网上搜索。

计算方法实验报告

数值计算方法上机实验报告

班级：

姓名：

学号：

实验一舍入误差与数值稳定性

一、目的与要求：

1、通过上机编程，复习巩固以前所学程序设计语言；

2、通过上机计算，了解舍入误差所引起的数值不稳定性。

3、通过上机计算，了解运算次序对计算结果的影响，从而尽量避免大数吃小数的现象。

二、实验内容：

●通过正反两个实例的计算，了解利用计算机进行数值计算中舍入误差所引起的数值不稳定性，深入理解初始小的舍入误差可能造成误差积累从而对计算结果的巨大影响。

●通过实际编程，了解运算次序对计算结果的影响，了解实数运算符合的结合律和分配律在计算机里不一定成立。

三、程序与实例

1对n=0,1,2,…,20计算定积分

=

算法1利用递推公式

=

-5

n=1,2,…,20

取

ln6-ln5

0.182322

算法2利用递推公式

n=20,19,…,1

注意到

取

算法：

#include

#include

voidmain（）

{inti,j,k;

longdoublea[21],b[21];

a[20]=0.0087301587;

b[0]=0.18232155;

for（k=0;k<20;k++）

{b[k+1]=-5*b[k]+1.0/（k+1）;

}

for（i=20;i>0;i--）

{a[i-1]=-1.0/5*a[i]+1.0/（5*i）;

}

printf（"nI（A）nI（B）\n"）;

for（j=0;j<21;j++）

printf（"%2d%1.8lf%2d%1.8lf\n",j,b[j],j,a[j]）;

}

运行结果：

说明：

从计算结果可以看出，算法一是不稳定的，而算法二是稳定的。

实验二拉格朗日插值与牛顿插值

一、目的与要求：

熟悉拉格朗日插值多项式和牛顿插值多项式，注意其不同特点；

二、实验内容：

通过拉格朗日插值多项式的实例计算，了解这种求解方法，分析其的优缺点。

三、程序与实例

算法

1．输入x

y

（i=0,1,2,,n）,令L（x

）=0;

2．对=0,1,2,,n计算

l

l

（x）=

L

L

+l

（x）y

程序与实例

例1给定函数表

x

0

1

2

3

4

5

y

-7

-4

5

26

65

128

试用Lagrange插值法求一个三次插值多项式

并由此求

的近似值。

拉格朗日插值：

#include

voidmain（）

{

inti,n;

floatx[100],fx[100];

floatxk,px,fz,fm;

cout<<"输入数据个数：

";

cin>>n;

for（i=0;i

{

cout<<"输入第"<

";

cin>>x[i]>>fx[i];

}

cout<<"已知的数据是:

";

for（i=0;i

cout<<"（"<

cout<

intflag=1;

while（flag!

=0）

{

cout<<"输入所求f（x）的x值：

";

cin>>xk;

px=0;

for（i=0;i

{

fz=1;

fm=1;

for（intj=0;j

if（j!

=i）

{

fz=fz*（xk-x[j]）;

fm=fm*（x[i]-x[j]）;

}

px=（px+（fx[i]*fz）/fm）;

}

cout<<"所求f（x）="<

cout<<"继续请按1,退出请按0:

";

cin>>flag;

}

}

运行结果：

实验三常微分方程数值解法

一、目的与要求：

熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉法

会编制上述方法的计算程序

针对实习题编制程序，并上机计算其所需要的结果；

二、实验内容：

熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉法,体会其解法的功能。

三、程序与实例

例1设初值问题

，取

，试用Euler方法、后退的Euler方法和梯形公式求解。

程序源代码：

#include

#include

#defineN6

voidmain（）

{inti;

doublex=0.0;

doublea[N],b[N],c[N],d[N],X[N];

a[0]=1;b[0]=1;c[0]=1;d[0]=1;X[0]=0.0;

for（i=0;i<5;i++）

{a[i+1]=0.9*a[i]+0.1*x+0.1;

b[i+1]=（b[i]+0.1*x+0.11）/1.1;

c[i+1]=（0.95*c[i]+0.1*x+0.105）/1.05;

d[i+1]=x+exp（x）;

x=x+0.1;

X[i+1]=x;

}

printf（"xEuler方法后退Euler方法梯形公式准确解\n"）;

for（i=0;i

printf（"%0.1f%1.6f%1.6f%1.6f%1.6f\n",X[i],a[i],b[i],c[i],d[i]）;

}

运行结果：

实验四方程求根

一、目的与要求：

通过对二分法和牛顿迭代法作编程练习和上机运算，进一步体会它们在方程求根中的不同特点；

比较二者的计算速度和计算精度。

二、实验内容：

通过对二分法和牛顿迭代法作编程练习和上机运算，进一步体会它们在方程求根中的不同特点

三、程序与实例

例1：

用二分法求方程

在区间

的根。

要求误差不超过0.01。

程序源代码：

#include"stdio.h"

#include"math.h"

floata[100],b[100],x[100],f[100];

floatfunction（floatx）

{

floatf;

f=x*x*x+x*x-3*x-3;

returnf;

}

voidmain（）

{inti=0,j,k=0,l=0,n=0;

floatx1,x2,x0,fx1,fx2,fx0;

x1=1;x2=2;

a[i++]=x1;b[l++]=x2;

fx1=function（x1）;

fx2=function（x2）;

do

{

x0=（x1+x2）/2.0;

x[n++]=x0;

fx0=function（x0）;

f[k]=fx0;

if（fx0*fx1<0）

{

x2=x0;

fx2=fx0;

f[k++]=fx2;

}

else

{

x1=x0;

fx1=fx0;

f[k++]=fx1;

}

b[l++]=x2;

a[i++]=x1;

}while（fabs（x2-x1）>0.01）;

printf（"kabxfx\n"）;

for（j=0;j<7;j++）

printf（"%d%f%f%f%f\n",j,a[j],b[j],x[j],f[j]）;

printf（"Therootis%f\n",x0）;

}

运行结果：

实验五：

高斯消去法求解线性方程组

一、目的与要求：

通过对高斯消去法编程练习和上机运算，进一步体会它们在求解线性方程组中的特点，求出精确解。

二、实验内容：

用高斯消去法求出线性方程组的精确解。

三、程序与实例：

2.给定线性方程组

=

，已知精确解

=

。

程序源代码：

#include

#include

voidmain（）

{

inti,j,k,n,l,h,flag;

floata[100][100],b[100],m[100][100],x[100],y[100],z[100],c[100],d,t;

printf（"Inputn:

\n"）;

scanf（"%d",&n）;

printf（"inputa[i][j]:

\n"）;

for（i=0;i

{

for（j=0;j

{

scanf（"%f",&a[i][j]）;

}

}

printf（"inputb[i]:

\n"）;

for（i=0;i

{

scanf（"%f",&b[i]）;

}

/*******运算过程*******/

flag=1;

for（k=0;k

{

/*---选主元素---*/

d=a[k][k];

h=k;

for（l=k;l

{

if（fabs（a[l][k]）>fabs（d））

{

d=a[l][k];

h=l;

}

}

if（h!

=k）

{

for（j=k;j

{

t=a[h][j];a[h][j]=a[k][j];a[k][j]=t;

}

t=b[k];b[k]=b[h];b[h]=t;

}

if（a[k][k]==0）

flag=0;

/*---具体运算---*/

for（i=k+1;i

{

if（flag==0）

break;

m[i][k]=a[i][k]/a[k][k];

for（j=k;j

{

a[i][j]=a[i][j]-m[i][k]*a[k][j];

}

b[i]=b[i]-m[i][k]*b[k];

}

}

/*******回代过程*******/

if（flag!

=0）

{

x[n-1]=b[n-1]/a[n-1][n-1];

for（i=n-2;（i>0）||（i==0）;i--）

{

y[i]=0;

for（j=i+1;j

{

y[j]=y[j-1]+a[i][j]*x[j];

}

x[i]=（b[i]-y[j-1]）/a[i][i];

}

/*******输出结果*******/

for（i=0;i

{

printf（"x[%d]=%f\n",i,x[i]）;

}

}

else

printf（"ohmygod!

\n"）;

}

运行结果：

实验六解线性方程组的迭代法

一、目的与要求：

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

会编制雅可比及高斯—塞德尔迭代法的程序；

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

二、实验内容：

会编制雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例

例1用迭代法解方程组

=

1.JACOBI迭代法：

#include

#include

intfunction（floaty[3],floatx[3]）;/*判断是否收敛，满足精度函数申明*/

floatx[3]={0,0,0},z;/*定义初始向量x*/

inti,j,k,n=3;

main（）

{

floata[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}},b[3]={7.2,8.3,4.2};

floaty[3],sum;

intflag;

for（k=0;k<100;k++）/*迭代的次数*/

{

for（i=0;i

{

sum=0;

for（j=0;j

{

if（j!

=i）

sum=sum+a[i][j]*x[j];

}

y[i]=（b[i]-sum）/a[i][i];/*算出该迭代时的y[i]*/

}

for（i=0;i

{

printf（"x%d=%-10.6f",i+1,y[i]）;

}

printf（"\n"）;

flag=function（y,x）;/*调用函数function（）*/

if（flag==1）/*结束循环*/

break;

}

}

intfunction（floaty[3],floatx[3]）/*判断是否收敛，满足精度函数的定义*/

{

intflag=0;/*标志主函数中的循环是否要结束*/

z=fabs（y[0]-x[0]）;

for（i=0;i

if（z

z=fabs（y[i]-x[i]）;

if（z<10e-6）

{

flag=1;

printf（"迭代次数k=%d\n",k+1）;/*输出得到最后结果迭代的次数*/

printf（"最后结果:

\n"）;

for（i=0;i

printf（"x%d=%-10.6f",i+1,y[i]）;/*输出方程组的解*/

printf（"\n"）;

}

else

for（i=0;i

x[i]=y[i];

return（flag）;

}

运行结果：

2.GAUSS-SEIDEL迭代法源程序：

#include

#include

#include

#include

#defineMAX_n100

#definePRECISION0.0000001

#defineMAX_Number1000

voidVectorInput（floatx[],intn）

{inti;

for（i=1;i<=n;++i）

{printf（"x[%d]=",i）;

scanf（"%f",&x[i]）;

}

}

voidMatrixInput（floatA[][MAX_n],intm,intn）

{inti,j;

printf（"\n===BegininputMatrixelements===\n"）;

for（i=1;i<=m;++i）

{printf（"Input_Line%d:

",i）;

for（j=1;j<=n;++j）

scanf（"%f",&A[i][j]）;

}

}

voidVectorOutput（floatx[],intn）

{inti;

for（i=1;i<=n;++i）

printf（"x[%d]=%.9f",i,x[i]）;

}

intIsSatisfyPricision（floatx1[],floatx2[],intn）

{inti;

for（i=1;i<=n;++i）

if（fabs（x1[i]-x2[i]）>PRECISION）return1;

return0;

}

intJacobi_（floatA[][MAX_n],floatx[],intn）

{floatx_former[MAX_n];

inti,j,k;

printf（"\nInputvectorx0:

\n"）;

VectorInput（x,n）;

k=0;

do{

for（i=1;i<=n;++i）

{printf（"x[%d]=%.9f",i,x[i]）;

x_former[i]=x[i];

}

printf（"\n\n"）;

for（i=1;i<=n;++i）

{

x[i]=A[i][n+1];

for（j=1;j<=n;++j）

if（j!

=i）x[i]-=A[i][j]*x[j];//Onlydefference:

herex[j]replacex_former[j]

if（fabs（A[i][i]）>PRECISION）

x[i]/=A[i][i];

else

return1;

}

++k;

}while（IsSatisfyPricision（x,x_former,n）&&k

if（k>=MAX_Number）

return1;

else

{printf（"\nG-S%dtimes!

",k）;

return0;

}

}

voidmain（）

{intn;

floatA[MAX_n][MAX_n],x[MAX_n];

printf（"\nInputn="）;

scanf（"%d",&n）;

if（n>=MAX_n-1）

{printf（"\n\007nmust<%d!

",MAX_n）;

exit（0）;

}

MatrixInput（A,n,n+1）;

if（Jacobi_（A,x,n））

printf（"\nG-SFailed!

"）;

else

{printf（"\nOutputSolution:

"）;

VectorOutput（x,n）;

}

printf（"\n\n\007Pressanykeytoquit!

\n"）;

getch（）;

}

运行结果：