数值分析实验报告.docx

数值分析实验报告.docx

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

数值分析实验报告

《数值分析》

实验报告

学院：

计算机科学与软件学院

姓名：

XXX

班级：

计算机XX班

学号：

XXXXXX

实验一：

舍入误差与数值稳定性

实验目的：

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

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

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

实验内容：

用两种不同的顺序计算

，分析其误差的变化。

实验流程图：

实验源程序:

#include

#include

voidmain（）

{inti;

floats1=0,s2=0,d1,d2;

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

s1=s1+1.0f/（i*i）;

for（i=10000;i>=1;i--）

s2=s2+1.0f/（i*i）;

d1=（float）（fabs（1.644834-s1））;

d2=（float）（fabs（1.644834-s2））;

printf（"正向求和结果为%f\n误差为%f\n\n",s1,d1）;

printf（"反向求和结果为%f\n误差为%f\n\n",s2,d2）;

if（d1

printf（"正向求和误差小于负向求和误差\n"）;

elseif（d1==d2）

printf（"正向求和误差等于负向求和误差\n"）;

else

printf（"正向求和误差大于负向求和误差\n"）;

}

实验结果：

实验分析：

第一次做数值实验，又一次使用C语言编程，没有了刚学习C语言的艰难，能够将实验步骤转换成流程图并编写出完整的实验代码，在经过多次调试、改正后得到正确的程序和结果。

这个实验较简单，计算误差时如果输入数据有误差，而在计算过程中舍入误差不增长，则称此算法是稳定的，否则称此算法是数值不稳定的，减少运算次数可以减小舍入误差。

在运算中，如果参加运算的数的数量级相差很大，而计算机位数有限，如不注意运算次序就可能出现大数“吃掉”小数的现象，进而影响计算结果的可靠性，所以计算过程中要注意运算次序，避免出现这种现象。

实验二：

拉格朗日插值法和牛顿插值法

实验目的：

分别用拉格朗日差值和牛顿插值解决数学问题，并比较各方法的优略。

1、拉格朗日插值

实验内容：

按下列数据

x

-3.0

-1.0

1.0

2.0

3.0

y

1.0

1.5

2.0

2.0

1.0

作二次插值，并求x

=-2，x

=0，x

=2.75时的函数近似值。

实验流程图：

实验源程序:

#include

floatlagrange（int,float,float[],float[]）;

voidmain（）

{

floata,x[50],y[50],l;

intm,n;

printf（"题目:

按下列数据\n"）;

printf（"x:

-3.0-1.01.02.03.0\n"）;

printf（"y:

1.01.52.02.01.0\n"）;

printf（"作二次插值，并求x=-2,x=0,x=2.75时的函数近似值.\n\n"）;

printf（"输入插值次数:

"）;

scanf（"%d",&n）;

printf（"输入计算次数:

"）;

scanf（"%d",&m）;

for（inti=0;i<=m;i++）

{

for（intj=0;j<=n;j++）

{printf（"输入要计算的第%d个数的第%d个节点及其函数值:

",i+1,j+1）;

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

}

printf（"输入要计算的x的值:

"）;

scanf（"%f",&a）;

l=lagrange（n,a,x,y）;

printf（"%f\n",l）;

}

}

floatlagrange（intn,floata,floatx[],floaty[]）

{

floatl=0,w;

inti,j;

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

{

w=1;

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

{

if（i!

=j）

w=w*（（a-x[j]）/（x[i]-x[j]））;

}

l=l+w*y[i];

}

returnl;

}

实验结果：

2、牛顿插值

实验内容：

按下列数据

x

0.30

0.42

0.50

0.58

0.66

0.72

y

1.04403

1.08462

1.11803

1.15603

1.19817

1.23223

作五次插值，并求x

=0.46，x

=0.55，x

=0.60时的函数近似值。

实验源程序:

#include

#defineM6

floatw（floatX,intn,floatx[]）;

floatquotient（intk,inti,floatf[][M],floatx[],floaty[]）;

floatnewton（floatX,intn,floatf[][M],floatx[],floaty[]）;

voidmain（）

