计算方法 实验 线性方程组的数值解法.docx

计算方法 实验 线性方程组的数值解法.docx

《计算方法 实验 线性方程组的数值解法.docx》由会员分享，可在线阅读，更多相关《计算方法 实验 线性方程组的数值解法.docx（17页珍藏版）》请在冰豆网上搜索。

计算方法实验线性方程组的数值解法

姓名：

学号：

班级：

线性方程组的数值解法

实习目的

（1）通过实习进一步掌握高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的基本思想；

（2）通过实习进一步掌握高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的计算步骤，并能灵活应用；

（3）通过对高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的调试练习，进一步体会各种算法的特点；

（4）通过对上机调试运行，组不培养解决实际问题的编程能力。

实习要求

（1）熟悉TurboC的编译环境；

（2）实习前复习高斯消去法、列主元高斯消去法、直接三角分解法、雅克比迭代法、高斯-赛德尔迭代法以及追赶法的基本思想和过程；

（3）实习前复习高斯消去法、列主元高斯消去法、直接三角分解法、雅克比迭代法、高斯-赛德尔迭代法以及追赶法的计算步骤。

实习设备

（1）硬件设备：

单机或网络环境下的微型计算机一台；

（2）软件设备：

DOS3.3以上操作系统，TurboC2.0编译器。

实习内容

实习一高斯消去法

（1）用高斯消去法求解线性方程组：

（2）要求：

①将线性方程组写成用矩阵表示的形式，即Ax=b的形式。

②输出系数矩阵的原始元素和经高斯消去法消去后的矩阵元素。

③经高斯消去法消去后的矩阵是一个什么形式的矩阵？

④请写出程序的运行结果。

#include

#defineN4

intGauss（floata[N][N],floatb[N]）

{

inti,j,k,flag=1;

floatt;

for（i=0;i

{

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

{

flag=0;

break;

}

else

{

for（j=i+1;j

{

t=-a[j][i]/a[i][i];

b[j]=b[j]+t*b[i];

for（k=i;k

a[j][k]=a[j][k]+t*a[i][k];

}

}

}

return（flag）;

}

voidzg_matric（floata[N][N],floatb[N]）

{

inti,j;

for（i=0;i

{

for（j=0;j

printf（"%10f",a[i][j]）;

printf（"%10f",b[i]）;

printf（"\n"）;

}

printf（"\n"）;

}

voidmain（）

{

staticfloata[N][N]={{1,-4,6,-5},{-5,21,-33,32},{6,-26,43,-48},{5,-24,45,-64}};

floatb[N]={8,-32,25,-10};

floatx[N]={0,0,0,0};

inti,j,flag;

zg_matric（a,b）;

flag=Gauss（a,b）;

if（flag==0）

printf（"Gaussmethoddoesnotrun."）;

else

{

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

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

{

x[i]=b[i];

for（j=i+1;j

x[i]=x[i]-a[i][j]*x[i];

x[i]=x[i]/a[i][i];

}

for（i=0;i

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

}

}

（3）思考题：

①高斯消去法的基本思想是什么？

答：

高斯消去法是最古老的求解线性代数方程组的方法之一，是消去法的一种特殊形式。

②高斯消去法由哪两个过程组成？

答：

高斯消去法包括消元与回代两个过程。

实习二列主元高斯消去法求解线性方程组

（1）用列主元高斯消去法求解线性方程组：

（2）要求：

①输出系数矩阵的原始元素和经列主高斯消去法消去后的矩阵元素。

②经列主高斯消去法消去后的矩阵是一个什么形式的矩阵？

③请写出输出后的运行结果。

④若线性方程组的准确解为：

，请通过计算绝对误差证明，在一般情况下用经列主高斯消去法求解线性方程组是稳定的。

#include

#include

#defineN3

voidColGauss（floata[N][N],floatb[N]）

{

floatt,max_fab;

inti,j,k,l;

for（i=0;i

{

l=i;

max_fab=fabs（a[i][i]）;

for（j=i+1;j

if（fabs（a[j][i]>max_fab））

{

l=j;

max_fab=fabs（a[l][i]）;

}

if（i<1）

{

t=b[i];

b[i]=b[k];

b[k]=t;

for（j=i;j

{

t=a[i][j];a[i][j]=a[l][j];a[l][j]=t;

}

}

for（j=i+1;j

{

t=-a[j][i]/a[i][i];

b[j]=b[j]+t*b[i];

for（k=i;k

a[j][k]=a[j][k]+t*a[i][k];

}

}

}

