北航研究生数值分析编程大作业1

数值分析大作业

一、算法设计方案

1、矩阵初始化

矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，

在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。这样所需要的存储单元数大大减少，从而极大提高了运算效率。

2、利用幂法求出5011λλ，

幂法迭代格式：

0111111n

k k k k k

k T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩

非零向量 当1210/-≤-k k βββ时，迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。

3、利用反幂法求出ik s λλ,

反幂法迭代格式：

0111111n

k k k k k

k T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩

非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：

（1）作分解LU A =

对于n k ,...,2,1=执行

[]

[]

s k n r k k k i c c c c c n s k k k j c c

c c k s k

s k t k s k r i t t s t i k s k i k s k i j

s j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1

（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y)

)

1,...,2,1(/)(:/:)

,...,3,2(:,1),min(1.1.11),1max(,1--=-

===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i n

s n n t

i r i t t s t i i i

使用反幂法，直接可以求得矩阵按模最小的特征值s λ。

求与数)39,...,2,1(401

5011=-+=k k k λλλμ最接近的特征值ik λ，对矩阵I

A k μ-实行反幂法，即可求出对应的k k ik k k μλλβλ+==,/1。

4、求出A 的条件数和行列式

根据max 2()s

cond A λλ=，其中分子分母分别对应按模最大和最小的特征值。 det()A 的计算：由于A LU =,其中L 为下三角矩阵，且对角线元素为1，故det()1L =，所以有A LU U ==，又U 为上三角矩阵，故det()U 为对其对角线上各元素的乘积，最后可得det()det()A U =。

二、程序源代码

（1）定义所需要的函数：

#include

#include

#include

#define N 501

#define R 2

#define S 2

int min(int a,int b); // 求最小值

int max(int a,int b,int c); // 求最大值

double Fan_two(double x[N]);//计算二范数

void FenjieLU(double (*C)[N]);//解线性方程组的LU分解过程

void Solve(double (*C)[N], double *b,double *x);//解线性方程组的求解过程double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//幂法

double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//反幂法

};

（2）程序的主函数，Main.cpp代码如下：

void main()

{

double C[R+S+1][N];

double u[N];

double y[N];

double miu[39];

double C1[R+S+1][N];

double bta = 1.0;

double Namda1,Namda501,NamdaS;

double Namda[39];

double CondA2;

double detA = 1.0;

double D = 1.0e-12;

int i, j, k;

FILE * fp;

fp = fopen("Namda.txt","w");

//对数组进行初始化//

int i, j;

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

{

u[i] = 1;