數(shù)值分析大作業(yè)一

數(shù)值分析大作業(yè)一冪法反冪法求矩陣特征值2014.11尤萬(wàn)龍2014.11尤萬(wàn)龍PT1400133數(shù)值分析大作業(yè)一題目：設(shè)有501*501的實(shí)對(duì)稱(chēng)矩陣A其中，，b=0.16，c=-0.064。矩陣A的特征值，并且有 求的值求A的與數(shù)最接近的特征值。求A的（譜范數(shù)）條件數(shù)和行列式。算法設(shè)計(jì)由于C語(yǔ)言中數(shù)組角標(biāo)都是從0開(kāi)始的，所以在數(shù)組MatrixC[5][501]中檢索A的帶內(nèi)元素aij的方法是：A的帶內(nèi)元素aij=C中的元素ci-j+2,j。由于λ1?λ2?…?λ501，所以在以所有特征值建立的數(shù)軸上，λ1、λ501位于數(shù)軸的兩端，兩者之一必為按模最大。利用冪法，可以求出來(lái)按模最大的特征值，但該值可能為λ1也可能為λ501；上步求出按模最大的特征值λM后，將原矩陣平移λM,再利用冪法求一次平移后矩陣的按模最大的特征值，即是數(shù)軸上另一端點(diǎn)值λM′。比較λM與λM+λM′的大小，大的為矩陣A的最大特征值，小的為A矩陣的最小的特征值；利用反冪法，求矩陣A的按模最小的特征值。但是反冪法中要用到線性方程組的求解，而原矩陣A又是帶狀矩陣，采用LU分解。所以在這之前要定義一個(gè)LU分解子程序，將A矩陣分解為單位下三角矩陣L和上三角矩陣U的乘積。先利用循環(huán)求出k從1到39變化的uk的值。當(dāng)循環(huán)次數(shù)為39時(shí)，在每次循環(huán)時(shí)都將壓縮后的矩陣A的第三行減去相應(yīng)的，然后調(diào)用LU分解的子程序，利用反冪法求出與uk最接近的特征值，該特征值等于利用反冪法求出的值與uk的和。A的譜范數(shù)條件數(shù)等于按模最大的特征值的絕對(duì)值與按模最小的特征值的絕對(duì)值之比，按模最大的特征值與按模最小的特征值已分別在前面求出。A的行列式的值就是矩陣A進(jìn)行LU分解后U的對(duì)角線元素的乘積。先把帶狀線性矩陣A[501][501]轉(zhuǎn)存為一個(gè)矩陣c[5][501].源程序/*headerfilecomheader.h*/#ifndefCOMHEADER_H_INCLUDED#defineCOMHEADER_H_INCLUDED#include<iostream>#include<iomanip>#include<cmath>#include<fstream>//#include"matrixfuc.h"usingnamespacestd;constintcolumn_size=501;constdoubleErr=1.0e-12;#endif/*matrixfunctionclassdeclarematrixfuc.h*/#pragmaonce#include"comheader.h"classmatrixfuc{public: matrixfuc(void); ~matrixfuc(void); doublematrix_init(); doublepower_eig(); doubleinvpower_eig(doubleMatB[5][column_size]); voidmatrix_LU(doubleMatB[5][column_size]); double matA[5][column_size];};/*realizationofmatrixfucmatrixfuc.cpp*/#include"matrixfuc.h"intint_max2(inta,intb);intint_min2(inta,intb);intint_max3(inta,intb,intc);matrixfuc::matrixfuc(void){}matrixfuc::~matrixfuc(void){}/*初始化矩陣*/doublematrixfuc::matrix_init(void){ inti,j; for(i=0;i<5;i++) for(j=0;j<column_size;j++) { if(i==0||i==4) { matA[i][j]=-0.064; } elseif(i==1||i==3) { matA[i][j]=0.16; } elseif(i==2) { matA[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1)) -0.64*exp(0.1/(j+1)); } } matA[0][0]=matA[0][1]=matA[1][0]=0; matA[3][500]=matA[4][500]=matA[4][499]=0; return(matA[i][j]);}doublematrixfuc::power_eig(){ doubleu[column_size],unitu[column_size]; doublesum,lengthu,beta,temp; inti,j; beta=0.0; for(i=0;i<column_size;i++) u[i]=1; /*迭代初始向量*/ do { temp=beta; sum=0; for(i=0;i<column_size;i++) sum=sum+u[i]*u[i]; lengthu=sqrt(sum); for(i=0;i<column_size;i++) unitu[i]=u[i]/lengthu; for(i=1;i<=column_size;i++) { sum=0; for(j=int_max2(i-2,1);j<=int_min2(i+2,column_size);j++) sum=sum+matA[i-j+2][j-1]*unitu[j-1]; u[i-1]=sum; } sum=0; for(i=0;i<column_size;i++) sum=sum+unitu[i]*u[i]; beta=sum; }while(fabs((beta-temp)/beta)>=Err); return(beta);}doublematrixfuc::invpower_eig(doubleMatB[5][501]){ doubleu[column_size],unitu[column_size],unitu_u[column_size]; doublesum,lengthu,beta,temp; doubleMatD[5][column_size]; inti,j; beta=0.0; for(i=0;i<column_size;i++) u[i]=1; /*迭代初始向量*/ do { for(intx=0;x<5;x++) { for(inty=0;y<column_size;y++) { MatD[x][y]=MatB[x][y]; } } //matrix_init(); temp=beta; sum=0; for(i=0;i<column_size;i++) sum=sum+u[i]*u[i]; lengthu=sqrt(sum); for(i=0;i<column_size;i++) { unitu[i]=u[i]/lengthu; unitu_u[i]=unitu[i]; } matrix_LU(MatD); for(i=2;i<=column_size;i++) { sum=0; for(j=int_max2(i-2,1);j<=i-1;j++) sum+=MatD[i-j+2][j-1]*unitu_u[j-1]; unitu_u[i-1]-=sum; } u[column_size-1]=unitu_u[column_size-1]/MatD[2][column_size-1]; for(i=column_size-1;i>=1;i--) { sum=0; for(j=i+1;j<=int_min2(i+2,column_size);j++) sum+=MatD[i-j+2][j-1]*u[j-1]; u[i-1]=(unitu_u[i-1]-sum)/MatD[2][i-1]; } sum=0; for(i=0;i<column_size;i++) sum=sum+unitu[i]*u[i]; beta=sum; }while(fabs((temp-beta)/temp)>=Err); beta=1/beta; return(beta);}inlinevoidmatrixfuc::matrix_LU(doubleMatB[5][column_size]){ constints=2; constintr=2; doublesum; intk,i,j; for(k=1;k<=column_size;k++) { for(j=k;j<=int_min2(k+s,column_size);j++) { sum=0; for(i=int_max3(1,k-r,j-s);i<=k-1;i++) sum+=MatB[k-i+s][i-1]*MatB[i-j+s][j-1]; MatB[k-j+s][j-1]-=sum; } for(j=k+1;j<=int_min2(k+r,column_size);j++) { sum=0; for(i=int_max3(1,j-r,k-s);i<=k-1;i++) sum+=MatB[j-i+s][i-1]*MatB[i-k+s][k-1]; MatB[j-k+s][k-1]=(MatB[j-k+s][k-1]-sum)/MatB[s][k-1]; } }}intint_max2(inta,intb) { return(a>b?a:b);/*求兩個(gè)數(shù)字中最大值的子程序*/}inlineintint_min2(inta,intb){ return(a<b?a:b);/*求兩個(gè)數(shù)字中最小值的子程序*/}inlineintint_max3(inta,intb,intc){ intt; if(a>b) t=a;/*求三個(gè)數(shù)字中最大值的子程序*/ elset=b; if(t<c)t=c; return(t);}/***********************************************************Thisisaprogramforsolvingeigenvalueandeigenvectormatrix_eig.cpp***********************************************************/#include"comheader.h"#include"matrixfuc.h"/*******************************************************MainFunction********************************************************/intmain(){ doublea1,a2; // matrixfucmatrix_fuc; matrix_fuc.matrix_init(); a1=matrix_fuc.power_eig(); if(a1<0) { cout<<"矩陣A最小的特征值lambda1："<<endl; //Res<<"矩陣A最小的特征值lambda1："<<endl; } elseif(a1>=0) cout<<"矩陣A最大的特征值lambda501："<<endl; cout<<setiosflags(ios::scientific)<< setprecision(12)<<a1<<endl; for(inti=0;i<column_size;i++) matrix_fuc.matA[2][i]=matrix_fuc.matA[2][i]-a1; a2=matrix_fuc.power_eig()+a1; if(a2<0) cout<<"矩陣A最小的特征值lambda1："<<endl; elseif(a2>=0) cout<<"矩陣A最大的特征值lambda501："<<endl; cout<<setiosflags(ios::scientific)<<setprecision(12)<<a2<<endl; /*利用反冪法計(jì)算矩陣A的按模最小特征值*/ doublea3; matrix_fuc.matrix_init(); a3=matrix_fuc.invpower_eig(matrix_fuc.matA); cout<<"矩陣A按模最小的特征值lambdas:"<<endl; cout<<setiosflags(ios::scientific)<<setprecision(12)<<a3<<endl; matrix_fuc.matrix_init(); doubleu[501]={0}; doubleMatB[5][column_size]; cout<<"與數(shù)uk最接近的特征值："<<endl; for(inti=0;i<39;i++) { for(intz=0;z<column_size;z++) { u[i]=a1+(i+1)*(a2-a1)/40; } for(intz=0;z<column_size;z++) { matrix_fuc.matA[2][z]=matrix_fuc.matA[2][z]-u[z]; } for(intx=0;x<5;x++) { for(inty=0;y<column_size;y++) { MatB[x][y]=matrix_fuc.matA[x][y]; } } u[i]=matrix_fuc.invpower_eig(MatB)+u[i]; cout<<"lambda"<<"[i"<<i+1<<"]"<<" ";