北航数值分析大作业第一题幂法与反幂法
北航数值分析大作业第一题幂法与反幂法
-标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
《数值分析》计算实习题目
第一题：
1. 算法设计方案
（1）1λ，501λ和s λ的值。
1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。
2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。
（2）与140
k λλμλ-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
#include
void init_a();//初始化A
double get_an_element(int,int);//取A 中的元素函数
double powermethod(double);//原点平移的幂法
double inversepowermethod(double);//原点平移的反幂法
int presolve(double);//三角LU分解
int solve(double [],double []);//解方程组
int max(int,int);
int min(int,int);
double (*u)[502]=new double[502][502];//上三角U数组
double (*l)[502]=new double[502][502];//单位下三角L数组
double a[6][502];//矩阵A
int main()
{
int i,k;
double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;
double lambda[40];
init_a();//初始化A
lambdat1=powermethod(0);
lambdat2=powermethod(lambdat1);
lambda1=lambdat1 lambda501=lambdat1>lambdat2?lambdat1:lambdat2; presolve(0); lambdas=inversepowermethod(0); det=1; for(i=1;i<=501;i++) det=det*u[i][i]; for (k=1;k<=39;k++) { mu[k]=lambda1+k*(lambda501-lambda1)/40; presolve(mu[k]); lambda[k]=inversepowermethod(mu[k]); } printf("------------所有特征值如下------------\n"); printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501); printf("λs=%1.11e\n",lambdas); printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas)); printf("detA=%1.11e \n",det); for (k=1;k<=39;k++) {printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u; delete []l;//释放堆内存 return 0; } void init_a()//初始化A { int i; for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064; for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16; for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); } double get_an_element(int i,int j)//从A中节省存储量的提取元素方法 { if (fabs(i-j)<=2) return a[i-j+3][j];else return 0; } double powermethod(double offset)//幂法 { int i,x1; double u[502],y[502]; double beta=0,prebeta=-1000,yita=0; for (i=1;i<=501;i++) u[i]=1,y[i]=0;//设置初始向量u[] for (int k=1;k<=10000;k++) { yita=0; for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]); for (i=1;i<=501;i++) y[i]=u[i]/yita; for (x1=1;x1<=501;x1++) {u[x1]=0;for (int x2=1;x2<=501;x2++) u[x1]=u[x1]+((x1==x2) (get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2];} prebeta=beta;beta=0; for (i=1;i<=501;i++) beta=beta+ y[i]*u[i]; if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta- beta)/beta),k);break;}; //输出中间过程,包括偏移量,误差,迭代次数 } return (beta+offset); } double inversepowermethod(double offset)//反幂法 { int i; double u[502],y[502]; double beta=0,prebeta=0,yita=0; for (i=1;i<=501;i++) u[i]=1,y[i]=0; //设置初始向量u[] for (int k=1;k<=10000;k++) { yita=0; for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]); for (i=1;i<=501;i++) y[i]=u[i]/yita; solve(u,y); prebeta=beta;beta=0; for (i=1;i<=501;i++) beta=beta+ y[i]*u[i]; beta=1/beta; if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;}; //输出中间过程,包括偏移量,误差,迭代次数 } return (beta+offset); } int presolve(double offset)//三角LU分解 { int i,k,j,t; double sum; for (k=1;k<=501;k++) for (j=1;j<=501;j++) { u[k][j]=l[k][j]=0; if (k==j) l[k][j]=1; } //初始化LU矩阵 for (k=1;k<=501;k++) { for (j=k;j<=min(k+2,501);j++) { sum=0; for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++) sum=sum+l[k][t]*u[t][j]; u[k][j]=((k==j)(get_an_element(k,j)-offset):get_an_element(k,j))-sum; } if (k==501) continue; for (i=k+1;i<=min(k+2,501);i++) { sum=0; for (t=max(1,max(i-2,k-2));t<=(k-1);t++) sum=sum+l[i][t]*u[t][k]; l[i][k]=(((i==k)(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k]; } } return 0; } int solve(double x[],double b[])//解方程组 { int i,t; double y[502]; double sum; y[1]=b[1]; for (i=2;i<=501;i++) { sum=0; for (t=max(1,i-2);t<=i-1;t++) sum=sum+l[i][t]*y[t]; y[i]=b[i]-sum; } x[501]=y[501]/u[501][501]; for (i=500;i>=1;i--) { sum=0; for (t=i+1;t<=min(i+2,501);t++) sum=sum+u[i][t]*x[t]; x[i]=(y[i]-sum)/u[i][i]; } return 0; } int max(int x,int y) 《数值分析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 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC 文件。 3）程序清单，生成M 文件。 解答： >> A=rand(5) %随机产生5*5矩阵 求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A 矩阵和A 的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=??? ???? ???? ?? ???6286.13929.03467.15979.08044.03929.00944.12144.18506.13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613.09173.04187.1 《数值分析》计算实习题目 第一题： 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 《数值分析》实验8 一．实验名称：反幂法 二、实验目的： （1） 掌握求矩阵按模最小特征值及其对应的特征向量的规范化反幂法； （2） 掌握原点移位法。 三、实验要求 （1） 按照题目要求完成实验内容 （2） 写出相应的实验原理与C 语言程序 （3） 给出实验结果，结果分析 （4） 写出相应的实验报告 四、实验题目： 用反幂法求以下矩阵的指定特征值及其特征向量（迭代终止误差为1e-3）： 41411014110?? ? ? ??? 接近2的特征值及其特征向量. 程序： #include 数值分析之幂法及反幂法C 语言程序实例 1、算法设计方案： ①求1λ、501λ和s λ的值： s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。 1λ、501λ：已知矩阵A 的特征值满足关系 1n λλ<< ，要求1λ、及501λ时，可 按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。 b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m B A I λ=+，对矩阵B 用反幂法 求得B 的按模最小特征值2m λ。 c ． 321m m m λλλ=- 则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。 ②求和A 的与数5011 140 k k λλμλ-=+最接近的特征值 ik λ（k=0，1，…39）： 求矩阵A 的特征值中与k μ最接近的特征值的大小，采用原点平移的方法： 先求矩阵 B=A-k μI 对应的按模最小特征值k β，则k β+k μ即为矩阵A 与k μ最接近的特征值。 重复以上过程39次即可求得ik λ（k=0，1，…39）的值。 ③求A 的（谱范数）条件数2cond()A 和行列式det A ： 在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。 求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()s cond A λλ= ，max λ和s λ分别为模最大特征值与模最小特征值。 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： 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 =-， 幂法和反幂法的matlab实现 幂法求矩阵主特征值及对应特征向量 摘要 矩阵特征值的数值算法，在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题，所以说研究利用数学软件解决求特征值的问题是非常必要的。实际问题中，有时需要的并不是所有的特征根，而是最大最小的实特征根。称模最大的特征根为主特征值。 幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，它最大的优点是方法简单，特别适用于大型稀疏矩阵，但有时收敛速度很慢。 用java来编写算法。这个程序主要分成了四个大部分：第一部分为将矩阵转化为线性方程组；第二部分为求特征向量的极大值；第三部分为求幂法函数块；第四部分为页面设计及事件处理。其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法，再经过一系列的验证和迭代得到结果。 关键字：主特征值；特征向量；线性方程组；幂法函数块 POWER METHOD FOR FINDING THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF THE MATRIX ABSTRACT Numerical algorithm for the eigenvalue of matrix, in science and engineering technology, a lot of problems in mathematics are attributed matrix characteristic value problem, so that studies using mathematical software to solve the eigenvalue problem is very necessary. In practical problems, sometimes need not all eigenvalues, but the maximum and minimum eigenvalue of real. The characteristic value of the largest eigenvalue of the modulus maximum. Power method is a calculation of main features of the matrix values (matrix according to the characteristics of the largest value) and the corresponding eigenvector of iterative method. It is the biggest advantage is simple method, especially for large sparse matrix, but sometimes the convergence speed is very slow. Using java to write algorithms. This program is divided into three parts: the first part is the matrix is transformed into linear equations; the second part for the sake of feature vector of the maximum; the third part is 1．实验目的： 1熟练掌握C 语言程序设计，编程求解问题。 2.运用幂法求解住特征值和特征向量。 2．实验内容： 例题： 用幂法求 A= ??????????0.225.05.025.00.10.15.00.10.1 的特征值和特征向量。 完整代码以及截图如下: #include "stdio.h" #include "math.h" #define M 3 void main() { float fan(),max(),e1,e2,r1,r2; void au(),ex(),print_x(),std(); static float a[M][M]={{1.0,1.0,0.5},{1.0,1.0,0.25},{0.5,0.25,2.0}}; static float u0[M],u1[M],maxn0,maxn1; int i; printf("*********************************\n"); printf("****** 幂法*********\n"); printf("******求特征值与特征向量*********\n"); printf("*********************************\n\n"); printf("input precision e1,e2:"); scanf("%f,%f",&e1,&e2); printf("\ninput u(%d):",M); for (i=0;i 幂法反幂法求解矩阵大小特征值及其对应的特征向量 ————————————————————————————————作者：————————————————————————————————日期： 数值计算解矩阵的按模最大最小特征值及对应的特征向量 一.幂法 1. 幂法简介： 当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。矩阵A 需要满足的条件为： (1) 的特征值为A i n λλλλ,0||...||||21 ≥≥≥> (2) 存在n 个线性无关的特征向量，设为n x x x ,...,,21 1.1计算过程： i n i i i u x x αα,1 ) 0()0(∑==，有对任意向量不全为0，则有 1 11111221 12111 1 1 11 1 011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k n i i k i i n i i i k )(k (k))(k αλλλλλα++++=+=+++≈? ? ????+++======∑∑ 可见，当||1 2 λλ越小时，收敛越快；且当k 充分大时，有1)11 11)11111λαλαλ=??????==+++(k )(k k (k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。 2 算法实现 . ,, 3,,1 , ).5() 5(,,,,||).4();max(,).3() (max(;0,1).2(,).1()() () (停机否则输出失败信息转置若转否则输出若计算最大迭代次数，误差限，初始向量输入矩阵βλβεβλβλε←+←<<-←←= ←←k k N k y x Ay x x abs x y k N x A k k k 3 matlab 程序代码 数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1) 目标：使用带双步位移的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**+==。 数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序 矩阵的特征值与特征向量的计算 摘要 物理，力学，工程技术中的很多问题在数学上都归结于求矩阵特征值的问题，例如振动问题（桥梁的振动，机械的振动，电磁振动等）、物理学中某些临界值的确定问题以及理论物理中的一些问题。矩阵特征值的计算在矩阵计算中是一个很重要的部分，本文使用幂法和反幂法分别求矩阵的按模最大，按模最小特征向量及对应的特征值。 幂法是一种计算矩阵主特征值的一种迭代法，它最大的优点是方法简单，对于稀疏矩阵比较合适，但有时收敛速度很慢。其基本思想是任取一个非零的初始向量。由所求矩阵构造一向量序列。再经过所构造的向量序列求出特征值和特征向量。 反幂法用来计算矩阵按模最小特征向量及其特征值，及计算对应于一个给定近似特征值的特征向量。本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。然后经过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。 关键词：矩阵；特征值；特征向量；冥法；反冥法 THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIX ABSTRACT Physics, mechanics, engineering technology in a lot of problems in mathematics are attributed to matrix eigenvalue problem, such as vibration (vibration of the bridge, mechanical vibration, electromagnetic vibration, etc.) in physics, some critical values determine problems and 一、问题的描述及算法设计 （一）问题的描述 我所要做的课题是：对称矩阵的条件数的求解设计 1、求矩阵A 的二条件数 问题 A=?? ?? ? ?????----210121012 2、设计内容： 1）采用幂法求出A 的 错误!未找到引用源。. 2）采用反幂法求出A 的错误!未找到引用源。. 3）计算A 的条件数 ⅡA Ⅱ2* ⅡA -1Ⅱ2=cond2（A ）=错误!未找到引用源。/错误!未找到引用源。.（精度要求为10-6） 3、设计要求 1）求出ⅡA Ⅱ2。 2）并进行一定的理论分析。 （二）算法设计 1、幂法算法 （1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1. （2）计算 v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k （3）若| m k = m 1-k |<ε，则停止计算（m k 作为绝对值最大特征值1λ，u )(k 作为相应的特征向量）否则置k=k+1，转（2） 2、反幂法算法 （1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1. （2）对A 作LU 分解，即A=LU （3）解线性方程组 Ly )(k =u )1(-k ,Uv )(k =y )(k （4）计算 m k =max(v )(k ), u )(k = v )(k / m k （5）若|m k =m 1-k |<ε，则停止计算（1/m k 作为绝对值最小特征值n λ，u )(k 作 为相应的特征向量）;否则置k=k+1，转（3）. 二、算法的流程图（一）幂法算法的流程图 《数值分析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 当前位置：第7章>>第1节>>7.1.3 逆幂法 逆幂法是求实方阵按模最小特征值及相应的特征向量的一种反迭代方法。 1. 求A按模最小的特征值 设非奇异矩阵A的n 个特征值为，其相应的特征向量为e ，则的特征值为 其相应的特征向量仍为。 按模最大的特征值的倒数则为矩阵A按模最小的特征值。 利用乘幂法求按模最大的特征值。 任取初始非零初始向量，作迭代序列 它等价于 （7.5） 我们可以通过反迭代过程，即解方程组. 求得. 当k 充分大时，则有 在实际计算中，为了减少运算量，先将矩阵A作三角分解A=LR 然后再求解方程组 2．求在附近的特征值 设与最接近的特征值为即有 作矩阵，它的特征值和相应的特征向量为 若用逆幂法于矩阵，则有 则可求出矩阵的按模最小的特征值和相应的特征向量为 于是得A在附近的特征值和相应的特征向量为 （7.6） 例3 用逆幂法求矩阵在3.4附近的特征值和相应的特征向量 解对进行三角分解得： 用半次迭代法，取，则 得 再解 得 再解 得 于是 练习7.1 1.用乘幂法求矩阵按模最大特征值与特征向量 . 乘幂法的计算公式 设矩阵A的n个特征值按模的大小排列为： 其相应的特征向量为且它们是线性无关的。 先任取非零初始向量，作迭代序列 首先将表示为 所以 为了得出计算和的公式，下面分三种情况讨论 1．为实根，且。 当不为0，k充分大时，则有 所以（7.1）2．为实根，且。 当不为0，k充分大时，则有 （7.2）于是得 从而有 （7.3）3．,且。当k充分大时，则有 在实际应用幂法时，可根据迭代向量个分量的变化情况判断属于那种情况。 若迭代向量各分量单调变化，且有关系式，则属于第1种情况； 若迭代向量各分量不是单调变化，但有关系式，则属于第2种情况； 一、算法设计方案 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 ξ???? ? ?=?????? ， 北京航空航天大学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)北航数值分析大作业一
幂法及反幂法
北航数值分析大作业第一题幂法与反幂法
实验8 反幂法
数值分析之幂法及反幂法C语言程序实例
北航数值分析报告第三次大作业
幂法和反幂法的matlab实现
数值分析幂法c语言实现
幂法反幂法求解矩阵大小特征值及其对应的特征向量
数值分析幂法与反幂法-matlab程序
北航数值分析大作业第二题精解
北航数值分析大作业第二题
北航数值分析报告大作业第八题
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序
数值分析试验幂法与反幂法matlab
北航数值分析课程第一次大作业讲解
乘幂反幂法
北航数值分析大作业3
BUAA数值分析大作业三
- 北航数值分析大作业一
- 北航数值分析大作业-第三题
- 北航数值分析报告第三次大作业
- 北航数值分析大作业3
- 北航数值分析大作业第二题(fortran)
- 北航数值分析大作业一(纯原创,高分版)
- 北航数值分析大作业 第一题 幂法与反幂法
- 北航数值分析大作业第二题
- 北航数值分析大作业第二次
- 北航数值分析大作业二纯高分版
- 北航数值分析报告大作业第三题(fortran)
- 2015年北航数值分析考试重点NEW
- 北航数值分析大作业第一题幂法与反幂法
- 北航数值分析课程第一次大作业讲解
- 北航数值分析大作业第二题
- 北航数值分析编程大作业一
- 北航数值分析大作业二
- 北航数值分析大作业第二题精解
- BUAA数值分析大作业三
- 北航数值分析第二次大作业--QR分解