voidzg_matric（floata[N][N],floatb[N]）

{

inti,j;

for（i=0;i

{

for（j=0;j

printf（"%10f",a[i][j]）;

printf（"%10f",b[i]）;

printf（"\n"）;

}

printf（"\n"）;

}

voidmain（）

{

staticfloata[N][N]={{0.02,4.0,1.38},{2.0,3.32,0.575},{4.0,5.0,4.0}};

floatb[N]={9.0,16.0,1.0};

floatx[N]={0,0,0};

inti,j;

zg_matric（a,b）;

ColGauss（a,b）;

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

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

{

x[i]=b[i];

for（j=i+1;j

x[i]=x[i]-a[i][j]*x[j];

x[i]=x[i]/a[i][i];

}

for（i=0;i

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

}

（3）思考题：

①列主高斯消去法的基本思想是什么？

答：

列主高斯消去法在高斯消去法每次进行消元之前先选取绝对值最大的元素作为主元素（位于主对角线上的在消去过程中作除数的元素称为主元素，简称主元）进行消去。

②在高斯消去法的基础上，引进经列主高斯消去法的两个主要原因是什么？

答：

当

时，消元过程就无法进行。

另外，即使

，但其绝对值很小时，用它做除数，根据数值运算中“用绝对值很小的数做除数，舍入误差会增大，而且严重影响计算结果的精度”的原则，这种方法在一定程度上也具有局限性。

为了克服这一局限性，产生了列主高斯消元法。

实习三柯朗分解法

（1）用柯朗分解法求解线性方程组：

（2）要求：

①输出系数矩阵的原始元素和经柯朗分解法后的L、U矩阵元素。

②经柯朗分解法后的L、U矩阵各是一个什么形式的矩阵？

③请写出程序的运行结果。

#include

#defineN4

voidLu（floata[N][N],floatL[N][N],floatU[N][N]）

{

inti,j,k;

for（k=0;k

{

for（i=k;i

{

L[i][k]=a[i][k];

for（j=0;j<=k-1;j++）

L[i][k]-=（L[i][j]*U[j][k]）;

}

for（j=k+1;j

{

U[k][j]=a[k][j];

for（i=0;i<=k-1;i++）

U[k][j]-=（L[k][i]*U[i][j]）;

U[k][j]/=L[k][k];

}

}

}

voidsolve（floatb[N],floatL[N][N],floatU[N][N],floatx[N]）

{

intj,k;

floaty[N];

for（k=0;k

{

y[k]=b[k];

for（j=0;j<=k-1;j++）

y[k]=y[k]-L[k][j]*y[j];

y[k]=y[k]/L[k][k];

}

for（k=N-1;k>=0;k--）

{

x[k]=y[k];

for（j=k+1;j

x[k]=x[k]-U[k][j]*x[j];

}

}

voidzg_matric（floata[N][N],floatb[N]）

{

inti,j;

for（i=0;i

{

for（j=0;j

printf（"%10f",a[i][j]）;

printf（"%10f",b[i]）;

printf（"\n"）;

}

printf（"\n"）;

}

voidmain（）

{

staticfloata[N][N]={{1,2,3,4},{1,3,4,5},{2,1,4,4},{2,3,2,5}};

floatb[N]={2,3,3,2};

floatx[N]={0,0,0,0};

floatU[N][N]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};

floatL[N][N];

inti;

zg_matric（a,b）;

Lu（a,L,U）;

solve（b,L,U,x）;

for（i=0;i

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

}

（4）思考题：

①柯朗分解法的基本思想是什么？

答：

②在柯朗分解法汇中使用主对角线上的元素做分母计算时，如果分母绝对值很小（甚至可能为零），仍用它做除数，就会导致舍入误差增大（甚至不能保证柯朗分解法顺利进行），严重影响计算结果的精度。

为解决这一问题，可以采取什么措施？

答：

实习四追赶法

（1）用追赶法求解线性方程组Ax=b：

（2）要求：

请写出程序的运行结果。

#include

#defineN4

voidzhuigan（floata[N][N],floatb[N],floatr[N],floaty[N]）

{

floatq;

inti;

r[0]=a[1][0]/a[0][0];

y[0]=b[0]/a[0][0];

for（i=1;i

{

q=a[i][i]-r[i-1]*a[i-1][i];

r[i]=a[i+1][i]/q;

y[i]=（b[i]-y[i-1]*a[i-1][i]）/q;

}

y[N-1]=（b[N-1]-y[N-2]*a[N-2][N-1]）/（a[N-1][N-1]-r[N-2]*a[N-2][N-1]）;

}

voidsolve（floatr[N],floaty[N],floatx[N]）

{

intk;

x[N-1]=y[N-1];

for（k=N-2;k>=0;k--）

x[k]=y[k]-r[k]*x[k+1];

}

voidmain（）

{

floata[N][N]={{1,3,0,0},{2,7,3,0},{0,2,7,3},{0,0,2,7}};

floatb[N]={-2,-8,-6,5};

floatr[N],y[N];

floatx[N];

inti;

zhuigan（a,b,r,y）;

solve（r,y,x）;

for（i=0;i

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

}

（3）思考题：

①追赶法的基本思想是什么？

答：

追赶法是高斯消去法的一种简化形式，同样分为消元和回代两个过程。

②追赶法适合怎样形式的问题？

答：

系数矩阵为三对角矩阵的方程组解法。

实习五雅克比迭代法与高斯-赛德尔迭代法

（1）用雅克比迭代法与高斯-赛德尔迭代法求解线性方程组Ax=b：

（2）要求：

①请写出程序的运行结果。

②请写出迭代次数。

雅克比迭代法：

#include

#include

#defineN3

#defineMAX_N100

#defineeps1e-6

intyacobi（floata[N][N],floatb[N],floatx[N]）

{

floatd,s,max;

floaty[N];

inti,j,k,flag;

k=0;

for（i=0;i

x[i]=0.0;

while

（1）

{

max=0.0;

k++;

for（i=0;i

{

s=0.0;

for（j=0;j

{

if（j==i）continue;

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

}

y[i]=（b[i]-s）/a[i][i];

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

if（max

}

if（max

{

flag=1;

break;

}

if（k>=MAX_N）

{

flag=0;

break;

}

for（i=0;i

x[i]=y[i];

}

return（flag）;

}

voidzg_matric（floata[N][N],floatb[N]）

{

inti,j;

for（i=0;i

{

for（j=0;j

printf（"%10f",a[i][j]）;

printf（"%12f",b[i]）;

printf（"\n"）;

}

printf（"\n"）;

}

voidmain（）

{

floata[N][N]={{6,2,-1},{1,4,-2},{-3,1,4}};

floatb[N]={-3,2,4};

floatx[N];

inti,k;

zg_matric（a,b）;

k=yacobi（a,b,x）;

if（k==1）

for（i=0;i

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

else

printf（"TheMethodisdisconvergent!

\n"）;

}

高斯-赛德尔迭代法：

#include

#include

#defineN3

#defineMAX_N100

#defineeps1e-6

intseidel（floata[N][N],floatb[N],floatx[N]）

{

floatd,s,max,temp;

inti,j,k,flag;

k=0;

for（i=0;i

x[i]=0.0;

while

（1）

{

max=0.0;

k++;

for（i=0;i

{

s=0.0;

temp=x[i];

for（j=0;j

{

if（j==i）continue;

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

}

x[i]=（b[i]-s）/a[i][i];

d=fabs（x[i]-temp）;

if（max

}

if（max

{

flag=1;

break;

}

if（k>=MAX_N）

{

flag=0;

break;

}

}

return（flag）;

}

voidzg_matric（floata[N][N],floatb[N]）

{

inti,j;

for（i=0;i

{

for（j=0;j

printf（"%10f",a[i][j]）;

printf（"%12f",b[i]）;

printf（"\n"）;

}

printf（"\n"）;

}

voidmain（）

{

floata[][N]={{6,2,-1},{1,4,-2},{-3,1,4}};

floatb[N]={100,-200,345};

floatx[N];

inti,k;

zg_matric（a,b）;

k=seidel（a,b,x）;

if（k==1）

for（i=0;i

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

else

printf（"TheMethodisdisconvergent!

\n"）;

}

（3）思考题：

①雅克比迭代法与高斯-赛德尔迭代法的基本思想是什么？

答：

②观察右端项对迭代收敛是否有影响。

答：