数值分析大作业c语言版,数值分析大作业

以下为《数值分析大作业》的无排版文字预览，完整格式请下载

下载前请仔细阅读文字预览以及下方图片预览。图片预览是什么样的，下载的文档就是什么样的。

数值分析大作业 二

姓名：苗某某 单位：航天五院502所

算法设计方案

本次大作业主要涉及到的算法有：

矩阵A的拟上三角化；

矩阵A的QR分解；

带双步位移QR分解求解矩阵A的全部特征值；

列主元高斯消去法求解线性方程组。

其中，列主元高斯消去法用于求解实特征值所对应的特征向量。算法的具体应用和算法描述如下。

求矩阵A的拟上三角化矩阵A(n-1)

矩阵A的拟上三角化算法步骤如下：

记，并记的第r列到第n列元素为。

对于执行如下算法：

若全为零，则令，转(5)；否则转(2)。

计算

令

计算

继续

当此算法执行完后，就得到矩阵A的拟上三角化矩阵A(n-1)

A(n-1)的QR分解以求Q、R、RQ

记，并记。

对于执行

(1) 若全为零，则令，转(5)；否则转(2)；

(2) 计算

(3) 令

(4)计算

(5)继续

经过上述运算之后，就得到和上三角矩阵

然后利用和相乘可计算。

求矩阵A的全部特征值Lamta(i)

(1) 调用矩阵A拟上三角化算法得到的拟上三角矩阵，给定精度水平和迭代的最大次数100；

(2) 记，令k=1，m=10;

(3) 判断是否，是，则为特征值，m=m-1，goto(4)，否则goto(5);

(4) 如果m=1，则为特征值，goto(10)；若m=0，goto(10)；如果m>1，goto(3);

(5) 如果m=2，则得到A的两个特征值和, 即右下角二阶子阵对应的特征方程的两个根，goto(11);否则goto(6);

(6) 如果，则得到A的两个特征值两个特征值和，即右下角二阶子阵对应的特征方程的两个根，令m=m-2，goto(4);否则goto(7);

(7)如果k>100，则求解失败，程序退出;否则goto(8);

(8)记计算

调用QR分解算法对进行QR分解得到矩阵Q；

Ak+1 = QkTAkQk

(9)置k=k+1,goto(3);

(10)A的所有特征值已经计算完毕，停止计算。

此算法执行完后即求得矩阵A的所有特征值。

求矩阵A相应实特征值的特征向量

利用列主元高斯消去法求解Lamta(i)所对应的特征向量，需要注意的地方是，由于特征向量为非0向量，所以求解回带过程时取x10 = 1。

1. 消元过程

(1)选行号，使 。

(2)交换与。

(3)对于计算

2. 回代过程

x10 = 1

经过上述过程解线性方程组所得的解向量即为对应于特征值Lamta(i)的特征向量。

源程序：

#include

#include

#include

#define DIMENSION 10

#define Failure 0

#define Success 1

#define Standard 1

#define Ephasilong pow(10,12)

#define CompareTimes 100

void intialMatrixA(double A[DIMENSION][DIMENSION]);

void intialMatrix_Lamta(double Lamta[DIMENSION][2]);

void printMatrix(double matrix[DIMENSION][DIMENSION]);

void printEigenValue(double Lamta[DIMENSION][2]);

void intialMatrixQ_toUnitE(double Q[DIMENSION][DIMENSION]);

void intialMatrixR_toZero(double R[DIMENSION][DIMENSION]);

//拟上三角化子程序

int allIsZero(double A[DIMENSION][DIMENSION], int rowStartIndex,

int rowEndIndex, int clomn);

double caculate_dr(double A[DIMENSION][DIMENSION],

int rowStart, int rowEnd, int clomn);

int sgn(double a);

void getArray_U(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double c, int index);

void getArray_P_Q(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double h, double P[DIMENSION], double Q[DIMENSION]);

double get_T(double P[DIMENSION], double U[DIMENSION], double h);

void getArray_W(double Q[DIMENSION], double t,

double U[DIMENSION], double W[DIMENSION]);

void getMatrixA(double A[DIMENSION][DIMENSION], double W[DIMENSION],

double U[DIMENSION], double P[DIMENSION]);

void simulateUpTriangel(double A[DIMENSION][DIMENSION]);

//QR分解法所用子程序

void getArray_QR_U(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double c, int index, int M_Dimension);

void getArray_Omega(double Q[DIMENSION][DIMENSION], double U[DIMENSION],

double Omega[DIMENSION], int M_Dimension);

void getArray_Q(double Q[DIMENSION][DIMENSION], double Omega[DIMENSION],

double U[DIMENSION], double h, int M_Dimension);

void getArray_P(double A[DIMENSION][DIMENSION], double U[DIMENSION], double h,

double P[DIMENSION], int M_Dimension);

void getMatrixAr(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double P[DIMENSION], int M_Dimension);

void matrix_QR_disassemble(double A[DIMENSION][DIMENSION],

double Q[DIMENSION][DIMENSION],

double R[DIMENSION][DIMENSION], int M_Dimension);

//带双步位移QR方法求解矩阵A的全部特征值所用子程序

int checkPricision(double number);

int doubleDisplacement_QR(double A[DIMENSION][DIMENSION], double Lamta[DIMENSION][2]);

void getEigenVector(double A[DIMENSION][DIMENSION], double Lamta);

int rmain_Guass(int dimension, double coefficient[DIMENSION][DIMENSION],

double rightValue[DIMENSION], double X[DIMENSION]);

//主程序

int main()

