北航数值分析大作业一Word下载.doc
- 文档编号:3643854
- 上传时间:2023-05-02
- 格式:DOC
- 页数:10
- 大小:306.23KB
北航数值分析大作业一Word下载.doc
《北航数值分析大作业一Word下载.doc》由会员分享,可在线阅读,更多相关《北航数值分析大作业一Word下载.doc(10页珍藏版)》请在冰点文库上搜索。
使用位移的方式求得另一特征值即可。
2求A的与数最接近的特征值。
题目可看成求以为偏移量后,按模最小的特征值。
即以为偏移量做位移,使用反幂法求出按模最小特征值后,加上,即为所求。
3求A的(谱范数)条件数和行列式detA。
矩阵A为非奇异对称矩阵,可知,
(1-1)
其中为按模最大特征值,为按模最小特征值。
detA可由LU分解得到。
因LU均为三角阵,则其主对角线乘积即为A的行列式。
二算法实现
1幂法
使用如下迭代格式:
(2-1)
终止迭代的控制理论使用,
实际使用
(2-2)
由于不保存A矩阵中的零元素,只保存主对角元素a[501]及b,c值。
则上式中简化为:
(2-3)
2反幂法
(2-4)
其中,解方程求出。
求解过程中使用LU分解,由于A为5对角矩阵,选择追赶法求取LU分解。
求解过程如下:
追赶法求LU分解的实现:
(2-5)
由上式推出分解公式如下:
(2-6)
推导出回代求解公式如下:
(2-7)
(2-8)
3及A行列式求解
(2-9)
由式(2-5)可得:
(2-10)
三源程序
#include<
stdio.h>
math.h>
doubleep=1e-12,b=0.16,c=-0.064;
intj=0;
doublepower(doublea[501]);
//幂法
doubleinv_power(doublea[501]);
//反幂法
doubledet(doublea[501]);
//求det
intmain() //主程序
{
inti,k;
doubleA[501],B[501],beta_1,beta_501,beta_s,beta_k;
doublemu;
for(i=0;
i<
501;
i++)
A[i]=(1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));
beta_1=power(A);
//第一问
printf("
λ1\t=%.12e\t迭代次数:
%d\n"
beta_1,j);
i++)//位移
B[i]=A[i]-beta_1;
beta_501=power(B)+beta_1;
printf("
λ501\t=%.12e\t迭代次数:
beta_501,j);
beta_s=inv_power(A);
λs\t=%.12e\t迭代次数:
beta_s,j);
for(k=1;
k<
=39;
k++)//第二问
{
mu=beta_1+k*(beta_501-beta_1)/40;
for(i=0;
B[i]=A[i]-mu;
beta_k=inv_power(B)+mu;
printf("
λi%d\t=%.12e\t迭代次数:
k,beta_k,j);
}
cond(A)2=%.12e\n"
beta_1/beta_s);
//第三问
detA\t=%.12e\n"
det(A));
}
doublepower(doublea[501])//幂法
inti=0,N=5000;
doubleb=0.16,c=-0.064;
doubleu[501],y[501];
doublem=1,beta;
u[i]=1;
j=0;
while(j<
N)
{
for(i=0;
{
y[i]=u[i]/fabs(m);
}
u[0]=a[0]*y[0]+b*y[1]+c*y[2];
u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];
u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];
u[500]=c*y[498]+b*y[499]+a[500]*y[500];
for(i=2;
499;
{u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];
}
beta=0;
for(i=0;
{
if(fabs(u[i])>
=fabs(beta))
beta=u[i];
}
if(beta<
0)
if(fabs(fabs(beta)-fabs(m))/fabs(beta)<
ep)
break;
if(fabs(beta-m)/fabs(beta)<
break;
m=beta;
j++;
}
returnbeta;
doubleinv_power(doublea[501])//反幂法
doublep[501],r[501],t[501],q[501],u[501],y[501];
doublebeta,m=1;
inti,N=1000;
p[0]=a[0];
t[0]=b/p[0];
r[1]=b;
p[1]=a[1]-r[1]*t[0];
q[0]=c/p[0];
q[1]=c/p[1];
t[1]=(b-r[1]*q[0])/p[1];
for(i=2;
r[i]=b-c*t[i-2];
p[i]=a[i]-c*q[i-2]-r[i]*t[i-1];
q[i]=c/p[i];
t[i]=(b-r[i]*q[i-1])/p[i];
}
u[i]=1;
j=0;
y[i]=u[i]/fabs(m);
u[0]=y[0]/p[0];
u[1]=(y[1]-r[1]*u[0])/p[1];
for(i=2;
u[i]=(y[i]-c*u[i-2]-r[i]*u[i-1])/p[i];
u[499]=u[499]-t[499]*u[500];
for(i=498;
i>
=0;
i--)
u[i]=u[i]-t[i]*u[i+1]-q[i]*u[i+2];
beta=0;
if(beta<
return1/beta;
doubledet(doublea[501])//求det
doubledet_A=1;
doublep[501],r[501],t[501],q[501];
inti;
}
det_A=det_A*p[i];
returndet_A;
四程序结果
五计算过程中的现象
使用作为终止迭代条件时,出现迭代无法终止的情况,通过调试发现按模最大特征值为负时,当k充分大后,迭代向量各分量不断变号,使得与异号,判别式不收敛。
因此将终止迭代条件修改为,程序实现如下:
if(beta<
if(fabs(beta-m)/fabs(beta)<
break;
从迭代次数可以看出与收敛较慢,由按模最大特征值与按模次大特征值的比值越小,收敛速度越慢,可知存在与和的模相近的特征值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业
![提示](https://static.bingdoc.com/images/bang_tan.gif)