计算方法实验指导书.docx
《计算方法实验指导书.docx》由会员分享，可在线阅读，更多相关《计算方法实验指导书.docx（11页珍藏版）》请在冰点文库上搜索。
计算方法实验指导书
《计算方法》实验指导书
自动化学院
自动化实验室
实验一全选主元高斯消去法求解线性方程组
一、实验目的：
用全选主元高斯消去法求解n阶线性方程组
其中
，
，
二、实验方法说明：
用全选主元高斯消去法分两步进行。
第一步：
消去过程
对于k从0开始，一直到n-2作以下三步：
从系数矩阵A的第k行、第k列开始的子阵中选取绝对值最大的元素，并将它交换到主元素的位置上。
归一化。
即
，
消去。
即
，
，
第二步：
回代过程
，
三、实验内容
已知4阶方程组的系数矩阵为
，
求该方程组的解。
实验二一般实矩阵的QR分解
一、实验目的：
用Gram-Schmidt正交化方法对n阶方阵A进行QR分解，其中：
二、实验方法说明：
若：
=
，
=
，…，
=
是n维空间中n个线性无关的相量，则有：
这时矩阵A可以分解为正交矩阵Q和上三角矩阵R的乘积。
其中：
三、实验内容
用Gram-Schmidt正交化方法对
进行QR分解。
参考程序：
1、全选主元高斯消去法求解线性方程组：
#include"stdlib.h"
#include"stdio.h"
#include"math.h"
intcagaus（a,b,n,x）
intn;
doublea[],b[],x[];
{int*js,l,k,i,j,is,p,q;
doubled,t;
js=malloc（n*sizeof（int））;
l=1;
for（k=0;k<=n-2;k++）
{d=0.0;
for（i=k;i<=n-1;i++）
for（j=k;j<=n-1;j++）
{t=fabs（a[i*n+j]）;
if（t>d）{d=t;js[k]=j;is=i;}
}
if（d+1.0==1.0）l=0;
else
{if（js[k]!
=k）
for（i=0;i<=n-1;i++）
{p=i*n+k;q=i*n+js[k];
t=a[p];a[p]=a[q];a[q]=t;
}
if（is!
=k）
{for（j=k;j<=n-1;j++）
{p=k*n+j;q=is*n+j;
t=a[p];a[p]=a[q];a[q]=t;
}
t=b[k];b[k]=b[is];b[is]=t;
}
}
if（l==0）
{free（js）;printf（"fail\n"）;
return（0）;
}
d=a[k*n+k];
for（j=k+1;j<=n-1;j++）
{p=k*n+j;a[p]=a[p]/d;}
b[k]=b[k]/d;
for（i=k+1;i<=n-1;i++）
{for（j=k+1;j<=n-1;j++）
{p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[（n-1）*n+n-1];
if（fabs（d）+1.0==1.0）
{free（js）;printf（"fail\n"）;
return（0）;
}
x[n-1]=b[n-1]/d;
for（i=n-2;i>=0;i--）
{t=0.0;
for（j=i+1;j<=n-1;j++）
t=t+a[i*n+j]*x[j];
x[i]=b[i]-t;
}
js[n-1]=n-1;
for（k=n-1;k>=0;k--）
if（js[k]!
=k）
{t=x[k];x[k]=x[js[k]];x[js[k]]=t;}
free（js）;
return
（1）;
}
main（）
{inti;
doublea[4][4]={{0.2368,0.2471,0.2568,1.2671},
{0.1968,0.2071,1.2168,0.2271},
{0.1581,1.1675,0.1768,0.1871},
{1.1161,0.1254,0.1397,0.1490}};
doublex[4],b[4]={1.8471,1.7471,1.6471,1.5471};
i=cagaus（a,b,4,x）;
if（i!
=0）
for（i=0;i<=3;i++）
printf（"x（%d）=%e\n",i,x[i]）;
}
参数说明：
a----双精度实型二维数组，体积为
。
存放方程组的系数矩阵，在本函数中要被破坏。
b-----双精度实型一维数组，体积为
。
存放方程组右端的常数向量，在本函数中要被破坏。
n----整型变量。
方程的阶数。
x----双精度实型一维数组，体积为
。
返回方程组的解向量。
2、用Gram-Schmidt正交化法进行QR分解
#include"stdio.h"
#include"math.h"
voidschmidt（a,b,abj,n）
doublea[],b[],abj[];
intn;
{
inti,j,k;
double*bpm,*bp;
bp=malloc（（n*n）*sizeof（double））;
bpm=malloc（n*sizeof（double））;
for（i=0;i<=n-1;i++）
{
for（j=0;j<=n-1;j++）
{bp[i*n+j]=a[n*i+j];
abj[i*n+j]=0.0;
}
bpm[i]=0.0;
}
for（i=0;i<=n-1;i++）
{bpm[0]=bpm[0]+bp[i*n]*bp[i*n];}
for（i=0;i<=n-1;i++）
{b[i*n]=bp[i*n]/sqrt（bpm[0]）;}
for（i=1;i<=n-1;i++）
{
for（j=0;j<=i-1;j++）
{
for（k=0;k<=n-1;k++）
{abj[j*n+i]=a[k*n+i]*b[k*n+j]+abj[j*n+i];}
for（k=0;k<=n-1;k++）
{bp[k*n+i]=bp[k*n+i]-b[k*n+j]*abj[j*n+i];}
}
for（j=0;j<=n-1;j++）
{bpm[i]=bpm[i]+bp[j*n+i]*bp[j*n+i];}
for（j=0;j<=n-1;j++）
{b[j*n+i]=bp[j*n+i]/sqrt（bpm[i]）;}
}
for（i=0;i<=n-1;i++）
{abj[i*n+i]=sqrt（bpm[i]）;}
return;
}
main（）
{inti,j,k;
doublea[3][3]={{0,1,1},{1,1,0},{1,0,1}},b[3][3],abj[3][3];
schmidt（a,b,abj,3）;
printf（"MatQis:
\n"）;
for（i=0;i<=2;i++）
{for（j=0;j<=2;j++）
{printf（"%e",b[i][j]）;}
printf（"\n"）;
}
printf（"\n"）;
printf（"MatRis:
\n"）;
for（i=0;i<=2;i++）
{for（j=0;j<=2;j++）
{printf（"%e",abj[i][j]）;}
printf（"\n"）;
}
}
参数说明：
a----双精度实型二维数组，体积为
。
存放A矩阵。
b----双精度实型二维数组，体积为
。
返回Q矩阵。
abj----双精度实型二维数组，体积为
。
返回R矩阵。
n----A的阶数。
升级会员 