北航数值分析大作业一Word文档格式.docx
- 文档编号:22760983
- 上传时间:2023-02-05
- 格式:DOCX
- 页数:15
- 大小:136.30KB
北航数值分析大作业一Word文档格式.docx
《北航数值分析大作业一Word文档格式.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业一Word文档格式.docx(15页珍藏版)》请在冰豆网上搜索。
分别为按模最大的特征值和按模最小的特征值
⑦本题中LU分解采用doolittle方法,分解后所得
矩阵即为对A进行初等行列式变换后所得矩阵,切变换过程中行列式值不变,因此
2、源程序
#include<
math.h>
stdio.h>
stdlib.h>
#defineN501//矩阵为501*501的矩阵
#defines2//上半带宽
#definer2//下半带宽3个宏定义为方便LU分解及求解方程组过程
/**********全局变量定义**********************/
doubleA[5][501];
doubleu[501],y[501];
doublelambda[41];
//lambda[0]为λ1及最小特征值,lambda[40]为λ501及最大特征值
doublelambda_s,lambda_m;
//按模最小(大)的特征值;
doubleLU[5][501];
/********子函数声明**************************/
voidinit_A(doubleA[][501]);
//初始化A矩阵
doublemodule_value_u(doublett[]);
//求u[501]的模值
voidinit_u(doublett[]);
//初始化迭代初始向量u0
doublepower_method(doubleoffset);
//带原点偏移的幂法,返回值为特征值
doubleinverse_power_method(doubleoffset);
//带原点偏移的反幂法,返回值为特征值,子函数中打印出偏移量,求得的特征值,误差,迭代次数
voidLU_decomposition(doublec[][501]);
//参数c为矩阵LU[5][501]首地址,程序进行完后,保存分解后的L和U缩减后的矩阵
voidsolve(doublec[][501],doubleb[],doublex[]);
////第1个参数为LU分解完的LU矩阵,第2个参数b为已知的右端值,第3个参数x为求得的解向量存储位置
intmax_2(inta,intb);
intmax_3(inta,intb,intc);
intmin_2(inta,intb);
voidmain()//主程序
{
inti;
doubleuk;
//偏移量
doubledet;
//A的行列式的值
init_A(A);
//初始化矩阵A
lambda_m=power_method(0);
//求按模最大的特征值
lambda[40]=power_method(lambda_m);
//求相反方向另一个端点的特征值
if(lambda_m>
lambda[40])//若大小反向,则交换两个元素中的值,得到λ1和λ501
{
lambda[0]=lambda[40];
lambda[40]=lambda_m;
}
else
lambda[0]=lambda_m;
lambda_s=inverse_power_method(0);
det=1;
for(i=0;
i<
N;
i++)
det=det*LU[s][i];
for(i=1;
40;
uk=lambda[0]+(lambda[40]-lambda[0])*i/40;
lambda[i]=inverse_power_method(uk);
printf("
-----------Theresultsareasfollows-------------\n"
);
λ[1]=%1.11e\nλ[501]=%1.11e\n"
lambda[0],lambda[40]);
λs=%1.11e\n"
lambda_s);
cond(A)=%1.11e\n"
fabs(lambda_m/lambda_s));
detA=%1.11e\n"
det);
for(i=1;
printf("
λ[i%d]=%1.11e"
i,lambda[i]);
if(i%3==0)
printf("
\n"
printf("
}
voidinit_A(doubleA[][501])//初始化A矩阵
inti,j;
5;
if(abs(i-2)==0)
{
for(j=0;
j<
501;
j++)
A[i][j]=(1.64-0.024*(j+1))*sin(0.2*j+0.2)-0.64*exp(0.1/(j+1));
}
elseif(abs(i-2)==1)
A[i][j]=0.16;
else
A[i][j]=-0.064;
voidinit_u(doublett[])//初始化迭代初始向量,自变量为初始向量的数组
tt[i]=1;
doublemodule_value_u(doublett[])//求u[501]的模值
inti=501;
doublet=0;
t=t+tt[i]*tt[i];
returnsqrt(t);
doublepower_method(doubleoffset)//带原点偏移的幂法,返回值为特征值
doubleb=0,b0=0,e;
doublem;
//求得的向量模值
inti,j,k=0;
//k表示迭代次数
init_u(u);
do
k++;
b0=b;
m=module_value_u(u);
for(i=0;
y[i]=u[i]/m;
u[i]=0;
for(j=max_2(0,i-s);
min_2(i+s+1,N);
{
u[i]=A[i-j+2][j]*y[j]+u[i];
if(i==j)
u[i]=u[i]-offset*y[j];
}
b=0;
b=b+y[i]*u[i];
}while(fabs(e=b-b0)>
1e-12);
b=b+offset;
λ=%fe=%ek=%d\n"
b,e,k);
//分别为特征值,迭代误差,迭代次数
return(b);
voidLU_decomposition(doublec[][501])//参数c为矩阵LU[5][501]首地址,程序进行完后,保存分解后的L和U缩减后的矩阵
intj,k,t;
for(k=0;
k<
k++)//求U矩阵的第k行以及L矩阵的第k列
for(j=k;
min_2(k+s,N-1)+1;
j++)//
for(t=max_3(0,k-r,j-s);
t<
k;
t++)//求U矩阵的行中各个元素的循环
c[k-j+s][j]=c[k-j+s][j]-c[k-t+s][t]*c[t-j+s][j];
}//每次计算完U的元素后才能计算L的元素,所以上下两个循环不能合并在一起
if(c[s][k]==0)
printf("
a[%d][%d]=0,Unabletoshowthesolutionofequations\n"
k,k);
exit
(1);
if(j<
min_2(k+s,N-1)+1)//因为矩阵L没有N行,因此这里要加一个判断,符合条件才能进行下面的循环语句
for(t=max_3(0,j-r+1,k-s);
t++)//求L矩阵的行中各个元素的循环
{
c[j-k+s+1][k]=c[j-k+s+1][k]-c[j-t+s+1][t]*c[t-k+s][k];
}
c[j-k+s+1][k]=c[j-k+s+1][k]/c[s][k];
voidsolve(doublec[][501],doubleb[],doublex[])//第1个参数为分解完的LU矩阵,第2个参数b为已知的右端值,第3个参数x为求得的解向量存储位置
{inti,t;
i++)//此循环求出向量y
x[i]=b[i];
for(t=max_2(0,i-r);
i;
t++)
x[i]=x[i]-c[i-t+s][t]*x[t];
for(i=N-1;
i>
=0;
i--)
for(t=i+1;
min_2(i+s,N-1)+1;
x[i]=x[i]/c[s][i];
doubleinverse_power_method(doubleoffset)//带原点偏移的反幂法,返回值为特征值
if(i!
=2)
LU[i][j]=A[i][j];
LU[i][j]=A[i][j]-offset;
LU_decomposition(LU);
//LU分解,只分解一次即可
y[i]=u[i]/module_value_u(u);
solve(LU,y,u);
b=1/b;
b=b+offset;
offset=%1.11eλ=%1.11ee=%1.11ek=%d\n"
offset,b,e,k);
returnb;
intmax_2(inta,intb)
if(a<
b)
a=b;
returna;
intmax_3(inta,intb,intc)
c)
a=c;
intmin_2(inta,intb)
if(a>
3、计算结果如下所示
在幂法和反幂法中偏移量、特征值、误差及计算迭代次数也在计算过程中打印出来
4、迭代过程中的现象、原因及解决方法
迭代的初始向量对计算结果可能产生极为重要的影响,如果选取初始向量不恰当,则可能得到错误的结果。
以最大(或最小)特征值——也即按模最大特征值为例:
分别取一下两组初始向量:
1.
2.
当初始向量为
时可以得到按模最大的特征向量为:
由以上可以看出,初始向量对于计算结果有着十分正要的影响。
原因分析:
初始向量可以表示成:
;
假设
,则采用幂法迭代k次后有:
,只有当
,时,才有
即幂法能够正常使用的前提是
,当初始向量中有若干0时,极有可能出现
的情况,此时有两种可能:
一种可能是由于舍入误差影响,迭代若干次后x1方向上分量不为零,相当于重新开始迭代并最终得到正确的结果,但这种情况下迭代次数较多,而另一种情况则是直接得出错误的结果。
为了保证
成立,最好取初始向量的每个分量都不为零。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业