?? +
字號:
//編程語言：JAVA
/**
* @(#)work.java
*
* work application《數值分析》計算實習題目第一題（運用冪法與反冪法）
*
* @author 503所 馬靈霞
* @version 1.00 2008/11/11
*/
import java.text.*;
import java.math.*;
import java.util.*;
public class work {
public static void main(String[] args) {
DecimalFormat myformat = new DecimalFormat();
myformat.applyPattern("0.000000000000E0000");
/*
*聲明全局變量matrixA的維數dimensionA,
*A的（譜范數）條件數condA，行列式detA,最大特征值eigenvalueMax,最小特征值eigenvalueMin
*模為最小的特征值absoluteEigenvalueMin,最大absoluteEigenvalueMax
*第二問benchmark[]=eigenvalueMin+k(eigenvalueMax-eigenvalueMin)/40,
*與benchmark[]最接近的特征值序列adjacenceBenchmark[]
*diagonalA[]存儲matrixA的對角線元素，constantB與constantC為matrixA的兩個元素，compressA[][]為matrixA的壓縮矩陣
**/
int dimensionA=501;
double condA=0,detA=0,eigenvalueMax=0,eigenvalueMin=0,absoluteEigenvalueMin=0,absoluteEigenvalueMax=0;
int dimensionBenchmark=40;
double benchmark[]=new double[dimensionBenchmark];
double adjacenceBenchmark[]=new double[dimensionBenchmark];
double diagonalA[]=new double[dimensionA];
double constantB=0.16;
double constantC=-0.064;
for(int i=0;i<=500;i++){
diagonalA[i]=(1.64-0.024*(i+1))*Math.sin(0.2*(i+1))-0.64*Math.exp(0.1/(i+1));
}//A的對角線元素的值
double compressA[][]=new double[5][dimensionA];
compressA[0][0]=0;
compressA[0][1]=0;
for(int i=2;i<=dimensionA-1;i++){
compressA[0][i]=constantC;
}
compressA[1][0]=0;
for(int i=1;i<=dimensionA-1;i++){
compressA[1][i]=constantB;
}
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=diagonalA[i];
}
for(int i=0;i<=dimensionA-2;i++){
compressA[3][i]=constantB;
}
compressA[3][dimensionA-1]=0;
for(int i=0;i<=dimensionA-3;i++){
compressA[4][i]=constantC;
}
compressA[4][dimensionA-2]=0;
compressA[4][dimensionA-1]=0;//將A壓縮存儲到compressA中
//對compressA進行LU分解，求A的行列式
for(int k=0;k<=dimensionA-1;k++){
for(int j=k;j<=Math.min(k+2,dimensionA-1);j++){
double sum1=0;
for(int t=Math.max(Math.max(0,k-2),j-2);t<=k-1;t++){
sum1=sum1+compressA[k-t+2][t]*compressA[t-j+2][j];
}
compressA[k-j+2][j]=compressA[k-j+2][j]-sum1;
}
for(int i=k+1;i<=Math.min(k+2,dimensionA-1);i++){
if(k>=dimensionA-1) break;
else{
double sum2=0;
for(int t=Math.max(Math.max(0,i-2),k-2);t<=k-1;t++){
sum2=sum2+compressA[i-t+2][t]*compressA[t-k+2][k];
}
compressA[i-k+2][k]=(compressA[i-k+2][k]-sum2)/compressA[2][k];
}
}
}// LU分解結束
detA=1;
for(int i=0;i<=dimensionA-1;i++){
detA=detA*compressA[2][i];
}//矩陣A的行列式求出
/*用冪法求矩陣matrixA的最大特征值eigenvalueMax與最小特征值eigenvalueMin
*聲明初始向量vectorU并賦初始值,
*局部變量absoluteU來存向量vectorU的模
*局部變量identityU來存vectorU單位化之后的向量
*former,current迭代產生的前后兩個值,用來判斷何時退出迭代
*局部變量e來存精度水平
*應還原compressA
**/
compressA[0][0]=0;
compressA[0][1]=0;
for(int i=2;i<=dimensionA-1;i++){
compressA[0][i]=constantC;
}
compressA[1][0]=0;
for(int i=1;i<=dimensionA-1;i++){
compressA[1][i]=constantB;
}
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=diagonalA[i];
}
for(int i=0;i<=dimensionA-2;i++){
compressA[3][i]=constantB;
}
compressA[3][dimensionA-1]=0;
for(int i=0;i<=dimensionA-3;i++){
compressA[4][i]=constantC;
}
compressA[4][dimensionA-2]=0;
compressA[4][dimensionA-1]=0;//compress還原結束
double vectorU[]=new double[dimensionA];
for(int i=0;i<=dimensionA-1;i++){
vectorU[i]=1;
}//聲明初始向量vectorU
double e=0.000000000001;// 精度水平e
double former=0,current=0,eigenvalue1=0,eigenvalue2=0;
//冪法迭代開始
do{
double absoluteU=0;
double identityU[]=new double[dimensionA];
double sum3=0;
for(int i=0;i<=dimensionA-1;i++){
sum3=sum3+vectorU[i]*vectorU[i];
}
absoluteU=Math.sqrt(sum3);
for(int i=0;i<=dimensionA-1;i++){
identityU[i]=vectorU[i]/absoluteU;
}
vectorU[0]=0;
for(int i=4;i>=2;i--){
vectorU[0]=vectorU[0]+compressA[i][0]*identityU[i-2];
}
vectorU[1]=0;
for(int i=4;i>=1;i--){
vectorU[1]=vectorU[1]+compressA[i][1]*identityU[i-1];
}
for(int k=2;k<=dimensionA-3;k++){
double sum4=0;
for(int i=4;i>=0;i--){
sum4=sum4+compressA[i][k]*identityU[i+k-2];
}
vectorU[k]=sum4;
}
vectorU[dimensionA-2]=0;
for(int i=3;i>=0;i--){
vectorU[dimensionA-2]=vectorU[dimensionA-2]+compressA[i][dimensionA-2]*identityU[i+dimensionA-4];
}
vectorU[dimensionA-1]=0;
for(int i=2;i>=0;i--){
vectorU[dimensionA-1]=vectorU[dimensionA-1]+compressA[i][dimensionA-1]*identityU[i+dimensionA-3];
}
double sum5=0;
for(int i=0;i<=dimensionA-1;i++){
sum5=sum5+identityU[i]*vectorU[i];
}
former=current;
current=sum5;
eigenvalue1=sum5;
absoluteEigenvalueMax=Math.abs(eigenvalue1);
}while((Math.abs(current-former)/Math.abs(current))>=e);//迭代結束
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=compressA[2][i]-eigenvalue1;
}//得出平移A后的壓縮矩陣
for(int i=0;i<=dimensionA-1;i++){
vectorU[i]=1;
}
former=current=0;
//對平移后的A迭代開始
do{
double absoluteU=0;
double identityU[]=new double[dimensionA];
double sum3=0;
for(int i=0;i<=dimensionA-1;i++){
sum3=sum3+vectorU[i]*vectorU[i];
}
absoluteU=Math.sqrt(sum3);
for(int i=0;i<=dimensionA-1;i++){
identityU[i]=vectorU[i]/absoluteU;
}
vectorU[0]=0;
for(int i=4;i>=2;i--){
vectorU[0]=vectorU[0]+compressA[i][0]*identityU[i-2];
}
vectorU[1]=0;
for(int i=4;i>=1;i--){
vectorU[1]=vectorU[1]+compressA[i][1]*identityU[i-1];
}
for(int k=2;k<=dimensionA-3;k++){
double sum4=0;
for(int i=4;i>=0;i--){
sum4=sum4+compressA[i][k]*identityU[i+k-2];
}
vectorU[k]=sum4;
}
vectorU[dimensionA-2]=0;
for(int i=3;i>=0;i--){
vectorU[dimensionA-2]=vectorU[dimensionA-2]+compressA[i][dimensionA-2]*identityU[i+dimensionA-4];
}
vectorU[dimensionA-1]=0;
for(int i=2;i>=0;i--){
vectorU[dimensionA-1]=vectorU[dimensionA-1]+compressA[i][dimensionA-1]*identityU[i+dimensionA-3];
}
double sum5=0;
for(int i=0;i<=dimensionA-1;i++){
sum5=sum5+identityU[i]*vectorU[i];
}
former=current;
current=sum5;
eigenvalue2=sum5;
}while((Math.abs(current-former)/Math.abs(current))>=e);//迭代結束
if(Math.signum(eigenvalue1)==1){
eigenvalueMax=eigenvalue1;
eigenvalueMin=eigenvalue1+eigenvalue2;
}
else{
eigenvalueMin=eigenvalue1;
eigenvalueMax=eigenvalue1+eigenvalue2;
}//判斷正負，得出最大/最小特征值
System.out.println("矩陣A的最大特征值為："+myformat.format(eigenvalueMax));
System.out.println("矩陣A的最小特征值為："+myformat.format(eigenvalueMin));
/*還原compressA,然后用反冪法求matrixA的特征值的最小模absoluteEigenvalueMin,以及第二問的與benchmark接近的特征值adjacenceBenchmark.
*其中absoluteEigenvalueMin,即為benchmark=0時的adjacenceBenchmark的絕對值
*聲明變量mid1[]存儲方程的中間解，processEigenvalue[]來存儲compressA-benchmark的模最小的特征值的值*/
double processEigenvalue[]=new double[40];
processEigenvalue[0]=0;
compressA[0][0]=0;
compressA[0][1]=0;
for(int i=2;i<=dimensionA-1;i++){
compressA[0][i]=constantC;
}
compressA[1][0]=0;
for(int i=1;i<=dimensionA-1;i++){
compressA[1][i]=constantB;
}
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=diagonalA[i];
}
for(int i=0;i<=dimensionA-2;i++){
compressA[3][i]=constantB;
}
compressA[3][dimensionA-1]=0;
for(int i=0;i<=dimensionA-3;i++){
compressA[4][i]=constantC;
}
compressA[4][dimensionA-2]=0;
benchmark[0]=0;
for(int i=1;i<=39;i++){
benchmark[i]=eigenvalueMin+i*(eigenvalueMax-eigenvalueMin)/40;
}//為benchmark賦值
//40次平移求特征值開始循環
for(int k=0;k<=39;k++){
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=compressA[2][i]-benchmark[k];
}//得出平移A后的壓縮矩陣
//對compressA進行LU分解
for(int n=0;n<=dimensionA-1;n++){
for(int j=n;j<=Math.min(n+2,dimensionA-1);j++){
double sum1=0;
for(int t=Math.max(Math.max(0,n-2),j-2);t<=n-1;t++){
sum1=sum1+compressA[n-t+2][t]*compressA[t-j+2][j];
}
compressA[n-j+2][j]=compressA[n-j+2][j]-sum1;
}
for(int i=n+1;i<=Math.min(n+2,dimensionA-1);i++){
if(n>=dimensionA-1) break;
else{
double sum2=0;
for(int t=Math.max(Math.max(0,i-2),n-2);t<=n-1;t++){
sum2=sum2+compressA[i-t+2][t]*compressA[t-n+2][n];
}
compressA[i-n+2][n]=(compressA[i-n+2][n]-sum2)/compressA[2][n];
}
}
}//LU分解結束
for(int i=0;i<=dimensionA-1;i++){
vectorU[i]=1;
}//初始向量賦值
former=current=0;
//反冪法求平移后矩陣模最小的特征值開始迭代
do{
double absoluteU=0;
double identityU[]=new double[dimensionA];
double sum6=0;
for(int i=0;i<=dimensionA-1;i++){
sum6=sum6+vectorU[i]*vectorU[i];
}
absoluteU=Math.sqrt(sum6);
for(int i=0;i<=dimensionA-1;i++){
identityU[i]=vectorU[i]/absoluteU;
}
//求下一次迭代需要的vectorU
double mid1[]=new double[dimensionA];
mid1[0]=identityU[0];
mid1[1]=identityU[1]-compressA[3][0]*mid1[0];
for(int i=2;i<=dimensionA-1;i++){
mid1[i]=identityU[i]-compressA[4][i-2]*mid1[i-2]-compressA[3][i-1]*mid1[i-1];
}
vectorU[dimensionA-1]=mid1[dimensionA-1]/compressA[2][dimensionA-1];
vectorU[dimensionA-2]=mid1[dimensionA-2]-compressA[1][dimensionA-1]*vectorU[dimensionA-1];
for(int i=dimensionA-3;i>=0;i--){
vectorU[i]=(mid1[i]-compressA[1][i+1]*vectorU[i+1]-compressA[0][i+2]*vectorU[i+2])/compressA[2][i];
}//vectorU求出
double sum7=0;
for(int i=0;i<=dimensionA-1;i++){
sum7=sum7+identityU[i]*vectorU[i];
}
former=current;
current=1/sum7;
}while((Math.abs(current-former)/Math.abs(current))>=e);//迭代結束
processEigenvalue[k]=current;
adjacenceBenchmark[k]=processEigenvalue[k]+benchmark[k];
if(k!=0){
System.out.println("矩陣A與第"+k+"個數"+myformat.format(benchmark[k])+"最接近的特征值為："+myformat.format(adjacenceBenchmark[k]));
}
//還原compressA
compressA[0][0]=0;
compressA[0][1]=0;
for(int i=2;i<=dimensionA-1;i++){
compressA[0][i]=constantC;
}
compressA[1][0]=0;
for(int i=1;i<=dimensionA-1;i++){
compressA[1][i]=constantB;
}
for(int i=0;i<=dimensionA-1;i++){
compressA[2][i]=diagonalA[i];
}
for(int i=0;i<=dimensionA-2;i++){
compressA[3][i]=constantB;
}
compressA[3][dimensionA-1]=0;
for(int i=0;i<=dimensionA-3;i++){
compressA[4][i]=constantC;
}
compressA[4][dimensionA-2]=0;
compressA[4][dimensionA-1]=0;//compress還原結束
}//40個特征值求出
absoluteEigenvalueMin=adjacenceBenchmark[0];
condA=Math.abs(absoluteEigenvalueMax/absoluteEigenvalueMin);
System.out.println("矩陣A的特征值的模的最小值為："+myformat.format(absoluteEigenvalueMin));
System.out.println("矩陣A的（譜范數）條件數為："+myformat.format(condA));
System.out.println("矩陣的行列式為："+myformat.format(detA));
}
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -