#include "stdafx.h"

#include <math.h>

// 矩陣特徵值

int Maximum(double a[], int n)

{

int p=0;

int i;

for (i=1;i<n;i++) {

if (abs(a[i])>abs(a[p]))

p = i;

}

return p;

}

void MiChengMethod()

{

printf("求矩陣特徵值: 冪乘法\n");

double A[3][3] = {{4,-1,0},{0,4,-1},{0,-1,4}};

double x0[3] = {0,0,1};

double x1[3];

double lamda0,lamda1;

int k,i,j;

int n=3;

double epsilon = 1e-3;

int p;

p = Maximum(x0,n);

lamda0 = x0[p];

for (i=0;i<n;i++)

x0[i] /= lamda0;       

k = 0;

for (;;) {

for (i=0;i<n;i++) {

x1[i] = 0;

for (j=0;j<n;j++)

x1[i] += A[i][j]*x0[j];        

}

p = Maximum(x1,n);

lamda1 = x1[p];

k++;       

if (k>1&&abs(lamda1-lamda0)<epsilon)

{

for (i=0;i<n;i++) x1[i]=x1[i]/lamda1;  // 特徵向量

break;

}

else

{

lamda0 = lamda1;   // 迭代

for (i=0;i<n;i++) x0[i]=x1[i]/lamda1;

}

}

PS(&k,1,"迭代次數");

PS(&lamda1,1,"按模特徵值lamda");

PS(x1,n,"特徵向量");

}

void

 PingYiYuanDianMethod()

{

printf("冪法的加速: 平移原點法\n");

double A[3][3] = {{4,-1,0},{0,4,-1},{0,-1,4}};

double x0[3] = {0,0,1};

double x1[3];

double lamda0,lamda1;

int k,i,j;

int n=3;

double epsilon = 1e-3;

int p;

double mu=2.9; //  平移量

for (i=0;i<n;i++)    // 平移

A[i][i]-=mu;

p = Maximum(x0,n);

lamda0 = x0[p];

for (i=0;i<n;i++)

x0[i] /= lamda0;

k = 0;

for (;;) {

for (i=0;i<n;i++) {

x1[i] = 0;

for (j=0;j<n;j++)

x1[i] += A[i][j]*x0[j];        

}

p = Maximum(x1,n);

lamda1 = x1[p];

k++;       

if (k>1&&abs(lamda1-lamda0)<epsilon)

{

for (i=0;i<n;i++) x1[i]=x1[i]/lamda1;  // 特徵向量

lamda1 += mu;

break;

}

else

{

lamda0 = lamda1;   // 迭代

for (i=0;i<n;i++) x0[i]=x1[i]/lamda1;

}

}

PS(&k,1,"迭代次數");

PS(&lamda1,1,"按模特徵值lamda");

PS(x1,n,"特徵向量");

}

void AitkenJiaSuMethod()

{

printf("冪法的加速: 埃特肯加速法\n");

double A[3][3] = {{4,-1,0},{0,4,-1},{0,-1,4}};

double x0[3] = {0,0,1};

double x1[3];

double lamda[3];

int k,i,j;

int n=3;

double epsilon = 1e-3;

int p;

double lam[2]={0}; // 加速

double maxv;

p = Maximum(x0,n);

lamda[0] = x0[p];

for (i=0;i<n;i++)

x0[i] /= lamda[0];     

k = 0;

for (;;)

{

for (i=0;i<n;i++) {

x1[i] = 0;

for (j=0;j<n;j++)

x1[i] += A[i][j]*x0[j];        

}

p = Maximum(x1,n);

maxv = x1[p];

k++;       

if (k<3)

lamda[k] = maxv;

else

{

lamda[0]=lamda[1];

lamda[1]=lamda[2];

lamda[2]=maxv;

}

if (k>2) // 加速公式

{

lam[0] = lam[1];

lam[1] = lamda[0]

-(lamda[1]-lamda[0])*(lamda[1]-lamda[0])/(lamda[2]-2*lamda[1]+lamda[0]);

}

if (k>3&&abs(lam[1]-lam[0])<epsilon)

break;

else

for (i=0;i<n;i++) x0[i]=x1[i]/maxv;

}

PS(&k,1,"迭代次數");

printf("按模特徵值lamda=%.7f\n",lam[1]);

}

void FanMiMethod()

