数值代数上机报告.docx
《数值代数上机报告.docx》由会员分享，可在线阅读，更多相关《数值代数上机报告.docx（33页珍藏版）》请在冰点文库上搜索。
数值代数上机报告
Doolittle分解
报告1
一、目的意义：
把矩阵A分解成一个下三角阵与一个上三角阵的乘积，即A=LR，其中L为下三角阵，R为上三角阵，这样原线性方程组就可以化为
的求解问题，方便求解。
二、算法：
1）输入系数矩阵A；
2）利用公式
和
交错进行，计算得出矩阵L和R；
3）回带到
中得出原线性方程组的解。
三、源程序：
#include
#include
#include
#include
#defineN100
main（）
{
inti,j,k,s,n;
printf（"请输入系数矩阵A的阶数：
n="）;
scanf（"%d",&n）;
floata[N][N]={0},L[N][N]={0},R[N][N]={0},sigma1,sigma2,b[N],y[N],x[N];
/*为L主对角线元素赋1*/
for（i=0;i{
L[i][i]=1;
}
printf（"请输入系数矩阵A：
\n"）;/*输入系数矩阵A*/
for（i=0;i{
for（j=0;jscanf（"%f",&a[i][j]）;
}
for（k=0;k{
for（j=k;j{
sigma1=0;
for（s=0;s<=k-1;s++）
sigma1+=L[k][s]*R[s][j];
R[k][j]=a[k][j]-sigma1;
}
for（i=k;i{
sigma2=0;
for（s=0;s<=k-1;s++）
sigma2+=L[i][s]*R[s][k];
L[i][k]=（a[i][k]-sigma2）/R[k][k];
}
}
printf（"\nA矩阵为：
\n"）;/*输出矩阵L、R*/
for（i=0;i{
for（j=0;jprintf（"%5.1f",a[i][j]）;
printf（"\n"）;
}
printf（"\nL矩阵为：
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.1f",L[i][j]）;
printf（"\n"）;
}
printf（"\nR矩阵为：
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.1f",R[i][j]）;
printf（"\n"）;
}
printf（"请输入b矩阵\n",i+1）;
for（i=0;iscanf（"%f",&b[i]）;
for（i=0;i{
sigma1=0;
for（k=0;k<=i-1;k++）
sigma1+=L[i][k]*y[k];
y[i]=b[i]-sigma1;
}
for（i=n-1;i>=0;i--）
{
sigma2=0;
for（k=i+1;ksigma2+=R[i][k]*x[k];
x[i]=（y[i]-sigma2）/R[i][i];
}
printf（"解得x为：
\n"）;
for（i=0;iprintf（"%5.1f",x[i]）;
printf（"\n"）;
}
四、计算结果与分析：
分析：
运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想
五、参考文献：
[1]刑志栋.矩阵数值分析.陕西：
陕西科学技术出版社，2005
[2]谭浩强.C语言程序设计.北京：
清华大学出版社，2005
报告2
cholesky分解
一、目的意义：
对称正定矩阵是实践中经常遇到的一种特殊矩阵类型矩阵，由于矩阵本身的流量好兴致，使得cholesky分解在存储和运算量上较一般消去法节省一半左右，且解的精度高，cholesky分解方法是目前计算机求解该类问题最有效的方法之一。
二、算法：
1）输入系数矩阵A；
2）
利用公式
交错进行，计算得出矩阵L和D；
3）回带到
中得出原线性方程组的解
三、源程序：
#include
#include
#include
#defineEPS1.0e-8
#defineN20
doublea[N][N],b[N],x[N];
intn;
intzhuyuan（introw）;/*选主元*/
voidhangjiaohuan（introw1,introw2）;/*行交换*/
voidxiaoyuan（introw）;/*消元*/
voidhuidai（）;/*回代*/
voidmain（）
{
printf（"请输入方程的维数n:
n="）;
scanf（"%d",&n）;
getchar（）;
printf（"输入%d行%d列系数矩阵A：
\n",n,n）;
for（inti=0;i{
for（intj=0;jscanf（"%lf",&a[i][j]）;
getchar（）;
}
printf（"\n输入线性方程组右端项b[%d]:
\n",n）;
for（i=0;i{
scanf（"%lf",&b[i]）;
}
getchar（）;
for（i=0;i{
doublerowmark=zhuyuan（i）;
if（rowmark==-1）
{
printf（"多解!
"）;
system（"pause"）;
return;
}
if（rowmark!
=i）
{
hangjiaohuan（i,rowmark）;
}
xiaoyuan（i）;
}
huidai（）;
printf（"\n线性方程组的解为:
\n"）;
for（i=0;i{
printf（"x%d=%lf",i+1,x[i]）;
printf（"\n"）;
}
printf（"\n"）;
system（"pause"）;
}
intzhuyuan（introw）
{
doubleelem=a[row][row];
introwmark=row;
for（inti=row+1;i{
if（elem{
elem=a[i][row];
rowmark=i;
}
}
if（fabs（elem）<=EPS）
{
return-1;
}
returnrowmark;
}
voidhangjiaohuan（introw1,introw2）
{
doubletmp;
tmp=b[row1];
b[row1]=b[row2];
b[row2]=tmp;
for（intj=0;j{
tmp=a[row1][j];
a[row1][j]=a[row2][j];
a[row2][j]=tmp;
}
}
voidxiaoyuan（introw）
{
for（inti=row+1;i{
intj=row;
doubletmp=a[i][j]/a[row][row];
b[i]-=tmp*b[row];
for（a[i][j++]=0;j{
a[i][j]-=tmp*a[row][j];
}
}
}
voidhuidai（）
{
x[n-1]=b[n-1]/a[n-1][n-1];
for（inti=n-2;i>=0;i--）
{
doublesum=0.0;
for（intj=i+1;j{
sum-=a[i][j]*x[j];
}
x[i]=（b[i]+sum）/a[i][i];
}
}
四、计算结果与分析：
分析：
运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想
五、参考文献：
[1]刑志栋.矩阵数值分析.陕西：
陕西科学技术出版社，2005
[2]谭浩强.C语言程序设计.北京：
清华大学出版社，2005
Jacobi迭代法
报告3
一、目的意义：
通过有限次的算术运算得到线性方程组的近似值
二、算法：
1）输入系数矩阵A；
2）输入迭代的初始向量；
3）利用公式
，计算得出原线性方程组的近似解。
三、源程序：
#include
#include
#include
#include
#defineEPS1e-6
#defineMAX100
#defineN100
float*Jacobi（floata[N][N],intn）
{
float*x,*y,s;
doubleepsilon;
inti,j,k=0;
x=（float*）malloc（n*sizeof（float））;
y=（float*）malloc（n*sizeof（float））;
printf（"请输入迭代的初始向量：
\n"）;
for（i=0;iscanf（"%f",&x[i]）;
while
（1）
{
epsilon=0;
k++;
for（i=0;i{
s=0;
for（j=0;j{
if（j==i）continue;
s+=a[i][j]*x[j];
}
y[i]=（a[i][n]-s）/a[i][i];
epsilon+=fabs（y[i]-x[i]）;
}
if（epsilon{
printf（"\n迭代次数为：
%d次\n",k）;
returny;
}
if（k>=MAX）
{
printf（"TheMethodisdisconvergent"）;
returny;
}
for（i=0;ix[i]=y[i];
}
}
main（）
{
inti,j,n;
floata[N][N];
float*x;
printf（"请输入系数矩阵的阶数：
n="）;
scanf（"%d",&n）;
printf（"请按行输入线性方程组的增广矩阵：
\n"）;
for（i=0;i{
for（j=0;jscanf（"%f",&a[i][j]）;
}
x=（float*）malloc（3*sizeof（float））;
x=Jacobi（a,n）;
printf（"\n原线性方程组的解为：
\n"）;
for（i=0;iprintf（"x[%d]=%f\n",i,x[i]）;
getch（）;
}
四、计算结果与分析：
分析：
运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想
五、参考文献：
[1]刑志栋.矩阵数值分析.陕西：
陕西科学技术出版社，2005
[2]谭浩强.C语言程序设计.北京：
清华大学出版社，2005
报告4
乘幂法
一、目的意义：
阶数超过四次的多项式的根一般不能用有限次运算求得，因而矩阵特征值的计算方法本质上总是采用迭代格式，而采用乘幂法可以计算最大的一个或按模最大的几个特征值和特征相应特征向量。
二、算法：
1）输入系数矩阵A；
2）取出事向量为
；
3）
利用乘幂法的迭代格式
计算得出原矩阵按模最大特征值和相应特征的向量。
三、源程序：
#include
#include
#defineN100
#defineM100
#defineepsilon0.00001
intmain（）
{
intn;
inti,j,k;
doublexmax,oxmax;
staticdoublea[N][N];
staticdoublex[N],nx[N];
printf（"输入矩阵维数n:
"）;
scanf（"%d",&n）;
if（n>N）
{
printf（"theinpurnislargerthanMAX_N,pleaseredefinetheN.\n"）;
return1;
}
if（n<=0）
{
printf（"pleaseinputanumberbetween1and%d.\n",N）;
return1;
}
//输入A矩阵
printf（"请输入矩阵A（%d,%d）:
\n",n,n）;
for（i=0;ifor（j=0;jscanf（"%lf",&a[i][j]）;
for（i=0;ix[i]=1;
oxmax=0;
for（i=0;i{
for（j=0;j{
nx[j]=0;
for（k=0;knx[j]+=a[j][k]*x[k];
}
xmax=0.0;
for（j=0;jif（fabs（nx[j]）>xmax）
xmax=fabs（nx[j]）;
for（j=0;jnx[j]/=xmax;
for（j=0;jx[j]=nx[j];
if（fabs（xmax-oxmax）{
printf（"按模最大特征值=%.4f\n",xmax）;
printf（"相应的特征向量为:
\n"）;
for（i=0;iprintf（"%.4f\n",nx[i]）;
break;
//return0;
}
oxmax=xmax;
}
return0;
}
四、计算结果与分析：
分析：
运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想
五、参考文献：
[1]刑志栋.矩阵数值分析.陕西：
陕西科学技术出版社，2005
[2]谭浩强.C语言程序设计.北京：
清华大学出版社，2005
报告5
逆幂法
一、目的意义：
逆幂法用于计算非奇异矩阵A的按模最小特征值和相应的特征向量
二、算法：
1）输入系数矩阵A和近似特征值；
2）将乘幂法用于
便可求出A的按模最小的特征值和相应的特征向量。
三、源程序：
#include
#include
#definee0.0001
#defineN100
voidmain（）
{
inti,j,k,n,p,s,kk=1,flag=1;
doubleq,m,lm,Lm,m0;
doubleA[N][N],B[N][N],L[N][N],R[N][N],l[N][N],r[N][N];
doublez0[N],z1[N],y0[N],y1[N],y[N];
printf（"请输入矩阵的阶数：
n="）;
scanf（"%d",&n）;
printf（"\n请输入矩阵A:
\n"）;
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
scanf（"%lf",&A[i][j]）;
printf（"\n请输入近似特征值：
"）;
scanf（"%lf",&lm）;
/*令B=A-lm*I,用Doolittle分解B=LR*/
for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
{
if（j==i）
B[i][j]=A[i][j]-lm;
else
B[i][j]=A[i][j];
}
for（k=1;k<=n;k++）
{
for（i=k;i<=n;i++）
{
for（q=0,p=1;pq=q+L[k][p]*R[p][i];
R[k][i]=B[k][i]-q;
}
for（i=k+1;i<=n;i++）
{
for（q=0,p=1;pq=q+L[i][p]*R[p][k];
L[i][k]=（B[i][k]-q）/R[k][k];
}
}
for（i=1;i<=n;i++）
for（j=i;j<=n;j++）
{
if（i==j）
L[i][j]=1;
else
L[i][j]=0;
}
for（i=2;i<=n;i++）
for（j=1;jR[i][j]=0;
/*第一步迭代，用Doolittle计算R*y1=y0中的y1*/
for（i=1;i<=n;i++）
y0[i]=1;
for（k=1;k<=n;k++）
{
for（i=k;i<=n;i++）
{
for（q=0,p=1;pq=q+l[k][p]*r[p][i];
r[k][i]=R[k][i]-q;
}
for（i=k+1;i<=n;i++）
{
for（q=0,p=1;pq=q+l[i][p]*r[p][k];
l[i][k]=（R[i][k]-q）/r[k][k];
}
}
for（i=1;i<=n;i++）
{
for（q=0,k=1;kq=q+l[i][k]*y[k];
y[i]=y0[i]-q;
}
for（i=n;i>=1;i--）
{
for（q=0,k=i+1;k<=n;k++）
q=q+r[i][k]*y1[k];
y1[i]=（y[i]-q）/r[i][i];
}
/*计算m,z1（存放在下列代码的z0中）*/
for（m=y1[1],i=2;i<=n;i++）
if（fabs（y1[i]）>fabs（y1[i-1]））
m=y1[i];
for（i=1;i<=n;i++）
z0[i]=y1[i]/m;
/*第n（>1）次迭代*/
kk=2;
while（flag==1）
{
/*用Doolittle解方程LR*y2=z1中的y2，其中B=LR*/
for（k=1;k<=n;k++）
{
for（i=k;i<=n;i++）
{
for（q=0,p=1;pq=q+l[k][p]*r[p][i];
r[k][i]=B[k][i]-q;
}
for（i=k+1;i<=n;i++）
{
for（q=0,p=1;pq=q+l[i][p]*r[p][k];
l[i][k]=（B[i][k]-q）/r[k][k];
}
}
for（i=1;i<=n;i++）
{for（q=0,k=1;kq=q+l[i][k]*y[k];
y[i]=z0[i]-q;
}
for（i=n;i>=1;i--）
{for（q=0,k=i+1;k<=n;k++）
q=q+r[i][k]*y1[k];
y1[i]=（y[i]-q）/r[i][i];
}
/*规范化*/
for（m=y1[1],i=2;i<=n;i++）
if（fabs（y1[i]）>fabs（y1[i-1]））m=y1[i];
for（i=1;i<=n;i++）
z1[i]=y1[i]/m;
/*判断：
若|Z[k]-Z[k-1]|>=e，则进行下一轮迭代；否则得解。
*/
for（s=0,i=1;i<=n;i++）
if（（fabs（z0[i]）-fabs（z1[i]））if（s==n）flag=0;
else
{kk++;
for（i=1;i<=n;i++）
{y0[i]=y1[i];
z0[i]=z1[i];
}
}
}
m0=1/m;
Lm=lm+m0;
printf（"最接近的特征值为：
\n"）;
printf（"%f\n\n",Lm）;
printf（"相应的特征向量为：
\n"）;
for（i=1;i<=n;i++）
printf（"%lf,",z1[i]）;
printf（"\n\n"）;
}
四、计算结果与分析：
分析：
运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想
五、参考文献：
[1]刑志栋.矩阵数值分析.陕西：
陕西科学技术出版社，2005
[2]谭浩强.C语言程序设计.北京：
清华大学出版社，2005
报告6
古典Jacobi算法
一、目的意义：
寻找相似变换，是对称矩阵A经过变换之后所得矩阵的非对角线元素的平方和减少，对角线元素的平方和增大，且保持对称性不变，不断的施行这种正交变换，最终是非对角元素的平方和任意的接近与零，对角线元素平方的取极大值。
二、算法：
1）输入系数矩阵A；
2）取
，形成一个相似矩阵序列
，k=1,2,…
3）选取
，使得元素
=0，这时
应满足
4）令
，
，因此可以得到：
，
三、源程序：
#include"stdio.h"
#include"math.h"
#defineP0.1
voidmain（）
{
d
升级会员 