{

floatx[M],y[M],f[M][M];

x[0]=0.30f;y[0]=1.04403f;

x[1]=0.42f;y[1]=1.08462f;

x[2]=0.50f;y[2]=1.11803f;

x[3]=0.58f;y[3]=1.15603f;

x[4]=0.66f;y[4]=1.19817f;

x[5]=0.72f;y[5]=1.23223f;

printf（"x:

0.300.420.500.580.660.72\n"）;

printf（"y:

1.044031.084621.118031.156031.198171.23223\n"）;

printf（"做五次插值,并求x=0.46,x=0.55,x=0.60时的函数近似值.\n\n"）;

for（inti=0;i

f[0][i]=y[i];

floatN[3];

N[0]=newton（0.46f,5,f,x,y）;

N[1]=newton（0.55f,5,f,x,y）;

N[2]=newton（0.60f,5,f,x,y）;

printf（"x=0.46时函数的近似值为%f\n",N[0]）;

printf（"x=0.55时函数的近似值为%f\n",N[1]）;

printf（"x=0.60时函数的近似值为%f\n",N[2]）;

}

floatw（floatX,intn,floatx[]）

{

floatw=1.0;

for（inti=0;i

w=w*（X-x[i]）;

returnw;

}

floatquotient（intk,inti,floatf[][M],floatx[],floaty[]）

{

if（k==0）

f[0][i]=y[i];

else

f[k][i]=（quotient（k-1,i,f,x,y）-quotient（k-1,i-1,f,x,y））/（x[i]-x[i-k]）;

returnf[k][i];

}

floatnewton（floatX,intn,floatf[][M],floatx[],floaty[]）

{

floatN;

N=f[0][0];

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

N=N+w（X,i,x）*quotient（i,i,f,x,y）;

returnN;

}

实验结果：

实验分析：

拉格朗日插值法和牛顿插值法用拉格朗日插值多项式计算函数近似值时，如需增加插值节点，那么原来算出的数据均不能利用，必须重新计算。

用牛顿差商插值多项式中各阶差商用相应查分代替，就可得到各种形式的等距节点插值公式。

牛顿插值法比拉格朗日插值法节省计算量，且便于程序设计，计算增加节点时，计算只需要增加一项，而且牛顿插值更容易计算高次插值。

实验三：

复化积分法

实验目的：

学会用复化积分法提高求积的精度，熟记并掌握复化梯形和复化辛卜生公式。

实验内容：

分别用复化梯形公式和复化辛卜生公式计算f（x）=sin（x）/x的积分，并与准确值比较判断精度。

实验流程图：

实验源程序：

#include

#include

voidmain（）

{

intm,n,k,i;

floatTn,Sn,d1,d2,a=0,b=0,c=0;

floatx[1000],y[1000],z[1000],w[1000];

{

printf（"将复化梯形区间划分:

"）;

scanf（"%d",&n）;

for（i=1;i

{

x[i]=float（1.0/n*i）;

y[i]=（float）（sin（x[i]）/x[i]）;

a=a+y[i];

}

Tn=（float）（（1.0+sin（1.0）/1.0+2*a）/2.0/n）;

printf（"复化梯形输出:

%f\n",Tn）;

printf（"将复化辛卜生区间划分为:

"）;

scanf（"%d",&m）;

}

for（k=1;k

{

z[k]=float（1.0/m*k）;

z[k+1]=float（1.0/m*（k+1））;

z[k+1/2]=float（z[k]+z[k+1]）/2;

w[k]=（float）（sin（z[k]）/z[k]）;

w[k+1/2]=（float）（sin（z[k+1/2]）/（（z[k]+z[k+1]）/2））;

b=b+w[k];

c=c+w[k+1/2];

}

Sn=（float）（（1.0+2*b+4*c+sin（1.0）/1.0）/6.0/m）;

printf（"复化辛卜生输出:

%f\n",Sn）;

{

d1=（float）（fabs（0.9460831-Tn））;

d2=（float）（fabs（0.9460831-Sn））;

printf（"复化梯形误差：

%f\n",d1）;

printf（"复化辛卜生误差：

%f\n",d2）;

if（d1>d2）

printf（"复化梯形求法精度低于复化辛卜生求法\n"）;

elseif（d1==d2）

printf（"复化梯形求法精度等于复化辛卜生求法\n"）;

else（d1

;printf（"复化梯形求法精度高于复化辛卜生求法\n"）;

}

}