{

printf("反冪法: 求矩陣最小特徵值\n");

double a[3][3] = {{1,-1,2},{-2,0,5},{6,-3,6}};

double b[3] = {1,1,1};

double l[3][3],u[3][3];

double y[3]={0};

int n=3;

int i,j;   

double x[3];

double lamda0,lamda1;

double epsilon=1e-3;

for (i = 0; i < n; i++) {

for (j = 0; j < n; j++) {

l[i][j] = u[i][j] = 0.;

}

}

// LU分解

int k = 0;;

int r;

for (k = 0; k < n; k++) {

// 先計算u

for (j = k; j < n; j++) {

u[k][j] = a[k][j];

for (r = 0; r < k && k > 0; r++) {

u[k][j] -= l[k][r]*u[r][j];

}

}

// 後計算l

l[k][k] = 1.;

for (i = k+1; i < n; i++) {

l[i][k] = a[i][k];

for (r = 0; r < k && k > 0; r++) {

l[i][k] -= l[i][r]*u[r][k];

}

l[i][k] /= u[k][k];

}

}

printf("l=\n");

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

printf("%f\t",l[i][j]);

printf("\n");

}

printf("u=\n");

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

printf("%f\t",u[i][j]);

printf("\n");

}

int p= Maximum(b,n);

lamda0=b[p];

for (i=0;i<n;i++) {

b[i]/=lamda0;

}

// 迭代

k=0;

for(;;)

{

// 回代求解

// 回代 求y

// Ly = b

y[0] = b[0];

for (i = 1; i < n; i++) {

y[i] = b[i];

for (j = 0; j < i; j++) {

y[i] = y[i] - l[i][j]*y[j];

}

}

// 求x

//Ux=y

x[n-1] = y[n-1]/u[n-1][n-1];

for (i = n-2; i >=0; i--) {

x[i] = y[i];

for (j = i+1; j < n; j++) {

x[i] = x[i] - u[i][j]*x[j];

}

x[i] /= u[i][i];

}

k++;

p=Maximum(x,n);

lamda1 =x[p];

printf("t%d=%f\n",k,lamda1);

printf("y%d= ",k);

for (i = 0; i < n; i++) {

printf("%.4f\t",x[i]/lamda1);

}

printf("\n");

if (abs(lamda1-lamda0)<epsilon)

{

break;

}

else

{

lamda0=lamda1;

for (i=0;i<n;i++) {

b[i]=x[i]/lamda1;

}

}

}

}

void FanMiMethod_2()

{

printf("反冪法: 已知某特徵值的近似值求更精確的特徵值\n");

double a[3][3] = {{4,-1,0},{0,4,-1},{0,-1,4}};

double b[3] = {0,0,1};

double l[3][3],u[3][3];

double y[3]={0};

int n=3;

int i,j;   

double x[3];

double lamda0,lamda1;

double epsilon=1e-3;

double lam=4.93;

for (i = 0; i < n; i++) {

a[i][i] -= lam;              // 減掉對角元

}

// LU分解, 先算u,後算l

int k = 0;;

int r;

for (k = 0; k < n; k++) {

// 先計算u, 行k

for (j = 0; j < n; j++) {

if (j<k)

u[k][j]=0.;

else

{

u[k][j] = a[k][j];

for (r = 0; r < k && k > 0; r++) {

u[k][j] -= l[k][r]*u[r][j];

}

}

}

// 後計算l,列k

for (i = 0; i < n; i++)

{

if (i<k) l[i][k]=0.;

else if(i==k) l[i][k]=1.;

else

{

l[i][k] = a[i][k];

for (r = 0; r < k && k > 0; r++) {

l[i][k] -= l[i][r]*u[r][k];

}

l[i][k] /= u[k][k];

}

}

}

printf("l=\n");

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

printf("%f\t",l[i][j]);

printf("\n");

}

printf("u=\n");

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

printf("%f\t",u[i][j]);

printf("\n");

}

int p= Maximum(b,n);

lamda0=b[p];

for (i=0;i<n;i++) {

b[i]/=lamda0;

}

// 迭代

k=0;

for(;;)

