北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案：
按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。
总体思路：
1）初始化矩阵
首先需要将需要求解的矩阵输入程序。为了防止矩阵在后面的计算中被破坏保存A[][]。
2）对给定的矩阵进行拟上三角化
为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。这里用函数
void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。
3）对拟上三角化后的矩阵进行QR分解
对拟上三角化的矩阵进行QR分解会大大减小计算量。用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。
4）对拟上三角化后的矩阵进行带双步位移的QR分解
为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。这里用函数intTwoStepDisplacement_QR()来实现。这是用来存储计算得到的特征值的二维数组。考虑到特征值可能为复数，因此将所有特征值均当成复数处理。此函数中，QR分解部分用子函数void QR_decompositionMk()实现。这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。
5）计算特征向量
得到特征值后，计算实特征值相应的特征向量。这里判断所得特征值的虚数部分是否为零。求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。因此先初始化矩阵Array，计算(A-λI)，再求解方程组。
方程组的求解采用列主元的高斯消去法，由函数intGauss_column(double a[][N],double b[],double X[])实现。由于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。首先用高斯消去将矩阵(A-λI)化为上三角阵，其最后一行必定全为零。因此在反代的过程中，只要令x[]的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。
6）输出有关结果
根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。程序中构造了两个输出函数专门来解决不同输出结果的问题，void print_lambda(complex lambda[]);void print_matrix(double mat[][N])。
二、程序源代码：
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
#define L 100 //定义最大迭代次数
#define N 10 //定义矩阵大小
#define err 1e-12 //定义误差限
//定义一个结构体来表示复数
typedefstruct complex{
double rpart;
double ipart;
};
FILE * pReadFile;
//主函数
int _tmain(intargc, _TCHAR* argv[])
{
inti,j,t;
double y[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N];
struct complex lambda[N];
//声明要调用的函数
void QuasiTriangularization(double a[][N]);
void QR_decomposition(double A[][N],double Q[][N],double R[][N]);
void QR_decompositionMk(double mk[][N],int m, double a[][N]);
void print_lambda(complex lambda[]);
void print_matrix(double mat[][N]);
void multi_matrix(double a[][N],double b[][N],double c[][N]);
intTwoStepDisplacement_QR(double a[][N],complex lambda[]);
intGauss_column(double a[][N],double b[],double X[]);
//矩阵初始化
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j){
a[i][j] = 1.5 * cos((i+1) + 1.2 * (j+1));
A[i][j] = a[i][j];
}
else{
a[i][j] = sin(0.5 * (i+1) + 0.2 * (j+1));
A[i][j] = a[i][j];
}
}
}
for (i = 0; i < N; i++){
y[i] = 0;
}
//对矩阵a[][]拟上三角化
QuasiTriangularization(a);
//打开文件result.txt
pReadFile = fopen( "result.txt", "w" );
//打印结果到文件result.txt中
fprintf(pReadFile,"拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
//printf("拟上三角化之后的矩阵A[%d][%d]:\n",N,N);
print_matrix(a);
//对拟上三角化后的矩阵a[][]，进行QR分解
QR_decomposition(a,Q,R);
fprintf(pReadFile,"Q[%d][%d]:\n",N,N);
//printf("Q[%d][%d]:\n",N,N);
print_matrix(Q);
fprintf(pReadFile,"R[%d][%d]:\n",N,N);
//printf("R[%d][%d]:\n",N,N);
print_matrix(R);
multi_matrix(R,Q,RQ);
fprintf(pReadFile,"RQ[%d][%d]:\n",N,N);
//printf("RQ[%d][%d]:\n",N,N);
print_matrix(RQ);
//双步位移QR分解求全部特征值
TwoStepDisplacement_QR(a,lambda);
//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量
for (t = 0; t < N; t++){
if (lambda[t].ipart == 0){
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
if (i == j)
B[i][j] = A[i][j] - lambda[1].rpart;
else
B[i][j] = A[i][j];
}
}
Gauss_column(B,y,X);
fprintf(pReadFile,"\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
//printf("\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n",lambda[t].rpart,N);
for (i = 0; i < N; i++){
fprintf(pReadFile,"%.12e ",i,X[i]);
//printf("X[%d]: %.12e ",i,X[i]);
}
}
}
fclose( pReadFile );
scanf("%d",&i);
return 0;
}//主函数
/************************************* 矩阵的拟上三角化
输入n阶方阵a[][]，将a[][]拟上三角化
无返回值
**************************************/ void QuasiTriangularization(double a[][N]){
intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],qr[N],wr[N];
for (r = 0; r < N-2; r++){
//判断a[i][r](i=r+2,r+3,...,n-2)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+2; i < N; i++){
sum = sum || a[i][r];
}
//如果不是全部都是零，则计算
if (sum){
//计算dr
sum = 0;
for (i = r+1; i < N; i++){
sum += a[i][r] * a[i][r];
}
dr = sqrt(sum);
//计算cr
if (a[r+1][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * a[r+1][r];
//取向量ur[]
for (i = 0; i < N; i++){
if (i < r+1)
ur[i] = 0;
else
if (i == r+1)
ur[i] = a[i][r] - cr;
else
ur[i] = a[i][r];
}
//计算向量qr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[i][j] * ur[j];
qr[i] = sum / hr;
}
//计算向量pr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r+1; j < N; j++)
sum += a[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算tr
sum = 0;
for (i = r+1; i < N; i++)
sum += pr[i] * ur[i];
tr = sum / hr;
//计算wr[]
for (i = 0; i < N; i++){
if (i < r+1)
wr[i] = qr[i];
else
wr[i] = qr[i] - tr * ur[i];
}
//计算新产生的矩阵a[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
a[i][j] = a[i][j] - wr[i] * ur[j] - ur[i] * pr[j];
}
}
}
/*************************************
矩阵的QR分解
将A[][]分解为Q[][]*R[][]
无返回值
**************************************/
void QR_decomposition(double A[][N],double Q[][N],double R[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],wr[N];
//取矩阵R1[][]为A[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
R[i][j] = A[i][j];
//取矩阵Q1[][]为单位矩阵
for (i = 0; i < N; i++)
for (j = 0; j < N; j++){
if (i == j)
Q[i][j] = 1;
else
Q[i][j] = 0;
}
//
for (r = 0; r < N-1; r++){
//判断R[i][r](i=r+1,r+2,...,N-1)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+1; i < N; i++){
sum = sum || R[i][r];
}
if (sum){
//计算dr
sum = 0;
for (i = r; i < N; i++){
sum += R[i][r] * R[i][r];
}
dr = sqrt(sum);
//计算cr
if (R[r][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * R[r][r];
//取向量ur[]
for (i = 0; i < N; i++){
if (i < r)
ur[i] = 0;
else
if (i == r)
ur[i] = R[i][r] - cr;
else
ur[i] = R[i][r];
}
//计算向量wr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r; j < N; j++)
sum += Q[i][j] * ur[j];
wr[i] = sum;
}
//计算新的矩阵Qr+1[][],存储在Q[][]里面
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
Q[i][j] = Q[i][j] - wr[i] * ur[j] / hr;
//计算向量pr[]
for (i = 0; i < N; i++){
sum = 0;
for (j = r; j < N; j++)
sum += R[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算新产生的R[][]
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
R[i][j] = R[i][j] - ur[i] * pr[j];
}
}
}
/*************************************
Mk[][]矩阵的QR分解及求Ak+1[][]的计算
形参m为每次降阶之后的值
无返回值
**************************************/
void QR_decompositionMk(double mk[][N],int m, double a[][N]){ intr,i,j;
double tr,hr,cr,dr,sum;
double ur[N],pr[N],qr[N],wr[N],vr[N],br[N][N],Cr[N][N];
//取矩阵br[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
br[i][j] = mk[i][j];
//取矩阵Cr[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
Cr[i][j] = a[i][j];
//
for (r = 0; r <= m-1; r++){
//判断br[i][r](i=r+1,r+2,...,m)是否全为零
sum = 0;//变量sum使用前清零
for (i = r+1; i <= m; i++){
sum = sum || br[i][r];
}
if (sum){
//计算dr
sum = 0;
for (i = r; i <= m; i++){
sum += br[i][r] * br[i][r];
}
dr = sqrt(sum);
//计算cr
if (br[r][r] > 0)
cr = -dr;
else
cr = dr;
//计算hr
hr = cr * cr - cr * br[r][r];
//取向量ur[]
for (i = 0; i <= m; i++){
if (i < r)
ur[i] = 0;
else
if (i == r)
ur[i] = br[i][r] - cr;
else
ur[i] = br[i][r];
}
//计算向量vr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += br[j][i] * ur[j];
vr[i] = sum / hr;
}
//计算新的矩阵Br+1[][],存储在Br[][]里面
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
br[i][j] = br[i][j] - ur[i] * vr[j];
//计算向量pr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += Cr[j][i] * ur[j];
pr[i] = sum / hr;
}
//计算向量qr[]
for (i = 0; i <= m; i++){
sum = 0;
for (j = r; j <= m; j++)
sum += Cr[i][j] * ur[j];
qr[i] = sum / hr;
}
//计算tr
sum = 0;
for (i = r; i <= m; i++)
sum += pr[i] * ur[i];
tr = sum / hr;
//计算向量wr[]
for (i = 0; i <= m; i++){
if (i < r)
wr[i] = qr[i];
else
wr[i] = qr[i] - tr * ur[i];
}
//计算新产生的矩阵Cr[][]
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
Cr[i][j] = Cr[i][j] - wr[i] * ur[j] - ur[i] * pr[j]; }
}
//将计算出的Cr[][]矩阵中的值赋给矩阵A[][]，得到新的矩阵。
for (i = 0; i <= m; i++)
for (j = 0; j <= m; j++)
a[i][j] = Cr[i][j];
}
/*******************************************
高斯列主元素消元法求矩阵实特征值的特征向量
等到的特征向量放在X[]里面
返回值为整型
********************************************/ intGauss_column(double a[][N],double b[],double X[])
{
inti,j,k,row_no;
double temp,a_ik,m_ik,sum = 0;
//消元过程
for (k = 0; k < N - 1; k++){
//选行号，并记录
row_no = k;
a_ik = abs(a[k][k]);
for(i = k + 1; i < N; i++){
if (abs(a[i][k]) >a_ik){
a_ik = abs(a[i][k]);
row_no = i;
}
}
//交换刚刚选择的最大行与第K行的所有元素
if (row_no != k){//如果不是第K行的元素最大，则交换行for (j = k; j < N; j++){
temp = a[k][j];
a[k][j] = a[row_no][j];
a[row_no][j] = temp;
}
temp = b[k];
b[k] = b[row_no];
b[row_no] = temp;
}
if(a[k][k] == 0){
printf("error");
return(0);
}
for(i = k+1; i < N; i++){
m_ik = a[i][k] / a[k][k];
a[i][k] = 0;
for(j = k+1; j < N; j++){
a[i][j] = a[i][j] - m_ik * a[k][j];
}
b[i] = b[i] - m_ik * b[k];
}
}
//回代过程，由于a的行列式值为零，所以有无穷多个解。所以不妨设x[N-1] := 1 X[N-1] = 1;//b[N-1] / a[N-1][N-1];
for(k = N-2; k >= 0; k--){
sum = 0;
for(j = k+1; j < N; j++)
sum += X[j] * a[k][j];
X[k] = (b[k] - sum) / a[k][k];
}
return 0;
}
/*******************************************
打印特征向量的值到result文件中
无返回值
********************************************/
void print_lambda(complex lambda[]){
int i;
fprintf(pReadFile,"已经求出所有的%d个特征值λ：\n",N);
//printf("已经求出所有的%d个特征值λ：\n",N);
for (i = 0; i < N; i++){
if (lambda[i].ipart != 0){
if (lambda[i].ipart> 0)
fprintf(pReadFile,"%.12e + %.12e i\n",lambda[i].rpart,lambda[i].ipart);
//printf("%.12e + %.12e i\n",lambda[i].rpart,lambda[i].ipart);
else
fprintf(pReadFile,"%.12e + %.12e i\n",lambda[i].rpart,-lambda[i].ipart);
//printf("%.12e - %.12e i\n",lambda[i].rpart,-lambda[i].ipart);
}
else
fprintf(pReadFile,"%.12e\n",lambda[i].rpart);
//printf("%.12e\n",lambda[i].rpart);
}
fprintf(pReadFile,"\n");
//printf("\n");
}
/*******************************************
打印参量传入矩阵到result文件中
无返回值
********************************************/
void print_matrix(double mat[][N]){
inti,j;
for (i = 0; i < N; i++){
for (j = 0; j < N; j++)
fprintf(pReadFile,"%.12e ",mat[i][j]);
//printf("%.12e ",mat[i][j]);
//printf("\n");
fprintf(pReadFile,"\n");
}
}
/*******************************************
带双步位移的QR分解
将所有求出的特征值传入参量lambda[]中
返回值为整型
********************************************/ intTwoStepDisplacement_QR(double a[][N],complex lambda[]){
inti,j,m,l,flag = 1,count = 0,k = 0;
double mk[N][N],detD,b,s,t,sum,delta;
struct complex s1,s2;
m = N - 1;
//(3)到(9)的循环
while (1)
{
k++;
//(3)到(7)的循环
while (1)
{
//(3)到(4)的循环
while (1)
{
if (flag){
//(3)
if (abs(a[m][m-1]) <= err){
lambda[count].rpart = a[m][m];
lambda[count].ipart = 0;
count++;
m--;
}
else
break;
}
flag = 1;
//(4)
if (m == 0){
lambda[count].rpart = a[m][m];
lambda[count].ipart = 0;
print_lambda(lambda);
return 0;
}
else
if (m == -1){
print_lambda(lambda);
return 0;
}
}
//(5)
detD = a[m-1][m-1] * a[m][m] - a[m-1][m] * a[m][m-1];
b = -(a[m-1][m-1] + a[m][m]);
delta = b * b - 4 * detD;
if (delta >= 0){
s1.rpart = -(b + sqrt(delta)) / 2;
s1.ipart = 0;
s2.rpart = -(b - sqrt(delta)) / 2;
s2.ipart = 0;
}
else{
s1.rpart = -b / 2;
s1.ipart = -sqrt(-delta) / 2;
s2.rpart = -b / 2;
s2.ipart = sqrt(-delta) / 2;
}
//(6)
if (m == 1){
lambda[count] = s1;
count++;
lambda[count] = s2;
print_lambda(lambda);
return 0;
}
//(7)
if (abs(a[m-1][m-2]) <= err){
m = m - 2;
lambda[count] = s1;
count++;
lambda[count] = s2;
count++;
}
else
break;
flag = 0;
}
//(8)
if (k == L){
fprintf(pReadFile,"在规定的迭代次数L=%d 内未能得到全部特征值",L);
//printf("在规定的迭代次数L=%d 内未能得到全部特征值",L);
//scanf("%d",&i);
return 0;
}
//(9)
s = a[m-1][m-1] + a[m][m];
t = a[m-1][m-1] * a[m][m] - a[m][m-1] * a[m-1][m];
//计算矩阵Mk
for (i = 0; i <= m; i ++){
for (j = 0; j <= m; j++){
sum = 0;
for (l = 0; l <= m; l++)
sum += a[i][l] * a[l][j];
if (i == j)
mk[i][j] = sum - s * a[i][j] + t;
else
mk[i][j] = sum - s * a[i][j];
}
}
QR_decompositionMk(mk,m,a);
}
}
/*******************************************
矩阵乘法
a[][]与b[][]的乘积传入c[][]中
无返回值
********************************************/ void multi_matrix(double a[][N],double b[][N],double c[][N]){ inti,j,t;
double sum;
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
sum = 0;
for (t = 0; t < N; t++)
sum += a[i][t] * b[t][j];
c[i][j] = sum;
}
}
}
三、运算结果：
1.拟上三角化之后的矩阵A[10][10]:
-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000 -7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000 -8.782436382927e-002 -9.255889387184e-001 6.032599440534e-001 1.518860956469e-001 -2.347878362416e+000 2.372370104937e+000 1.819290822208e+000 3.237804101546e-001 2.205798440320e-001
2.102692662546e+000 1.816138086098e-001 1.278839089990e+000 -6.380578124405e-001 -4.154075603804e-001
0.000000000000e+000 1.728274599968e+000 -1.171467642785e+000 -1.243839262699e+000 -6.399758341743e-001
-2.002833079037e+000 2.924947206124e-001 -6.412830068395e-001 9.783997621285e-002 2.557763574160e-001
0.000000000000e+000 0.000000000000e+000 -1.291669534130e+000 -1.111603513396e+000 1.171346824096e+000
-1.307356030021e+000 1.803699177750e-001 -4.246385358369e-001 7.988955239304e-002 1.608819928069e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 1.560126298527e+000 8.125049397515e-001
4.421756832923e-001 -3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 -7.707773755194e-001
-1.583051425742e+000 -3.042843176799e-001 2.528712446035e-001 -6.709925401449e-001 2.544619929082e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
-7.463453456938e-001 -2.708365157019e-002 -9.486521893682e-001 1.195871081495e-001 1.929265617952e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+000 -7.701801374364e-001 -4.697623990618e-001 4.988259468008e-001 1.137691603776e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0007.013167092107e-001 1.582180688475e-001 3.862594614233e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+0004.843807602783e-001 3.992777995177e-001
2.QR分解之后拟上三角化之后的矩阵A[10][10]的结果Q[10][10]:
-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-001 6.486200753651e-002 3.709718861896e-001
1.855847143605e-001 -1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002
-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-002 6.124472142963e-003 4.440505443139e-002 1.975907909728e-002 2.062836970533e-002
0.000000000000e+000 -8.810520554692e-001 -3.989762796959e-001 3.720308728479e-002 2.127794064090e-001
1.064463557221e-001 -9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002
0.000000000000e+0000.000000000000e+000 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001
-3.533456368496e-001 3.101438948264e-002 2.248676491598e-001 1.000601783527e-001 1.044622748702e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 9.892235468621e-001 -1.239411731211e-001
-6.200358589825e-002 5.442272839461e-003 3.945881637235e-002 1.755813350011e-002 1.833059462907e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 5.323610690264e-001
-6.733900344896e-001 5.910581205868e-002 4.285425323867e-001 1.906901343193e-001 1.990794495295e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
-6.0597********e-001 -9.165783032819e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+000 9.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002 -4.501694411183e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+000 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+000 -7.221591336735e-001 6.917269588876e-001
3.QR分解之后拟上三角化之后的矩阵A[10][10]的结果R[10][10]:
2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000 -
3.558787493835e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-001 3.353802899665e-001
0.000000000000e+000 -1.961603277854e+000 2.407523727633e-001 7.054714572823e-001 5.957204318279e-001 5.526978774676e-001 -3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002
0.000000000000e+000 0.000000000000e+000 2.404534601993e+000 1.706758096328e+000 -4.239566704091e-001 3.405332305815e+000 -1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-001
0.000000000000e+000 0.000000000000e+0000.000000000000e+000 1.577122080722e+000 6.399535133956e-001 3.468127872427e-001 -5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-002
0.000000000000e+000 0.000000000000e+0000.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000 -2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-002 1.996629079956e-001
0.000000000000e+000 0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
1.231641451542e+000 1.619701003419e-001 1.962638275504e-001 5.350035621760e-001 -1.509273424767e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000-7.753441914209e-00
1 -3.464514508821e-001 4.312226803504e-001 1.234643696237e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00
01.296312940612e+000 -4.288053318338e-001 2.737334158165e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00 00.000000000000e+000 -6.707396440648e-001 -4.842320121884e-001
0.000000000000e+000
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+00 00.000000000000e+0000.000000000000e+000 7.168323926323e-002
4.QR分解之后拟上三角化之后的矩阵A[10][10]的结果Q、R乘积矩阵为:
1.163074414164e+000
2.632670934508e+000 -1.772796003272e+000 -8.668899138521e-002
3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -
4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-001
1.836115060851e+000 1.144286420080e-001 -9.880381403133e-001 5.589725694767e-001 4.694190067101e-002 -
2.978478237007e-001 1.617130577649e-002 6.936977702522e-001 1.367670571405e-001 1.419099231519e-002
0.000000000000e+000 -2.118520153533e+000 -1.876189745783e+000 -5.407071940597e-001 1.171538359721e+000 -2.550323020223e+000 1.691577936540e+000 1.229951613262e+000 1.387947777212e+000 8.667502917242e-001
0.000000000000e+0000.000000000000e+000 -8.471995127808e-001 4.382910468318e-001 -1.008632199185e+000 -7.959374261495e-001 4.769258865577e-001 4.072683083890e-001 4.096390493527e-001 3.363378940862e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+000 -1.432244342447e+000 -5.742284908055e-001
1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000 6.556779598004e-001 -9.275250974463e-001 2.529079844053e-001 6.905949216976e-001 -2.374430675823e-002 -2.429781119781e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+000
4.698400884876e-001 -2.730776009527e-001 7.821296259798e-001 -9.580964936399e-002 7.846239841323e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000 e+000 1.287679058937e+000 -3.576058900348e-001 -4.116725408806e-003 3.914268216423e-001
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000 e+0000.000000000000e+000-3.628760503545e-001 7.398980975354e-0017.241608309576e-002
0.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000e+0000.000000000000
e+0000.000000000000e+0000.000000000000e+000 -5.176670596524e-002 4.958522909877e-002
5.矩阵的所有特征值λ为:
λ1 =5.650488993501e-002
λ2 =6.360627875745e-001
λ3 =9.355889078188e-001
λ4 =-9.805309562902e-001 + 1.139489127430e-001*i
λ5 =-9.805309562902e-001 + 1.139489127430e-001*i
λ6 =-1.484039822259e+000
λ7 =1.577548557113e+000
λ8 =-2.323496210212e+000 + 8.930405177200e-001*i
λ9 =-2.323496210212e+000 + 8.930405177200e-001*i
λ10=3.383039617436e+000
6.矩阵的所有实特征值对应的特征向量为:
矩阵与特征值λ1=5.650488993501e-002对应的特征向量为X1[10]：
-5.105003830620e+000 -4.886319842355e+000 9.505161576142e+000 -6.788331788207e-001 -9.604334110489e+000 -3.0457********e+000 1.574873885600e+001 -7.395037126270e+000 -7.109654943654e+000 1.000000000000e+000
矩阵与特征值λ2=6.360627875745e-001对应的特征向量为X2[10]：
4.745031936534e+000 3.157868541747e+000 1.729946912419e+001 -1.980049331453e+000
-3.187521973521e+001 7.794009603193e+000 -1.004255685844e+001 1.670757770479e+001
1.310524273052e+001 1.000000000000e+000
矩阵与特征值λ3=9.355889078188e-001对应的特征向量为X3[10]：
2.792418944526e+000 1.598236841510e+000 -5.207507040906e-001 -1.667886451700e+000
-1.225708535858e+001 7.241214790770e+000 -5.398214501428e+000 2.841008912972e+001 -1.216518754415e+0011.000000000000e+000
矩阵与特征值λ6=-1.484039822259e+000对应的特征向量为X6[10]：
-1.730784582110e+001 2.408192534965e+001 4.133852365114e-001 -8.572044074529e+000
9.287334657923e-002 -7.832726042765e-002 -6.374274025708e-001 -3.403204760828e-001
-3.784865409365e-001 1.000000000000e+000
矩阵与特征值λ7=1.577548557113e+000对应的特征向量为X7[10]：
6.217350824575e-001 -1.115111815225e-001 -2.483773580801e+000 -1.306860840419e+000
-3.815605442529e+000 8.117305509387e+000 -1.239170883674e+000 -6.800309586190e-001
2.691900144837e+000 1.000000000000e+000
矩阵与特征值λ10=3.383039617436e+000对应的特征向量为X10[10]：
-4.367398732970e-001 -9.063574014237e-001 -1.963340418326e+000 -1.082857447819e+000 -1.269297783865e+000 -1.0724********e+000 3.625108238697e-001 1.682905297743e+000
2.112744711308e+000 1.000000000000e+000
北航2010-2011年研究生数值分析期末模拟试卷1-3
数值分析模拟试卷1 一、填空（共30分，每空3分） 1 设??? ? ??-=1511A ，则A 的谱半径=)(a ρ______，A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________, ],,[321+++n n n n x x x x f ，＝________. 3 设?????≤≤-++≤≤+=2 1,121 0,)(2 323x cx bx x x x x x S ，是以0，1，2为节点的三次样条函数，则b=________,c=________. 4 设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x q ，则 ?=1 )(dx x xq k ________，=)(2 x q ________. 5 设???? ??????=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三角阵，当 其对角线元素)3,2,1(=i L ii 满足条件________时，这种分解是唯一的. 二、（14分）设4 9,1,41,)(2102 3 === =x x x x x f , （1）试求)(x f 在]4 9,41[上的三次Hermite 插值多项式)(x H 使满足 2,1,0),()(==i x f x H i i ，)()(11x f x H '='. （2）写出余项)()()(x H x f x R -=的表达式. 三、（14分）设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3 2 41+ =+， （1） 证明R x ∈?0均有? ∞ →=x x n x lim （? x 为方程的根）； （2） 取40=x ，用此迭代法求方程根的近似值，误差不超过，列出各次迭代值； （3）此迭代的收敛阶是多少？证明你的结论. 四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式 有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？
北航数值分析大作业一
《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求 出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include 《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i 得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-， 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include 数值分析第二次大作业 史立峰 SY1505327 一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。 北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日 一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。 《数值分析A》 一、算法设计方案 整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。 1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。 2.对经过拟上三角化处理的矩阵进行QR分解。 3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。 二、源程序代码 #include 习 题 1. 指出有效数49×102，0.0490，490.00的绝对误差限、相对误差限和有效数字位数． 2. 将 3.142作为π的近似值，它有几位有效数字，相对误差限和绝对误差限各为多少？ 3. 要使101的近似值x * 的相对误差限不超过4102 1?×，问查开方表时x * 需要保留几位有效数字？ 4. 已知近似数x * 有两位有效数字，试估计其相对误差限． 5. 设x * 为x 的近似数, 证明n x * 的相对误差大约为x * 相对误差的n 1倍． 6. 某矩形的长和宽大约为100cm 和50cm, 应该选用最小刻度为多少cm 的测量工具, 才能保证计算出的面积误差（绝对值）不超过0.15cm 2． 7. 已知三角形面积c ab S sin 2 1=，测量a , b , c 时产生的相对误差为)(*a e r ，)(*b e r ，)(*c e r ，其中2 ,0*π< 航空航天大学 《数值分析》计算实习报告 第一大题 学院：自动化科学与电气工程学院 专业：控制科学与工程 学生姓名： 学号： 教师： 电话： 完成日期： 2015年11月6日 航空航天大学 Beijing University of Aeronautics and Astronautics 实习题目： 第一题 设有501501?的实对称矩阵A ， ??? ???? ?????????=5011A a b c b c c b c b a 其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1 .0-==???=--=c b i e i i a i i 。矩阵A 的特征值为)501,,2,1(???=i i λ，并且有 ||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤ 1.求1λ，501λ和s λ的值。 2.求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 3.求A 的(谱数)条件数2)A (cond 和行列式detA 。 说明： 1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时，应使矩阵A 的所有零元素都不储存。 3.打印以下容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求1λ，501λ和s λ的值。 由于||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤，可知绝对值最大特征值必为1λ和501 λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则 501λ=λ。将矩阵A 进行一下平移： I -A A'λ= （1） 对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。 s λ为按模最小特征值，||min ||501 1i i s λλ≤≤=，可对A 使用反幂法求得。 2、求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,...,2,1(=k k i λ。 计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换： )39,,2,1k -A A k k ==（I μ （2） 对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。 3、求A 的(谱数)条件数2)(A cond 和行列式detA 。 由矩阵A 为非奇异对称矩阵可得： | | )(min max 2λλ=A cond （3） 其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。 在进行反幂法求解时，要对A 进行LU 分解得到。因L 为单位下三角阵，行 列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。 《数值分析A》计算实习题目第一题 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有 对角线上元素的乘积。 二．源程序(VS2010环境下，C++语言) #include 《数值分析》计算实习作业 （第二题） 算法设计方案： 1、对矩阵A 赋值，取计算精度ε=1×10-12； 2、对矩阵A 进行拟上三角化，得到A (n-1)，并输出A (n-1)； 对矩阵A 的拟上三角化，通过直接调用子函数inftrianglize(A)来实现；拟上三角化得到的矩阵A (n-1)输出至文件solution.txt 中。 3、对A (n-1)进行QR 分解并输出Q 、R 及RQ 矩阵； QR 分解通过直接调用子函数QRdescom(A,Q,R, n)实现。 4、运用QR 方法求所有的特征值，并输出； (1)初始时令m=n ，在m>2的条件下执行； (2)判断如果|A mm-1|<ε，则得到一个特征值，m=m-1，转(4)；否则转(3)； (3)判断如果|A m-1m-2|<ε，则得到两个特征值，m=m-2，转(4)； (4)判断如果m ≤2，转(6)；否则转(5)； (5)执行相似迭代，转(2)； k k T k k k k k k k k k k Q A Q A R Q M I D A D tr A M ==+-=+1)2)det(( (6)求出最后的一个或两个特征值； (7)输出全部的特征值至文件solution.txt 中。 5、输出QR 分解法迭代结束之后的A (n-1)至文件solution.txt 中； 6、通过反幂法求出所有实特征值的特征向量并输出。 首先令B=(A-λi I)，其中λi 是实特征值；反幂法通过调用子函数Bpowmethod(B,x1)实现，最终λi 对应的特征向量就是x1；最后将所有的实特征值的特征向量输出。 一、算法设计方案 1.使用牛顿迭代法，对原题中给出的i x i 08.0=，j y j 05.05.0+=， (010 ,020i j ≤≤≤≤)的11*21组j i y x ,分别求出原题中方程组的一组解，于是得到一组和i i y x ,对应的j i t u ,。 2.对于已求出的j i t u ,，使用分片二次代数插值法对原题中关于u t z ,,的数表进行插值得到 ij z 。于是产生了z=f(x,y)的11*21个数值解。 3.从k=1开始逐渐增大k 的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的σ,k 。当7 10-<σ时结束计算，输出拟合结果。 4.计算)5,,2,1,8,,2,1)(,(),,(* ***???=???=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0* *+==。 二、算法实现方案 1、求(,)f x y ： （1）Newton 法解非线性方程组 0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79 t u v w x t u v w y t u v w x t u v w y +++-=??+++-=? ? +++-=??+++-=?， 其中，t, u, v ,w 为待求的未知量，x, y 为代入的已知量。 设(,,,)T t u v w ξ=，给定精度水平12110ε-=和最大迭代次数M ，则解该线性方程组的迭代格式为： *(0)(0)(0)(0)(0)(k+1) ()()1()(,,,)()()0,1,T k k k t u v w F F k ξξξ ξξξ-?=?'=-??= ? 在附近选取初值， 迭代终止条件为()(1) () 1/k k k ξξ ξε-∞ ∞ -≤，若k M >时仍未达到迭代精度，则迭代计算失 败。 其中，雅可比矩阵 0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u) + v + w - y - 1.07()0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79F ξ???? ? ?=?????? ， 北航数值分析计算实习 报告一 Company number：【0089WT-8898YT-W8CCB-BUUT-202108】 北京航空航天大学 《数值分析》计算实习报告 第一大题 学 院：自动化科学与电气工程学院 专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日 北京航空航天大学 Beijing University of Aeronautics and Astronautics 实习题目： 第一题 设有501501?的实对称矩阵A ， 其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==???=--=c b i e i i a i i 。矩阵A 的特征值为)501,,2,1(???=i i λ，并且有 1.求1λ，501λ和s λ的值。 2.求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,,2,1(???=k k i λ。 3.求A 的(谱范数)条件数2)A (cond 和行列式detA 。 说明： 1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。 2.选择算法时，应使矩阵A 的所有零元素都不储存。 3.打印以下内容： （1）全部源程序； （2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。 4.采用e 型输出实型数，并且至少显示12位有效数字。 一、算法设计方案 1、求1λ，501λ和s λ的值。 由于||min ||,501 150121i i s λλλλλ≤≤=≤???≤≤，可知绝对值最大特征值必为1λ和501λ其中之 一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。将矩阵A 进行一下平移： I -A A'λ= （1） 对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为 'λ+λ。 s λ为按模最小特征值，||min ||501 1i i s λλ≤≤=，可对A 使用反幂法求得。 2、求A 的与数40 1 5011λλλμ-+=k k 最接近的特征值)39,...,2,1(=k k i λ。 计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。因此对A 进行平移变换： )39,,2,1k -A A k k ==（I μ （2） 对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。 北京航空航天大学2020届研究生 《数值分析》实验作业 第九题 院系：xx学院 学号： 姓名： 2020年11月 Q9:方程组A.4 一、 算法设计方案 （一）总体思路 1.题目要求∑∑=== k i k j s r rs y x c y x p 00 ),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。 ),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。 2.),(* * j i y x f 与1使用相同方法求得，),(* * j i y x p 由计算得出的p(x,y)直接带入),(* * j i y x 求得。 1. ),(i i y x 与),(i i y x f 的数表的获得 对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。再将 ),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。 2.乘积型最小二乘曲面拟合 2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下： ????? ??=k i i k x x x x B 0000 ， ????? ??=k j j k y y y y G 0000 数表矩阵如下： ???? ? ? ?=),(),(),(),(0000j i i j y x f y x f y x f y x f U 记C=[rs c ]，则系数rs c 的表达式矩阵为： 11-)(-=G G UG B B B C T T T ）（ 通过求解如下线性方程，即可得到系数矩阵C 。 UG B G G C B B T T T =)()( 2.2计算),(),,(* ***j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值 ),(**j i y x f 的计算与),(j i y x f 相同。将),(**j i y x 代入原方程组，求解响应) ,(* *ij ij u t 进行分片双二次插值求得),(**j i y x f 。),(* *j i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。 二、 源程序 ********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, m dimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) dimension X1(8), Y1(5), FXY1(8, 5), PXY1(8, 5), UX1(8, 5), TY1(8, 5) 北航数值分析作业第一题： 一、算法设计方案 1.要求计算矩阵的最大最小特征值，通过幂法求得模最大的特征值，进行一定 判断即得所求结果； 2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝 对值满足反幂法的要求，故通过反幂法即可求得； 3.反幂法计算时需要方程求解中间过渡向量，需设计Doolite分解求解； 4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对角线相 乘即为矩阵的Det。 算法编译环境：vlsual c++6.0 需要编译函数：幂法，反幂法，Doolite分解及方程的求解 二、源程序如下： #include 《数值分析A》计算实习题目二 姓名 学号 数值分析第二次大作业 一、算法设计方案 首先构造矩阵A，利用Householder矩阵对矩阵A作相似变换，把A化为拟上三角矩阵A（n-1），算法如课本P61。 使用QR分解法对矩阵A（n-1）进行QR分解，算法如课本P59, 进而求出所得矩阵的Q、R、RQ矩阵。 然后对A（n-1）进行带双步位移的QR分解求矩阵的全部特征值，采用以下几步进行： 第一步：判断是否a m,m-1(k)<=ε ,若不是,则进入第四步。若是，则a m,m-1(k)为特征值，m=m-1，若m=1，则进入第二步，若m=2进入第三步，否则转第四步。 第二步：m=1，则a11(k)为特征值，转向结束步。 第三步：m=2,则可以求出A的两个特征值s1和s2，转向结束步。 第四步：判断是否a m-1,m-2(k)<=ε，若不是，进入第五步。若是，则得到A的两个特征值s1和s2，令m=m-2，若m=1，进入第二步，若m=2进入第三步，否则进入第一步。 第五步：判断是否达到循环上限，若达到，则结束，否则进入第六步。 第六步：对A进行双步位移QR分解，这里的算法如课本P64，分解后转向计数步。 计数步：对循环次数进行计数，并转向第一步。 结束步：显示所求得的特征值。 最后对实特征值利用列主元高斯消元法求解其对应的特征向量，算法如课本p17. 二、源程序代码 #include 数值分析第二题 梁进明 SY0906529 算法设计方案。 一．矩阵的QR 分解。把矩阵A 分解为一个正交矩阵Q 与一个上三角矩阵R 的乘积，称为矩阵A 的 正三角分解，简称QR 分解。QR 分解的算法如下： 记1A A =，并记[]r ij n n Ar a ?=，令1Q I =（n 阶单位矩阵） 对于r=1,2,…,n-1执行 (1) 若(1,2,...,)r ir a i r r n =++全为零，则令1r r Q Q +=，1r r A A +=转(5);否则转(2) (2) 计算 2r d = ()sgn()r r rr r c a d =-(若() 0r rr a =,则取r r c d =) 2() r r r r rr h c c a =- (3) 令()()()1,(0,...,0,,,...,)r r r T n r rr r r r nr u a c a a R +=-∈ (4) 计算 11//r r r T r r r r r T r r r r T r r r r Q u Q Q u h p A u h A A u p ωω++==-==- (5) 继续 当此算法执行完后就得到正交矩阵n Q Q =和上三角矩阵n R A =且有A QR =。 二．矩阵的 拟上三角化。对实矩阵A 的拟上三角化具体算法如下： 记(1) A A =，并记()r A 的第r 列到第n 列的元素为(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于1,2,...,2r n =-执行 （1） 若() (2,3,...,)r ir a i r r n =++全为零，则令(1) ()r r A A +=,转(5);否则转(2)。 （2） 计算北航数值分析大作业第一题幂法与反幂法
北航数值分析报告第三次大作业
北航数值分析大作业第二题精解
北航数值分析大作业第二题
北航数值分析报告大作业第八题
北航数值分析第二次大作业--QR分解
数值分析
北航数值分析计算实习报告一
北航数值分析课程第一次大作业讲解
北航数值分析大作业第二次
北航数值分析大作业3
北航数值分析计算实习报告一
BUAA数值分析大作业三
北航数值分析作业第一题题解
北航数值分析第二次大作业--QR分解
北航数值分析大作业 第二题 QR分解