实验结果：

实验分析：

许多实际问题常常需要计算积分才能求解，梯形公式、牛顿-柯特斯公式和辛卜生公式是常用的求积公式，为避免可选取适当多的节点，即选取相对高阶的牛顿-柯特斯公式，但由稳定性分析可知，高次插值会产生龙格现象，故采用分段低次插值，然后利用积分的区间可加性得积分，即所谓的分段低次合成的复化求积公式。

复化求积法能改善求积精度，复化的梯形法和辛卜生法当步长h—〉0时，均收敛到所求的积分值。

若将步长h减半（即等分数n加倍），则梯形法、辛卜生法和牛顿-柯特斯法的误差分别减至原有误差的1/4、1/16和1/64。

实验四：

改进欧拉方法解初值问题

实验目的：

学会用改进的欧拉公式求解初值。

实验内容：

用改进欧拉方法解初值问题y’=x+y;y（0）=1。

0

实验流程图：

实验源程序：

#include

#include

#defineMAX50

voidmain（）

{

doublea=0,b=1,h=0.1,y_p,y_c,x[MAX],y[MAX],z[MAX];

intn,i;

n=int（（b-a）/h）;

x[0]=0;

y[0]=1;

z[0]=-x[0]-1+2*exp（x[0]）;

for（i=0;i

{

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

y_p=y[i]+h*（x[i]+y[i]）;

y_c=y[i]+h*（x[i+1]+y_p）;

y[i+1]=（y_p+y_c）/2;

z[i+1]=-x[i+1]-1+2*exp（x[i+1]）;

}

printf（"X的值\t欧拉方法求得的解;\t准确解\t\t\t误差为\n"）;

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

printf（"%.1f\t%f\t\t%f\t\t%f\n",x[i],y[i],z[i],z[i]-y[i]）;

}

实验结果：

实验分析：

这个实验主要利用常微分方程求解，常微分方程的初值问题的数值解法的特点是：

都采用步进式，即求数值解是按节点的顺序逐步推进求得y0,y1,…yn。

这类算法，关键在于建立从已知信息y0,y1，…yi计算yi+1的递推公式。

其中梯形公式虽然提高了精度，但其算法复杂，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数f的值，计算量很大，而且往往难以预测，所以用改进欧拉公式，先用欧拉公式求得一个初步的近似值，再用梯形公式将它校正一次。

这样建立的预测-校正系统通常称为改进的欧拉公式。

设计计算方案有六个环节：

预测、改进、计算、校正、改进、计算。

实验五二分法、简单迭代法、牛顿迭代法

实验目的：

比较迭代法，二分法和牛顿迭代法计算方程的根的准确度。

实验内容：

分别用下列方法求f（x）=x3-3x-1=0在x0=2附近的根。

根的准确值为x*=1.87938524…，要求准确到四位有效数字，并对比各种算法的计算量。

（1）二分法；

（2）简单迭代法；（3）牛顿迭代法

实验流程图：

实验源程序：

#include

#include

voidERFENFA（）

{

doublex0,x1=1.0,x2=2.0,fx0,fx1,fx2;

intcount=0;

fx1=x1*x1*x1-3*x1-1;

fx2=x2*x2*x2-3*x2-1;

do{

x0=（x1+x2）/2;

fx0=x0*x0*x0-3*x0-1;

if（fx0==0）{

printf（"失败"）;break;}

if（fx0*fx1<0）{

x2=x0;fx2=fx0;}

else{x1=x0;fx1=fx0;}

count++;

}

while（fabs（x1-x2）>5e-4）;

printf（"二分法求方程的根为x=%f\n共二分了%d次\n",x0,count）;

}

voidDIEDAIFA（）

{

doublex0,x1=2.0;

intcount=0;

do

{

x0=x1;

x1=exp（log（3*x0+1）/3）;

count++;

}

while（fabs（x1-x0）>5e-4）;

printf（"简单迭代法求方程的根为x=%f\n共迭代了%d次\n",x1,count）;

}

voidNewton（）

{

doublex1,x2=2.0,f1,f2;

intcount=0;

do

{x1=x2;

f1=x1*x1*x1-3*x1-1;

f2=3*x1*x1-3;

x2=x1-f1/f2;

count++;

}

while（fabs（x1-x2）>5e-4）;

printf（"牛顿迭代法求方程的根为x=%f\n共迭代了%d次\n",x2,count）;

}

voidmain（）

{

ERFENFA（）;

DIEDAIFA（）;

Newton（）;

}

实验结果：

实验分析：

二分法、迭代法和牛顿迭代法都是方程求根的常用方法，二分法是逐步搜索方法的一种改进，在有根区间[a，b]上无限地进行二分过程，这些区间最终必将收缩于一点，该点就是所求的根。

迭代法是一种逐次逼近法，其基本思想是将隐式方程归结为一组显式的计算公式，就是说，迭代过程实质上就是一个逐步显式化的过程。

牛顿迭代法可以加速迭代速度。

总的来说，二分法算法简单，易于操作，但是精度较低，计算量大;简单迭代法精度较二分法有所提高但在迭代时需要选择合适的迭代公式，以得到较快的收敛速度；牛顿迭代法精度最高，计算量小，但需要保证在迭代过程中分母不能为零。

实验六高斯列主元消去法和直接三角分解法（LU分解）

实验目的：

了解解线性方程组的直接法，掌握直接三角分解法求解方程组的解。

实验内容：

用直接三角分解法（LU分解）求方程组的解

系数矩阵：

10787常向量：

10

75658

861096

759107

精确解为：

（-60，102，-27，16）

实验流程图：

实验源程序：

#include

#include

#defineMAX10

voidmain（）

{

doubleA[4][5]={{10,7,8,7,10},{7,5,6,5,8},{8,6,10,9,6},{7,5,9,10,7}};

doubles,t,detA,max;

inti,j,k,m,n=4;

detA=1;

//消元计算

for（k=0;k

{

max=fabs（A[k][k]）;

m=k;

for（i=k+1;i

if（fabs（A[i][k]）>max）

{

max=fabs（A[i][k]）;

m=i;

}

if（max==0）

{

detA=0;

break;

}

if（k!

=m）

{

for（i=k;i

{

t=A[k][i];

A[k][i]=A[m][i];

A[m][i]=t;

}

detA=-detA;

}

for（i=k+1;i

{

A[i][k]=A[i][k]/A[k][k];

for（j=k+1;j

A[i][j]=A[i][j]-A[i][k]*A[k][j];

}

detA=detA*A[k][k];

}

detA=detA*A[n-1][n-1];

//回代计算

A[n-1][n]=A[n-1][n]/A[n-1][n-1];

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

{

s=0;

for（j=i+1;j

s=s+A[i][j]*A[j][n];

A[i][n]=（A[i][n]-s）/A[i][i];

}

//输出结果

printf（"用直接三角分解法解的方程组的解为：

\n"）;

for（i=0;i

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

printf（"系数矩阵的行列式为：

%4.1f",detA）;

printf（"\n"）;

}

实验结果：

实验分析：

解线性方程组的直接方法就是经过有限步算术运算即可求得方程组精确解的方法，实际计算中由于舍入误差的存在和影响，这种方法也只能求得近似解，应用高斯消去法及其变形可以解决这个问题，将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到L,U的递推公式，而不需要任何中间步骤，这就是所谓的直接三角分解法，等价于求解方程组

（1）Ly=b，求y；

（2）Ux=y，求x。

L和U的计算公式：

。

U中元素ukj等于对应A中元素akj减去一个内积，此内积等于ukj左边L的同行元素与上边U的同列元素相应的乘积之和，L元素lik的计算与U中元素计算方法相同，但最后还需除以同列U的对角元ukk。

从代数角度看，直接分解法和高斯消去法本质上一样，由直接三角分解法计算公式，可采用“双精度累加”，以提高精度，如果不采用“双精度累加”，那么两种方法将给出相同的结果。