📄 +
字号:
//编程语言：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;
compressA[4][dimensionA-1]=0;//compressA还原结束
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 + -