{

// 回代求解

// 回代 求y

// Ly = b

y[0] = b[0];

for (i = 1; i < n; i++) {

y[i] = b[i];

for (j = 0; j < i; j++) {

y[i] = y[i] - l[i][j]*y[j];

}

}

// 求x

//Ux=y

x[n-1] = y[n-1]/u[n-1][n-1];

for (i = n-2; i >=0; i--) {

x[i] = y[i];

for (j = i+1; j < n; j++) {

x[i] = x[i] - u[i][j]*x[j];

}

x[i] /= u[i][i];

}

k++;

p=Maximum(x,n);

lamda1 =x[p];

printf("t%d=%f\n",k,lamda1);

//      printf("t%d=%f\n",k,1./lamda1+lam);

printf("y%d= ",k);

for (i = 0; i < n; i++) {

printf("%f\t",x[i]/lamda1);

}

printf("\n");

if (abs(lamda1-lamda0)<epsilon)

break;

else

{

lamda0=lamda1;

for (i=0;i<n;i++) {

b[i]=x[i]/lamda1;

}

}

}

printf("更精確的特徵值=%f\n",1./lamda1+lam);

}

// 將矩陣化成二次型的方法

void JacobiMethod()

{

printf("求實對稱矩陣特徵值的雅可比方法:\n");

double a[3][3] = {{-1,2,1},{2,-4,1},{1,1,-6}};

double b[3][3];

int i,j;

int p,q;

int n=3;

double epsilon = 5e-10;

double w;

double mu,lamda,s,c;

int k=0;

double B[3][3]={{1,0,0},{0,1,0},{0,0,1}};  // 特徵向量

double U[3][3]={{1,0,0},{0,1,0},{0,0,1}};  // 旋轉矩陣

double Tmp[3][3];

printf("輸入矩陣a%d=\n",k);

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

printf("%f\t",a[i][j]);

printf("\n");

}

for (;;) {

// 找到矩陣中非對角元最大值位置(p,q)

p = 0;

q = 1;

for (i=0;i<n;i++) {

for (j=0;j<n;j++) {

if (j!=i&&abs(a[i][j])>abs(a[p][q]))

{

p=i;q=j;

}

}

}

// 如果最大元充分小,結束

if (abs(a[p][q])<epsilon)

break;

lamda = 2*a[p][q];

mu = (a[p][p]-a[q][q]);

//      printf("tan=%f\n",lamda/mu);

if (abs(mu)<epsilon)

{

w=1.;

}

else

{

if (mu>0)

w=lamda/sqrt(lamda*lamda+mu*mu);

else

w=-lamda/sqrt(lamda*lamda+mu*mu);

}

s=w/sqrt(2*(1+sqrt(1-w*w)));

c=sqrt(1-s*s);

// 做一次旋轉

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

b[i][j]=a[i][j];

}

// p行q行

for (j=0;j<n;j++)

{

if (j!=p&&j!=q)

{

b[p][j] = a[p][j]*c+a[q][j]*s;

}

if (j!=p&&j!=q)

{

b[q][j]= -a[p][j]*s+a[q][j]*c;

}

}

// p列q列

for (i=0;i<n;i++)

{

if (i!=p&&i!=q)

{

b[i][p] = a[i][p]*c+a[i][q]*s;

}

if (i!=p&&i!=q)

{

b[i][q]= -a[i][p]*s+a[i][q]*c;

}

}

// 對角元

b[p][p]=a[p][p]*c*c+a[q][q]*s*s+2*a[p][q]*s*c;

b[q][q]=a[p][p]*s*s+a[q][q]*c*c-2*a[p][q]*s*c;

b[p][q]=b[q][p] = 0.;

//      b[p][q]=b[q][p] = (a[q][q]-a[p][p])*s*c+a[p][q]*(2*c*c-1);

// 特徵向量

for (i=0;i<n;i++) {

U[i][i]=1;

for (j=0;j<n;j++) {

if (j!=i)

{

U[i][j]=0;

}

}

}

U[p][p] = U[q][q]=c;

if (p<q)

{

U[p][q] = -s;

U[q][p] = s;

}

else

{

U[p][q] = s;

U[q][p] = -s;

}

for (i=0;i<n;i++) {

for (j=0;j<n;j++) {

Tmp[i][j] =0;

for (int k=0;k<n;k++) {

Tmp[i][j]+=B[i][k]*U[k][j];

}

}

}

k++;

// 如果旋轉後的值非常小, 置為0

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

{

if (i!=j&&abs(b[i][j])<5e-10)

{

b[i][j] = 0;

}

a[i][j]=b[i][j];

B[i][j]=Tmp[i][j];

}

}

}

printf("a%d=\n",k);

for (i = 0; i < n; i++)

{

printf("%.14f\t",a[i][i]);

}

printf("\n");

printf("B%d=\n",k);

for (i = 0; i < n; i++)

{

for (j = 0; j < n; j++)

{

printf("%.14f\t",B[i][j]);

}

printf("\n");

}

}