{

double A[DIMENSION][DIMENSION], Q[DIMENSION][DIMENSION],

R[DIMENSION][DIMENSION], RQ[DIMENSION][DIMENSION];

double Lamta[DIMENSION][2];

int i, j, p;

//初始化矩阵A和特征值向量

intialMatrixA(A);

intialMatrix_Lamta(Lamta);

//对矩阵A进行拟上三角化

simulateUpTriangel(A);

printf("\nThe A's simulateUpTriangel Matrix is: \n");

printMatrix(A);

//对矩阵A进行QR分解，并求解RQ

matrix_QR_disassemble(A, Q, R, DIMENSION);

printf("\nThe Q Matrix is: \n");

printMatrix(Q);

printf("\nThe R Matrix is: \n");

printMatrix(R);

printf("\nThe RQ Matrix is: \n");

//求解矩阵RQ

for(i=0; i

{

for(j=0; j

{

RQ[i][j] = 0;

for(p=i; p

{

RQ[i][j] += R[i][p] * Q[p][j];

}

}

}

printMatrix(RQ);

//带双步位移QR分解求解矩阵A的所有特征值，存入特征值数组

intialMatrixA(A);

doubleDisplacement_QR(A, Lamta);

printEigenValue(Lamta);

//求解实特征值所对应的特征向量

for(i=0; i

{

if(Lamta[i][1]==0)

{

intialMatrixA(A);

getEigenVector(A, Lamta[i][0]);

}

}

return 0;

}

//初始化矩阵A

void intialMatrixA(double A[DIMENSION][DIMENSION])

{

int i, j;

for(i=1; i<=DIMENSION; i++)

{

for(j=1; j<=DIMENSION; j++)

{

if(i == j)

{

A[i-1][j-1] = 1.5*cos(i+1.2*j);

}

else

{

A[i-1][j-1] = sin(0.5*i+0.2*j);

}

}

}

}

//初始化特征值数组

void intialMatrix_Lamta(double Lamta[DIMENSION][2])

{

int i;

for(i=0; i

{

Lamta[i][0]=0;

Lamta[i][1]=0;

}

}

//打印矩阵元素

void printMatrix(double matrix[DIMENSION][DIMENSION])

{

int i, j;

//printf("\nThe Matrix is: \n");

for(i=0; i

{

for(j=0; j

{

printf("%.12e\t", matrix[i][j]);

}

printf("\n");

}

printf("\n");

}

//打印显示特征值

void printEigenValue(double Lamta[DIMENSION][2])

{

int i;

for(i=0; i

{

printf("EigenValue[%d]: Re- %.12e\tIm-i%.12e\n", i, Lamta[i][0], Lamta[i][1]);

}

}

//将Q矩阵初始化成单位阵

void intialMatrixQ_toUnitE(double Q[DIMENSION][DIMENSION])

{

int i, j;

for(i=0; i

{

for(j=0; j

{

if(i==j)

{

Q[i][j]=1;

}

else

{

Q[i][j]=0;

}

}

}

}

//将矩阵R初始化成0矩阵

void intialMatrixR_toZero(double R[DIMENSION][DIMENSION])

{

int i, j;

for(i=0; i

{

for(j=0; j

{

R[i][j]=0;

}

}

}

//判断元素是否全为0，实则返回真，否则返回假

int allIsZero(double A[DIMENSION][DIMENSION], int rowStartIndex,

int rowEndIndex, int clomn)

{

int i;

for(i=rowStartIndex; i<=rowEndIndex; i++)

{

if(A[i][clomn] != 0)

{

return Failure;

}

}

return Success;

}

//计算Dr

double caculate_dr(double A[DIMENSION][DIMENSION],

int rowStart, int rowEnd, int clomn)

{

int i;

double sum=0;

for(i=rowStart; i<=rowEnd; i++)

{

sum = sum + A[i][clomn]*A[i][clomn];

}

return sqrt(sum);

}

//取符号函数

int sgn(double a)

{

return a>0 ? 1 : -1;

}

//得到向量U

void getArray_U(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double c, int index)

{

int i;

for(i=0; i

{

U[i] = 0;

}

U[index] = A[index][index-1] - c;

for(i=index+1; i

{

U[i] = A[i][index-1];

}

}

//得到向量P和Q

void getArray_P_Q(double A[DIMENSION][DIMENSION], double U[DIMENSION],

double h, double P[DIMENSION], double Q[DIMENSION])

{

int i, j;

double sum[2]={0, 0};

for(i=0; i

{

sum[0] = sum[1] = 0;

for(j=0; j

{

sum[0] = sum[0] + A[j][i] * U[j];

sum[1] = sum[1] + A[i][j] * U[j];

}

P[i] = sum[0] / h;

Q[i] = sum[1] / h;

}

}

//求t

double get_T(double P[DIMENSION], double U[DIMENSION], double h)

{

int i;

double sum = 0;

for(i=0; i

>>>>>>内容过长，仅展示头部和尾部部分文字预览，全文请查看图片预览。<<<<<<

Lamta= 1.***7113e+000 's EigenVetor:

6.2***e-001

-1.1***e-001

-2.***0842e+000

-1.***0439e+000

-3.8***e+000

8.1***e+000

-1.***3690e+000

-6.***6184e-001

2.***4883e+000

1.***0000e+000

Lamta= 3.***7436e+000 's EigenVetor:

-4.***2970e-001

-9.***4238e-001

-1.***8327e+000

-1.***7819e+000

-1.***3865e+000

-1.***4413e+000

3.***8697e-001

1.***7744e+000

2.***1308e+000

1.***0000e+000

[文章尾部最后500字内容到此结束,中间部分内容请查看底下的图片预览]

以上为《数值分析大作业》的无排版文字预览，完整格式请下载

下载前请仔细阅读上面文字预览以及下方图片预览。图片预览是什么样的，下载的文档就是什么样的。