反幂法Matlab程序设计
(Fortran编程)数值分析之反幂法
integer::N,i,k=1
real*8::P(N,N),L(N,N),U(N,N)
real*8::uk(N),u0(N),v(N),zc(N),y(N)
do i=1,N
u0(i)=1d0
!以上交换两行元素,然后开始消元
do I=K+1,N
TEMP=AB(I,K)/AB(K,K)
AB(I,:)=AB(I,:)-TEMP*AB(K,:)
end do
end do
AUP(:,:)=AB(1:N,1:N)
call MAX_ROU(v,N,m1)
uk=v/m1
write(21,98) k,v,uk,m1
98 format(I5,T15,3F16.10,T70,3F16.10,F16.10/)
ne SOLVE
real*8::AUP(N,N),BUP(N)
!A,B为增广矩阵
real*8::AB(N,N+1)
real*8::VTEMP1(N+1),VTEMP2(N+1)
AB(1:N,1:N)=A
AB(:,N+1)=B
!列主元消去法的核心部分
do K=1,N-1
end do
X(I)=X(I)/A(I,I)
end do
end subroutine UP_TRI
end module M_GAUSS
module POWER
use M_GAUSS
contains
subroutine SOLVE(P,L,U,N)
BUP(:)=AB(:,N+1)
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量-附Matlab程序
矩阵的特征值与特征向量的计算摘要物理,力学,工程技术中的很多问题在数学上都归结于求矩阵特征值的问题,例如振动问题(桥梁的振动,机械的振动,电磁振动等)、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分,本文使用幂法和反幂法分别求矩阵的按模最大,按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法,它最大的优点是方法简单,对于稀疏矩阵比较合适,但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值,及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词:矩阵;特征值;特征向量;冥法;反冥法THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACTPhysics, 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 theoretical physics in some of the problems. Matrix eigenvalue calculation is a very important part in matrix computation. In this paper, we use the power method and inverse power method to calculate the maximum of the matrix, according to the minimum characteristic vector and the corresponding characteristic value.Power method is an iterative method to calculate the eigenvalues of a matrix. It has the advantage that the method is simple and suitable for sparse matrices, but sometimes the convergence rate is very slow. The basic idea is to take a non - zero initial vector. Construct a vector sequence from the matrix of the matrix. Then the eigenvalues and eigenvectors are obtained by using the constructed vector sequence.The inverse power method is used to calculate the minimum feature vectors and their eigenvalues of the matrix, and to calculate the eigenvalues of the matrix. In this paper, we use the inverse power method to calculate the minimum eigenvalue of a matrix and its corresponding eigenvalues. The basic idea of calculating the minimum characteristic vector of a matrix is to transform it to the maximum characteristic vector of the modulus of the inverse matrix. Then, according to the model, the minimum feature vector of the original matrix is introduced.Key words: Matrix;Eigenvalue;Eigenvector;Iteration methods;目录1 引言 (1)2 相关定理。
幂法和反幂法的matlab实现
幂法和反幂法的matlab实现幂法求矩阵主特征值及对应特征向量摘要矩阵特征值的数值算法,在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题,所以说研究利用数学软件解决求特征值的问题是非常必要的。
实际问题中,有时需要的并不是所有的特征根,而是最大最小的实特征根。
称模最大的特征根为主特征值。
幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法,它最大的优点是方法简单,特别适用于大型稀疏矩阵,但有时收敛速度很慢。
用java来编写算法。
这个程序主要分成了四个大部分:第一部分为将矩阵转化为线性方程组;第二部分为求特征向量的极大值;第三部分为求幂法函数块;第四部分为页面设计及事件处理。
其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法,再经过一系列的验证和迭代得到结果。
关键字:主特征值;特征向量;线性方程组;幂法函数块POWER METHOD FOR FINDING THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF THEMATRIXABSTRACTNumerical algorithm for the eigenvalue of matrix, in science and engineering technology, alot 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 isthe exponentiation function block. The fourth part is the page design and eventprocessing .The basic process is a power law function block by calling the matrix is transformed into linear equations method, after a series of validation and iteration results.Power method for finding the eigenvalues and corresponding eigenvectors of the matrixKey words: Main eigenvalue; characteristic vector; linear equations; power function block、目录1幂法......................................................... . (1)1.1幂法的基本理论和推导 (1)1.2幂法算法的迭代向量规范化 (2)2概要设计........................................................ (3)2.1设计背景 (3)2.2运行流程........................................... . (3)2.3运行环境........................................... (3)3程序详细设计 (4)3.1矩阵转化为线性方程组……..………………………………………. .43.2特征向量的极大值 (4)3.3求幂法函数块............….....…………...…......…………………………3.4界面设计与事件处理..........................................................................4运行过程及结果................................................ (6)4.1 运行过程....................................... ..................………………………………………. .64.2 运行结果................................................ .. (6)4.3 结果分析.......................................... (6)5结论 (7)参考文献 (8)附录 (56)1 幂法设实矩阵nn ijaA ⨯=)(有一个完备的特征向量组,其特征值为n λλλ ,,21,相应的特征向量为nx x x ,,21。
matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量
竭诚为您提供优质文档/双击可除matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量篇一:幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1.幂法简介:当矩阵a满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵a需要满足的条件为:(1)|1||2|...|n|0,i为a的特征值xn(2)存在n个线性无关的特征向量,设为x1,x2,...,1.1计算过程:n对任意向量x,有x(0)(0)iui,i不全为0,则有i1x(k1)ax(k)...ak1x(0)aαiuiαiλik1uik1i1i1nnnk12k1λ1u1()a2u2()anun11k111u1k112|越小时,收敛越快;且当k充分大时,有可见,当|1 (k1)k111u1x(k1)x(k1)(k)x1(k),对应的特征向量即是。
kxx11u12算法实现(1).输入矩阵a,初始向量x,误差限,最大迭代次数n(2).k1,0;y(k)x(k)max(abs(x(k))(3).计算xay,max(x);(4).若||,输出,y,否则,转(5)(5).若kn,置kk1,,转3,否则输出失败信息,停机.3matlab程序代码function[t,y]=lpowera,x0,eps,n)%t为所求特征值,y 是对应特征向量k=1;z=0;%z相当于y=x0./max(abs(x0));%规范化初始向量x=a*y;%迭代格式b=max(x);%b相当于ifabs(z-b) t=max(x);return;endwhileabs(z-b)>epsz=b;y=x./max(abs(x));x=a*y;b=max(x);end[m,index]=max(a(matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量)bs(x));%这两步保证取出来的按模最大特征值t=x(index);%是原值,而非其绝对值。
幂法与反幂法
幂法与反幂法1 功能 幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。
反幂法用来计算矩阵按模最小的特征值及其特征向量,也可用来计算对应与一个给定近似特征值的特征向量。
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.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).3 Matlab 程序的实现3.1 幂法function [m,u]=pow(A,ep,N)%A 为矩阵;ep 为精度要求;N 为最大迭代次数;m 为绝对值最大的特征值;u 为对应最大特征值的特征向量。
N=100;ep=1e-6;n=length(A);u=ones(n,1);index=0;k=0;m1=0;while k<=Nv=A*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;if abs(m-m1)<epindex=1;break;endm1=m;k=k+1;end输入:A=[7 3 -2;3 4 -1;-2 -1 3];[m,u]=pow(A,1e-6) Enter结果:m = 9.6056u =1.00000.6056-0.39444.2 反幂法function[m ,u]=pow_inv(A,ep,N)%A为矩阵;ep为精度要求;N为最大迭代次数;m为绝对值最大的特征值;u为对应最大特征值的特征向量。
matlab反解方程
matlab反解方程标题:使用MATLAB反解方程的实用方法引言:在科学研究和工程应用中,解决非线性方程组是一个常见的问题。
尽管有许多数值方法可以用于求解这些方程,但MATLAB作为一种强大的计算工具,提供了一种更加直观和便捷的方式来反解方程。
本文将介绍如何使用MATLAB来反解方程,并提供一些实用的方法和技巧。
一、MATLAB基础知识在使用MATLAB反解方程之前,我们需要了解一些基础知识。
MATLAB是一种高级编程语言和数值计算环境,它具有强大的数值计算和图形处理能力。
我们可以使用MATLAB来进行数值计算、矩阵运算、符号计算、绘图等操作。
二、使用MATLAB反解方程的步骤1. 定义方程:首先,我们需要将要反解的方程用MATLAB的语法进行定义。
例如,我们要反解方程 f(x) = 0,可以使用MATLAB的符号计算工具箱来定义这个方程。
2. 求解方程:接下来,我们可以使用MATLAB提供的数值计算工具箱来求解方程。
MATLAB提供了许多求解非线性方程的函数,例如fzero、fsolve等。
这些函数可以帮助我们找到方程的根。
3. 分析结果:在求解方程之后,我们可以使用MATLAB的绘图和数据分析工具来分析结果。
例如,我们可以绘制方程的图像,或者计算方程的根的性质。
三、MATLAB反解方程的实例为了更好地理解如何使用MATLAB反解方程,我们将通过一个实例来演示。
假设我们要反解方程x^2 + 2x + 1 = 0。
我们可以使用MATLAB的符号计算工具箱来定义这个方程:syms xeqn = x^2 + 2*x + 1;然后,我们可以使用MATLAB的数值计算工具箱来求解方程:x = solve(eqn, x);我们可以绘制方程的图像来分析结果:fplot(eqn, [-5, 5]);通过MATLAB的求解和绘图功能,我们可以得到方程的根为x = -1。
同时,我们还可以观察到方程的图像是一个开口向上的抛物线。
幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介:当矩阵A 满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵A 需要满足的条件为:(1) 存在n 个线性无关的特征向量,设为n x x x ,...,,211.1计算过程:i n i i i u xx αα,1)0()0(∑==,有对任意向量不全为0,则有 可见,当||12λλ越小时,收敛越快;且当k 充分大时,有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k (k k )(k x x u x u x ,对应的特征向量即是)(k x 1+。
2 算法实现3 matlab 程序代码function [t,y]=lpowerA,*0,eps,N) % t 为所求特征值,y 是对应特征向量k=1;z=0; % z 相当于λy=*0./ma*(abs(*0)); % 规化初始向量*=A*y; % 迭代格式b=ma*(*); % b 相当于 βif abs(z-b)<eps % 判断第一次迭代后是否满足要求t=ma*(*);return ;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=*./ma*(abs(*));*=A*y;b=ma*(*);end[m,inde*]=ma*(abs(*)); % 这两步保证取出来的按模最大特征值t=*(inde*); % 是原值,而非其绝对值。
end4 举例验证选取一个矩阵A ,代入程序,得到结果,并与eig(A)的得到结果比拟,再计算A*y-t*y ,验证y 是否是对应的特征向量。
结果如下:结果正确,说明算法和代码正确,然后利用此程序计算15阶Hilb 矩阵,与eig(A)的得到结果比拟,再计算 A*y-t*y ,验证y 是否是对应的特征向量。
设置初始向量为*0=ones(15,1),结果显示如下可见,结果正确。
幂法和反幂法的matlab实现
幂法求矩阵主特征值及对应特征向量摘要矩阵特征值的数值算法,在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题,所以说研究利用数学软件解决求特征值的问题是非常必要的。
实际问题中,有时需要的并不是所有的特征根,而是最大最小的实特征根。
称模最大的特征根为主特征值。
幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法,它最大的优点是方法简单,特别适用于大型稀疏矩阵,但有时收敛速度很慢。
用java来编写算法。
这个程序主要分成了四个大部分:第一部分为将矩阵转化为线性方程组;第二部分为求特征向量的极大值;第三部分为求幂法函数块;第四部分为页面设计及事件处理。
其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法,再经过一系列的验证和迭代得到结果。
关键字:主特征值;特征向量;线性方程组;幂法函数块POWER METHOD FOR FINDING THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF THEMATRIXABSTRACTNumerical 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 the exponentiation function block. The fourth part is the page design and event processing .The basic process is a power law function block by calling the matrix is transformed into linear equations method, after a series of validation and iteration results.Power method for finding the eigenvalues and corresponding eigenvectors of the matrixKey words: Main eigenvalue; characteristic vector; linear equations; power function block、目录1幂法 (1)1.1幂法的基本理论和推导 (1)1.2幂法算法的迭代向量规范化 (2)2概要设计 (3)2.1设计背景 (3)2.2运行流程 (3)2.3运行环境 (3)3程序详细设计 (4)3.1矩阵转化为线性方程组 (4)3.2特征向量的极大值 (4)3.3求幂法函数块............….....…………...…......…………………………3.4界面设计与事件处理............….....…………...…......…………………………4 运行过程及结果 (6)4.1 运行过程.........................................................………………………………………. .64.2 运行结果 (6)4.3 结果分析 (6)5结论 (7)参考文献 (8)附录 (56)1 幂法设实矩阵n n ij a A ⨯=)(有一个完备的特征向量组,其特征值为n λλλ ,,21,相应的特征向量为n x x x ,,21。
matlab幂法求特征值和特征向量方法实现和函数表示
matlab幂法求特征值和特征向量方法实现和函数表示1. 引言在数值分析中,求解特征值和特征向量是一项重要而且经常出现的任务。
特征值和特征向量在矩阵和线性代数中有着广泛的应用,涉及到许多领域,如机器学习、信号处理、结构动力学等。
在matlab中,幂法是一种常用的求解特征值和特征向量的方法,同时也有对应的函数可以实现这一过程。
2. 幂法的原理幂法是一种迭代方法,它利用矩阵的特征值和特征向量的性质,通过不断地迭代计算,逼近矩阵的主特征值和对应的特征向量。
具体来说,假设A是一个n阶矩阵,它的特征值λ1>λ2≥...≥λn,并且对应着线性无关的特征向量v1,v2,...,vn。
如果选择一个任意的非零初始向量x0,并进行以下迭代计算:```x(k+1) = Ax(k) / ||Ax(k)||```其中,||.||表示向量的模长。
不断迭代计算后,x(k)将收敛到矩阵A的主特征向量v1上,并且相应的特征值即为A的主特征值λ1。
3. matlab实现幂法求解特征值和特征向量在matlab中,幂法的实现也非常简单。
可以使用自带的eig函数,该函数可以直接求解矩阵的特征值和特征向量。
使用方法如下:```[V,D] = eig(A)```其中,A为待求解的矩阵,V为特征向量矩阵,D为特征值矩阵。
利用eig函数,即可一步到位地求解矩阵的特征值和特征向量,非常简单方便。
4. 函数表示幂法求解特征值和特征向量的过程可以表示为一个matlab函数。
通过封装相关的迭代算法和收敛判据,可以方便地实现幂法的函数表示。
可以定义一个名为powerMethod的函数:```matlabfunction [lambda, v] = powerMethod(A, x0, maxIter, tol)% 初始化k = 1;x = x0;% 迭代计算while k <= maxItery = A * x;lambda = norm(y, inf);x = y / lambda;% 检查收敛性if norm(A * x - lambda * x) < tolbreak;endk = k + 1;endv = x;end```利用这个函数,就可以自己实现幂法求解特征值和特征向量的过程。
数值分析课程设计+幂法与反幂法MATLAB
一、问题的描述及算法设计(一)问题的描述本次课程设计我所要做的课题是:对称矩阵的条件数的求解设计 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).二、算法的流程图(一)幂法算法的流程图(二)反幂法算法的流程图三、算法的理论依据及其推导(一)幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法,也即,绝对值最大的特征值。
稍微修改该方法,也可以用来确定其他特征值。
幂法的一个很有用的特性是它不仅可以生成特征值,而且可以生成相应的特征向量。
实际上,幂法经常用来求通过其他方法确定的特征值的特征向量。
数值分析3.1幂法和反幂法
第三章 矩阵的特征值与特征向量
3.1 幂法与反幂法 3.2 Jacobi方法
3.3 QR方法
第三章 矩阵的特征值与特征向量
3.1幂法与反幂法
一、乘幂法 二、反幂法
三、带原点位移的反幂法
四、反幂法的特点
第三章 矩阵的特征值与特征向量
3.1幂法与反幂法
一、乘幂法
1、基本思想
2、算法(迭代公式) ◆一般算法
具体算法: (1)使用范数 2
1 X 1 yk , k 1 1 X 1
(2)使用范数
uk A yk 1
令
k
er u k er y k 1
T
T
k
lim k 1
留为作业自学
具体算法: (1)使用范数 2 1 X 1 yk , k 1 1 X 1
1 2 n
第三章 矩阵的特征值与特征向量
一、乘幂法 1、基本思想 设A有n个线性无关的特征向量 X 1 , X 2 ,, X n ,
AX j j X j , j 1,2,, n
3.1幂法与反幂法
★ 设 1为实数而且是单根: 1 2 n
u0 1 X 1 2 X 2 n X n
具体算法: 按取范数的不同, 迭代公式也不同。 (1)使用范数 2
任取初始向量u0 R n T k 1 u k 1 u k 1 u k 1 yk 1 k 1 (3.4) u k A yk 1 k yk 1T uk k 1,2,
T
精确结果:
X 1 (0,0.5,1) , 1 45
T
max( uk ) 表示 u k 的绝对值最大的分量。 (3)
数值分析之幂法及反幂法C语言程序实例
u[N-1]=y[N-1]/A[S][N-1];
for(i=N-1;i>0;i--)
{
temp=0;
for(t=i+1;t<=min(i+S,N);t++)
temp+=A[i-t+S][t-1]*u[t-1];
u[i-1]=(y[i-1]-temp)/A[S][i-1];
}
}
double Det_matrix()//求矩阵行列式值
if(fabs((value2-value1)/value1)<e)
break;
}
value_s=value1;
}
float Get_cond_A()//求矩阵条件数
{
float cond1;
cond1=fabs(value_abs_max/value_s);
return cond1;
}
void Value_translation_min()//偏移条件下反幂法求特征值
{
int i;
double det=1;
Init_matrix_A();
Resolve_LU();
for(i=0;i<N;i++)
det=det*A[2][i];
return det;
}
float Get_norm()//获得迭代向量模
{
int i;
float normal=0;
for(i=0;i<N;i++)
Init_u();
norm=Get_norm();
Get_yy(norm);
Back_substitution();
反幂法matlab程序
反幂法matlab程序
反幂法是一种求解矩阵最小特征值的方法,其原理是通过对原矩阵进行迭代求解,逐步逼近最小特征值的值和对应的特征向量。
以下是一个简单的反幂法的 MATLAB 程序:
% 反幂法求解矩阵最小特征值及对应特征向量
function [lambda, x] = inverse_power_iter(A, x0, max_iter, tol)
% 输入:A为矩阵,x0为初始向量,max_iter为最大迭代次数,tol为误差容限
% 输出:lambda为最小特征值,x为对应特征向量
n = size(A,1); % 矩阵维度
x = x0/norm(x0); % 归一化初始向量
for k = 1:max_iter
y = A\x; % 解Ax = y,速度较inv(A)*x快
lambda = x'*y; % 计算特征值
x = y/norm(y); % 归一化特征向量
if norm(A*x-lambda*x)<tol % 满足精度要求,停止迭代
return
end
end
disp('达到最大迭代次数,未满足精度要求!')。
matlab多项式负幂次拟合
matlab多项式负幂次拟合在MATLAB中进行多项式负幂次拟合,通常可以使用polyfit函数来实现。
这个函数可以用来拟合一个多项式模型,其形式为y = p1x^n + p2x^(n-1) + ... + pnx + c,其中n为多项式的次数。
如果需要进行负幂次拟合,可以通过对x取倒数来实现。
下面我将从几个方面来介绍如何在MATLAB中进行多项式负幂次拟合。
首先,你需要准备好你的数据集,包括自变量x和因变量y。
然后,你可以使用1./x来对x取倒数,得到新的自变量。
接下来,可以使用polyfit函数来进行负幂次拟合。
例如,如果你想拟合一个负二次多项式模型,可以使用以下代码:matlab.new_x = 1./x;p = polyfit(new_x, y, 2);这将会拟合一个负二次多项式模型,并返回多项式系数p。
另外,你也可以使用Curve Fitting Toolbox中的fit函数来进行拟合。
你可以选择Custom Equation类型,并输入你想要拟合的负幂次多项式模型。
然后,使用fit函数拟合你的数据,并得到拟合结果。
在进行多项式负幂次拟合时,需要注意拟合结果的解释和评估。
你可以使用拟合后的多项式模型来预测新的数据点,并通过观察拟合曲线与原始数据的拟合程度来评估拟合质量。
此外,还可以计算拟合后的残差来评估拟合的准确性。
总之,在MATLAB中进行多项式负幂次拟合,你可以使用polyfit函数或Curve Fitting Toolbox中的fit函数。
拟合结果的准确性需要通过对拟合曲线与原始数据的比较和残差分析来评估。
希望这些信息能够对你有所帮助。
幂法及其MATLAB程序
5.2 幂法及其MATLAB 程序5.2.2 幂法的MATLAB 程序用幂法计算矩阵A 的主特征值和对应的特征向量的MATLAB 主程序function [k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0;while ((k<=max1)&(state==1))Vk=A*V; [m j]=max(abs(Vk)); mk=m;tzw=abs(lambda-mk); Vk=(1/mk)*Vk;Txw=norm(V-Vk); Wc=max(Txw,tzw); V=Vk;lambda=mk;state=0;if (Wc>jd)state=1;endk=k+1;Wc=Wc;endif (Wc<=jd)disp('请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('请注意:迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:') endVk=V;k=k-1;Wc;例 5.2.2 用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量,精度510-=ε.并把(1)和(2)输出的结果与例5.1.1中的结果进行比较.(1)⎪⎪⎭⎫ ⎝⎛-=4211A ; (2)⎪⎪⎪⎭⎫ ⎝⎛=633312321B ;(3)⎪⎪⎪⎭⎫ ⎝⎛--=1124111221C ;(4)⎪⎪⎪⎭⎫ ⎝⎛---=20101350144D . 解 (1)输入MATLAB 程序>>A=[1 -1;2 4]; V0=[1,1]';[k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100),[V,D] = eig (A), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =33 3.00000173836804 8.691862856124999e-007Vk = V = wuV =-0.49999942054432 -0.70710678118655 0.44721359549996 -0.894428227562941.00000000000000 0.70710678118655 -0.89442719099992 -0.89442719099992Dzd = wuD =3 1.738368038406435e-006由输出结果可看出,迭代33次,相邻两次迭代的误差W c ≈8.69 19e-007,矩阵A 的主特征值的近似值lambda ≈3.000 00和对应的特征向量的近似向量V k ≈(-0.500 00,1.00000T ), lambda 与例5.1.1中A 的最大特征值32=λ近似相等,绝对误差约为1.738 37e-006,V k 与特征向量X =T22k T )1,21(- )0(2≠k 的第1个分量的绝对误差约等于0,第2个分量的绝对值相同.由wuV 可以看出,2λ的特征向量V (:,2) 与V k 的对应分量的比值近似相等.因此,用程序mifa.m 计算的结果达到预先给定的精度510-=ε.(2) 输入MATLAB 程序>>B=[1 2 3;2 1 3;3 3 6]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(B,V0,0.00001,100), [V,D] = eig (B), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = Dzd = wuD =3 9 0 9 0Vk = wuV =0.50000000000000 0.816496580927730.50000000000000 0.816496580927731.00000000000000 0.81649658092773V =0.70710678118655 0.57735026918963 0.40824829046386-0.70710678118655 0.57735026918963 0.408248290463860 -0.57735026918963 0.81649658092773(3) 输入MATLAB 程序>> C=[1 2 2;1 -1 1;4 -12 1];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(C,V0,0.00001,100), [V,D] = eig (C), Dzd=max(diag(D)), wuD=abs(Dzd-lambda),Vzd=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示请注意:迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =100 0.09090909090910 2.37758124193119Dzd = wuD =1.00000000000001 0.90909090909091Vk= Vzd = wuV =0.99999999999993 0.90453403373329 0.904534033733350.99999999999995 0.30151134457776 0.301511344577781.00000000000000 -0.30151134457776 -0.30151134457776由输出结果可见,迭代次数k 已经达到最大迭代次数max 1=100,并且lambda 的相邻两次迭代的误差Wc ≈2.377 58>2,由wuV 可以看出,lambda 的特征向量V k 与真值Dzd 的特征向量V zd 对应分量的比值相差较大,所以迭代序列发散.实际上,实数矩阵C 的特征值的近似值为i ,i ,010*********.000321=-==λλλ ,并且对应的特征向量的近似向量分别为X T1=1k (0.90453403373329,0.30151134457776,-0.30151134457776)T ,X =T 22k (-0.72547625011001,-0.21764287503300-0.07254762501100i, 0.58038100008801-0.29019050004400i )T ,X =T33k ( -0.72547625011001, -0.21764287503300 + 0.07254762501100i,0.58038100008801 + 0.29019050004400i)T0,0(21≠≠k k , 03≠k 是常数).(4)输入MATLAB 程序>> D=[-4 14 0;-5 13 0;-1 0 2]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(D,V0,0.00001,100), [V,Dt] =eig (D), Dtzd=max(diag(Dt)), wuDt=abs(Dtzd-lambda),Vzd=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下:k = lambda = Wc =19 6.00000653949528 6.539523793591684e-006Dtzd = wuDt =6.00000000000000 6.539495284840768e-006Vk = Vzd = wuV =0.79740048053564 0.79740048053564 0.797400480535640.71428594783886 0.56957177181117 0.79740021980618-0.24999918247180 -0.19935012013391 0.797403088133705.3 反幂法和位移反幂法及其MATLAB程序5.3.3 原点位移反幂法的MATLAB程序(一)原点位移反幂法的MATLAB主程序1用原点位移反幂法计算矩阵A的特征值和对应的特征向量的MATLAB主程序1 function [k,lambdan,Vk,Wc]=ydwyfmf(A,V0,jlamb,jd,max1)[n,n]=size(A); A1=A-jlamb*eye(n); jd= jd*0.1;RA1=det(A1);if RA1==0disp('请注意:因为A-aE的n阶行列式hl等于零,所以A-aE不能进行LU分解.')returnendlambda=0;if RA1~=0for p=1:nh(p)=det(A1(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意:因为A-aE的r阶主子式等于零,所以A-aE不能进行LU分解.')returnendendif h(1,i)~=0disp('请注意:因为A-aE的各阶主子式都不等于零,所以A-aE 能进行LU分解.')k=1;Wc =1;state=1; Vk=V0;while((k<=max1)&(state==1))[L U]=lu(A1); Yk=L\Vk;Vk=U\Yk; [mj]=max(abs(Vk));mk=m;Vk1=Vk/mk; Yk1=L\Vk1;Vk1=U\Yk1;[m j]=max(abs(Vk1));mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1);tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1);Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1);Txw=min(Txw1,Txw2); tzw=min(tzw1,tzw2);Vk=Vk2;mk=mk1; Wc=max(Txw,tzw);Vk=Vk2;mk=mk1;state=0;if(Wc>jd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif (Wc<=jd)disp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k 已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:')endhl,RA1endend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;例5.3.2 用原点位移反幂法的迭代公式(5.28),根据给定的下列矩阵的特征值n λ的初始值n λ~,计算与n λ对应的特征向量n X 的近似向量,精确到0.000 1. (1)⎪⎪⎪⎭⎫ ⎝⎛----210242011,2.0~2=λ;(2)⎪⎪⎭⎫ ⎝⎛-4211,001.2~2=λ;(3)⎪⎪⎪⎭⎫ ⎝⎛--3315358215211,8.26~3=λ.解 (1)输入MATLAB 程序>> A=[1 -1 0;-2 4 -2;0 -1 2];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.2,0.0001,10000)运行后屏幕显示结果 请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = hl =3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D =1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 00.7616 1.0000 -0.7616 0.3633 0 0.2384 00.4323 -0.3200 -0.4323 1.0000 0 0 1.6367(2)输入MATLAB 程序>> A=[1 -1;2 4];V0=[20,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.0001,100)运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = hl =2 2.0020 5.1528e-007 -1.0010 -0.0010Vk = V = D =1.0000 -1.0000 0.5000 2 0-1.0000 1.0000 -1.0000 0 3(3)输入MATLAB 程序>> A=[-11 2 15;2 58 3;15 3 -3];V0=[1,1,-1]';[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,8.26, 0.0001,100)运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambdan= Wc = hl =2 8.2640 6.9304e-008 -19.2600 -961.9924 -6.1256Vk = V = D =-0.7692 0.7928 0.6081 0.0416 -22.5249 0 00.0912 0.0030 -0.0721 0.9974 0 8.2640 0-1.0000 -0.6095 0.7906 0.0590 0 0 58.2609例 5.3.3 用原点位移反幂法的迭代公式(5.28),计算⎪⎪⎪⎭⎫ ⎝⎛-----=1026471725110A 的分别对应于特征值 1.001~11=≈λλ,.001 2~22=≈λλ, 001.4~33=≈λλ的特征向量1X ,2X ,3X 的近似向量,相邻迭代误差为0.001.将计算结果与精确特征向量比较. 解 (1)计算特征值 1.001~11=≈λλ对应的特征向量1X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]= ydwyfmf(A,V0,1.001, 0.001,100),[V,D]=eig(A);Dzd=min(diag(D)), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-1.00100000000000 5.98500100000000 -0.00299600100000k = lambda = RA1 =5 1.00200000000000 -0.00299600100000Vk = VD = wuV =-0.50000000000000 -0.40824829046386 0.81649658092773-0.50000000000000 -0.40824829046386 0.81649658092773-1.00000000000000 -0.81649658092773 0.81649658092773Wc = Dzd = wuD =1.378794763695562e-009 1.00000000000000 0.00200000000000 从输出的结果可见,迭代5次,特征向量1X 的近似向量1~X 的相邻两次迭代的误差Wc ≈1.379 e-009,由wuV 可以看出,1~X = Vk 与VD 的对应分量的比值相等.特征值1λ的近似值lambda ≈1.002与初始值=1~λ 1.001的绝对误差为0.001,而与 1λ的绝对误差为0.002,其中 =1X T )000000000001.000 , 000000000000.500- , 000000000000.500( -, =1~X T )000000000001.000 , 000000000000.500- , 000000000000.500(-. (2)计算特征值.001 2~22=≈λλ对应特征向量2X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001, 0.001,100) ,[V,D]=eig(A); WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-2.00100000000000 -8.01299900000000 0.00200099900000k = Wc = lambda = WD =2 3.131363162302120e-007 2.00200000000016 0.00200000000016Vk = VD = wuV =-0.24999999999999 0.21821789023599 -0.87287156094401 -0.49999999999999 0.43643578047198 -0.87287156094398 -1.00000000000000 0.87287156094397 -0.87287156094397 从输出的结果可见,迭代2次,特征向量2X 的近似向量2~X 的相邻两次迭代的误差Wc ≈3.131e-007,2~X 与2X 的对应分量的比值近似相等.特征值2λ的近似值lambda ≈2.002与初始值=2~λ 2.001的绝对误差约为0.001,而lambda 与2λ的绝对误差约为0.002,其中 =2~X T )00000000000000.1,99999999999499.0,99999999999249.0(---, =2X T ) 000000000001.000- ,000000000000.500- ,99999999999-0.249( . (3)计算特征值 001.4~33=≈λλ对应特征向量3X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,4.001, 0.001,100)[V,D]=eig(A);WD=lambda-max(diag(D)),VD=V(:,3),wuV=V(:,3)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-4.00100000000000 -30.00899900000000 -0.00600500099999 k = lambda = Wc = WD =2 4.00199999999990 1.996084182914842e-007 0.00199999999990Vk = VD = wuV =0.40000000000001 -0.32444284226153 -0.81110710565380 0.60000000000001 -0.48666426339229 -0.81110710565381 1.00000000000000 -0.81110710565381 -0.81110710565381 从输出的结果可见,迭代2次,特征向量3X 的近似向量3~X 的相邻两次迭代的误差Wc ≈1.996e-007,3~X 与3X 的对应分量的比值近似相等.特征值3λ的近似值 4.001~4.0022=≈λ与初始值lambda 的绝对误差近似为001.0,而lambda 与3λ的绝对误差约为0.002,其中 =3X (-0.400 000 000 000 00,-0.600 000 000 000 00,-1.000 000 000 000 00T ), =3~X T )000000000001.000 ,100000000000.600 ,10000000000.400(.(二)原点位移反幂法的MATLAB 主程序2用原点位移反幂法计算矩阵A 的特征值和对应的特征向量的MATLAB 主程序2function [k,lambdan,Vk,Wc]=wfmifa1(A,V0,jlamb,jd,max1)[n,n]=size(A); jd= jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1); lambda1=0;k=1;Wc =1;state=1; U=V0;while ((k<=max1)&(state==1))Vk=A1\U; [m j]=max(abs(Vk)); mk=m; Vk=(1/mk)*Vk;Vk1=A1\Vk;[m1 j]=max(abs(Vk1)); mk1=m1,Vk1=(1/mk1)*Vk1;U=Vk1,Txw=(norm(Vk1)-norm(Vk))/norm(Vk1);tzw=abs((lambda1-mk1)/mk1);Wc=max(Txw,tzw); lambda1=mk1;state=0;if (Wc>jd)state=1;endk=k+1;endif (Wc<=jd)disp('请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('请注意迭代次数k 已经达到最大迭代次数max1, 特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:') end[V,D] =eig(A,'nobalance'),Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例5.3.4 用原点位移反幂法的迭代公式(5.27),计算例题5.3.3,并且将这两个例题的计算结果进行比较.再用两种原点位移反幂法的MATLAB 主程序,求979999999990.999~1=λ对应的特征向量. 解 (1)计算特征值 1.001~11=≈λλ对应特征向量1X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,1.001,0.001,100)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =5 1.00200000000138 1.376344154436924e-006Vk’ = -0.50000000000000 -0.50000000000000 -1.00000000000000同理可得,另外与两个特征值对应的特征向量.(2)再用两种原点位移反幂法的MATLAB 主程序,求979999999990.999~1=λ对应的特征向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.99999999999997,0.001,100) 运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-0.99999999999997 6.00000000000045 0.00000000000010RA1 = 1.039168751049192e-013 k = 2 lambda = 1.00000000000000输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0, 0.99999999999997,0.001,100) 运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = 3 lambda = 1.00000000000000 Wc =5.412337245047640e-016Vk = 0.50000000000000 0.50000000000000 1.00000000000000 Wc = 4.317692037236759e-013 Vk =0.500000000000000.500000000000001.000000000000005.4 雅可比(Jacobi)方法及其MATLAB 程序5.4.3 雅可比方法的MATLAB 程序用雅可比方法计算对称矩阵A 的特征值和对应的特征向量的MATLAB 主程序function [k,Bk,V,D,Wc]=jacobite(A,jd,max1)[n,n]=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n); Aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(Aij);[m2 j]=max(m1);i=i(j);while ((k<=max1)&(state==1))k=k+1,aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(abs(aij));[m2 j]=max(m1);i=i(j),j,Aij=(Bk-diag(diag(Bk)));mk=m2*sign(Aij(i,j)),Wc=m2,Dk=diag(diag(Bk));Pk=P0;c=(Bk(j,j)-Bk(i,i))/(2*Bk(i,j)),t=sign(c)/(abs(c)+sqrt(1+c^2)),pii=1/( sqrt(1+t^2)), pij=t/( sqrt(1+t^2)),Pk(i,i)=pii;Pk(i,j)=pij;Pk(j,j)=pii; Pk(j,i)=-pij;Pk,B1=Pk'*Bk;B2=B1*Pk; Vk=Vk*Pk,Bk=B2,if (Wc>jd)state=1;elsereturnendPk;Vk;Bk=B2;Wc;endif (k>max1)disp('请注意迭代次数k 已经达到最大迭代次数max1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:')elsedisp('请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:')endWc;k=k; V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1]=eig(A,'nobalance')例5.4.2 用雅可比方法的MATLAB 程序计算矩阵⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=12101152302756135612A 的特征值i λ和对应的特征向量i X (4,3,2,1=i ).解 (1)保存名为jacobite.m 为M 文件;(2)输入MATLAB 程序>> A=[12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12];[k,B,V,D,Wc]=jacobite(A,0.001,100)(3)运行后屏幕显示如下:k = i = j = mk = Wc =1 2 1 -56 56c = t =-0.04464285714286 -0.95635313919972pii = pij =0.72270271801843 -0.69115901308510Pk =0.72270271801843 0.69115901308510 0 0 -0.69115901308510 0.72270271801843 0 0 0 0 1.00000000000000 00 0 0 1.00000000000000Vk =0.72270271801843 0.69115901308510 0 0 -0.69115901308510 0.72270271801843 0 00 0 1.00000000000000 00 0 0 1.00000000000000Bk =65.55577579518456 0 0.78579012788509 -0.72270271801843 -0.00000000000001 -46.55577579518456 3.51888247529217 -0.691159013085100.78579012788509 3.51888247529217 5.00000000000000 1.00000000000000 -0.72270271801843 -0.69115901308510 1.00000000000000 12.00000000000000 k =i = j = mk = Wc =2 3 2 3.51888247529217 3.51888247529217c = t =-7.32558932518824 -0.06793885568129pii = pij =0.99770011455446 -0.06778260409592Pk =1.00000000000000 0 0 00 0.99770011455446 0.06778260409592 00 -0.06778260409592 0.99770011455446 00 0 0 1.00000000000000Vk =0.72270271801843 0.68956942653035 0.04684855775127 0 -0.69115901308510 0.72104058455581 0.04898667221449 00 -0.06778260409592 0.99770011455446 00 0 0 1.00000000000000Bk =65.55577579518456 -0.05326290114092 0.78398290060672 -0.72270271801843 -0.05326290114093 -46.79484464383285 0 -0.757352030626270.78398290060672 0.00000000000000 5.23906884864829 0.95085155680318 -0.72270271801843 -0.75735203062627 0.95085155680318 12.00000000000000 k = i = j = mk = Wc =3 4 3 0.95085155680318 0.95085155680318c = t =-3.55519802380213 -0.13796227443116pii = pij =0.99061693994324 -0.13666776612460Pk =1.00000000000000 0 0 00 1.00000000000000 0 00 0 0.99061693994324 0.136667766124600 0 -0.13666776612460 0.99061693994324 Vk =0.72270271801843 0.68956942653035 0.04640897492032 0.00640268773403 -0.69115901308510 0.72104058455581 0.04852702732712 0.006694899061430 -0.06778260409592 0.98833863446096 0.136353445918420 0 -0.13666776612460 0.99061693994324 Bk =65.55577579518456 -0.05326290114092 0.87539690801061 -0.60877636330628 -0.05326290114093 -46.79484464383285 0.10350561019562 -0.750245751038800.87539690801061 0.10350561019562 5.10788720522532 -0.00000000000000 -0.60877636330628 -0.75024575103880 -0.00000000000000 12.13118164342297 k =i = j = mk = Wc =4 1 3 0.87539690801061 0.87539690801061c = t =-34.52598931799430 -0.01447880833914pii = pij =0.99989519853186 -0.01447729093877Pk =0.99989519853186 0 -0.01447729093877 00 1.00000000000000 0 00.01447729093877 0 0.99989519853186 00 0 0 1.00000000000000Vk =0.72329885394465 0.68956942653035 0.03594133368062 0.00640268773403 -0.69038403871280 0.72104058455581 0.05852805174080 0.006694899061430.01430846595712 -0.06778260409592 0.98823505512105 0.13635344591842-0.00197857901214 0 -0.13665344314206 0.99061693994324Bk =65.56845049923633 -0.05175883827808 -0.00000000000000 -0.60871256264964-0.05175883827809 -46.79484464383285 0.10426586517177 -0.75024575103880-0.00000000000000 0.10426586517177 5.09521250117356 0.00881343252823-0.60871256264964 -0.75024575103880 0.00881343252823 12.13118164342297 k = i = j = mk = Wc =5 4 2 -0.75024575103880 0.75024575103880c = t =39.27114962375084 0.01272992971264pii = pij =0.99991898429114 0.01272889838836Pk =1.00000000000000 0 0 00 0.99991898429114 0 -0.012728898388360 0 1.00000000000000 00 0.01272889838836 0 0.99991898429114Vk =0.72329885394465 0.68959505973603 0.03594133368062 -0.00237529014628-0.69038403871280 0.72106738763160 0.05852805174080 -0.002483695665250.01430846595712 -0.06604148348220 0.98823505512105 0.13720519702737-0.00197857901214 0.01260946237032 -0.13665344314206 0.99053668440964Bk =65.56845049923633 -0.05950288535679 -0.00000000000000 -0.60800441437674-0.05950288535680 -46.80439521951078 0.10436960328590 0.00000000000000-0.00000000000000 0.10436960328590 5.09521250117356 0.00748552889860-0.60800441437674 0.00000000000000 0.00748552889860 12.14073221910090 k =i = j = mk = Wc =6 4 1 -0.60800441437674 0.60800441437674c = t =-43.93694931878409 -0.01137847012503pii = pij =0.99993527149402 -0.01137773361366Pk =0.99993527149402 0 0 0.011377733613660 1.00000000000000 0 00 0 1.00000000000000 0-0.01137773361366 0 0 0.99993527149402Vk =0.72327906130899 0.68959505973603 0.03594133368062 0.00585436528595-0.69031109235777 0.72106738763160 0.05852805174080 -0.010338540582940.01274645560931 -0.06604148348220 0.98823505512105 0.13735911385404-0.01324851347145 0.01260946237032 -0.13665344314206 0.99045005670500Bk =65.57536865930122 -0.05949903382392 -0.00008516835377 -0.00000000000000-0.05949903382393 -46.80439521951078 0.10436960328590 -0.00067700797883-0.00008516835377 0.10436960328590 5.09521250117356 0.00748504437150-0.00000000000000 -0.00067700797883 0.00748504437150 12.13381405903603 k =i = j = mk = Wc =7 3 2 0.10436960328590 0.10436960328590c = t =-2.486337309269764e+002 -0.00201098208240pii = pij =0.99999797798167 -0.00201097801616Pk =1.00000000000000 0 0 00 0.99999797798167 0.00201097801616 00 -0.00201097801616 0.99999797798167 00 0 0 1.00000000000000…………………………………………………………………………请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:V1 =0.68990429476497 -0.03732423222484 0.00588594854431 -0.722913771734500.72058252860300 -0.05998661236737 -0.01028322161977 0.69069289931337-0.06802029759277 -0.98795368410472 0.13841044442471 -0.012779125692250.01288885768193 0.13768088498200 0.99030407443219 0.01325486405899D1 =-46.80463661419736 0 0 00 5.09541442877727 0 00 0 12.13382202426702 00 0 0 65.57540016115307k =10B =65.57540016045945 0.00000000000175 -0.00020481967566 0.000000148628360.00000000000175 -46.80463661419739 0.00000062739984 0.00000000000000-0.00020481967566 0.00000062739984 5.09541442947090 -0.000000000007370.00000014862836 -0.00000000000000 -0.00000000000737 12.13382202426704V =0.72291389811507 0.68990429521617 0.03732177568689 0.00588595055487-0.69069269613201 0.72058252932816 0.05998894273570 -0.010283223540620.01278247108107 -0.06802028564977 0.98795364164379 0.13841044446122-0.01325533307898 0.01288885601755 -0.13768084024946 0.99030407439520D =65.57540016045945 0 0 00 -46.80463661419739 0 00 0 5.09541442947090 00 0 0 12.13382202426704Wc =6.920584967017158e-0045.5 豪斯霍尔德(Householder)方法及其MATLAB程序5.5.1 豪斯霍尔德方法及其MATLAB程序求初等反射矩阵P,使得PX的第一个分量以外的其余的分量都为零的MATLAB主程序function [xigema,rou,miou,P,PX]=Householder(X)n=size(X);nX=norm(X,2);xigema=nX*sign(X(1));rou=xigema*(xigema+X(1));miou=[xigema,zeros(1,n-1)]'+X,E=eye(n,n); C=2*miou*(miou)';P=E-C/(norm(miou,2)^2); PX=P*X;例5.5.1设向量=X()T1,2,2,确定一个初等反射矩阵P,使得PX的后两个分量为零.解输入MATLAB程序>> X=[2 2 1]'; [xigema,rou,miou,P,PX]=Householder(X)运行后屏幕显示结果P = PX =-0.6667 -0.6667 -0.3333 -3.0000-0.6667 0.7333 -0.1333 0.0000-0.3333 -0.1333 0.9333 0.00005.5.2 矩阵约化为上豪斯霍尔德矩阵及其MATLAB程序用豪斯霍尔德变换将n阶矩阵A规约成上豪斯霍尔德矩阵的MATLAB主程序function [k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k))*sign(Ak(k+1,k)),uk= Ak(k+1:n,k)+ Sk*eye(n-k,1),ck=(norm(uk,2)^2)/2,Pk= eye(n-k,n-k)-uk*uk'/ck,Uk=[eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk],A1=Uk*Ak;Ak=A1,end例5.5.3 用初等反射矩阵正交相似约化实矩阵A 为上豪斯霍尔德矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=34 19- 37 78- 41- 31 11 72- 98 10.2- 78- 32-94- 21 12 1 0 1- 63- 72 1 5 2 3 17- 32 02 7 56- 51- 17 12- 34 52- 12A . 解 输入MATLAB 程序>> A=[12 -52 34 -12 17 -51;-56 7 2 0 32 -17;3 2 5 1 72 -63;-1 0 1 12 21 -94;-32 -78 -10.2 98 -72 11;31 -41 -78 37 -19 34];[k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)运行后屏幕显示结果k = Sk = ck =1 -71.6310 9.1423e+003uk = Pk =-127.6310 -0.7818 0.0419 -0.0140 -0.4467 0.43283.0000 0.0419 0.9990 0.0003 0.0105 -0.0102-1.0000 -0.0140 0.0003 0.9999 -0.0035 0.0034-32.0000 -0.4467 0.0105 -0.0035 0.8880 0.108531.0000 0.4328 -0.0102 0.0034 0.1085 0.8949Uk =1.0000 0 0 0 0 00 -0.7818 0.0419 -0.0140 -0.4467 0.43280 0.0419 0.9990 0.0003 0.0105 -0.01020 -0.0140 0.0003 0.9999 -0.0035 0.00340 -0.4467 0.0105 -0.0035 0.8880 0.10850 0.4328 -0.0102 0.0034 0.1085 0.8949Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.76430.0000 1.8892 5.7655 1.6556 72.7134 -63.9112-0.0000 0.0369 0.7448 11.7815 20.7622 -93.6963-0.0000 -76.8184 -18.3655 91.0066 -79.6101 20.71910.0000 -42.1447 -70.0897 43.7749 -11.6277 24.5846k = Sk = ck =2 87.6402 7.8464e+003uk = Pk =89.5295 -0.0216 -0.0004 0.8765 0.48090.0369 -0.0004 1.0000 0.0004 0.0002-76.8184 0.8765 0.0004 0.2479 -0.4126-42.1447 0.4809 0.0002 -0.4126 0.7736Uk =1.0000 0 0 0 0 00 1.0000 0 0 0 00 0 -0.0216 -0.0004 0.8765 0.48090 0 -0.0004 1.0000 0.0004 0.00020 0 0.8765 0.0004 0.2479 -0.41260 0 0.4809 0.0002 -0.4126 0.7736Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.7643-0.0000 -87.6402 -49.9272 100.7790 -76.9476 31.4002-0.0000 -0.0000 0.7219 11.8223 20.7005 -93.6570-0.0000 0.0000 29.4202 5.9564 48.8026 -61.06030.0000 0.0000 -43.8731 -2.8860 58.8230 -20.2818…………………………………………………………………………k = Sk = ck =4 -12.2088 195.0398uk = Pk =-15.9753 -0.3085 0.951211.6133 0.9512 0.3085Uk =1.0000 0 0 0 0 00 1.0000 0 0 0 00 0 1.0000 0 0 00 0 0 1.0000 0 00 0 0 0 -0.3085 0.95120 0 0 0 0.9512 0.3085Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.7643-0.0000 -87.6402 -49.9272 100.7790 -76.9476 31.40020.0000 -0.0000 -52.8292 -5.8754 21.3902 18.44030.0000 0.0000 0.0000 12.2088 40.2435 -106.81340.0000 0.0000 -0.0000 0.0000 64.7555 -34.09095.5.3 实对称矩阵的三对角化及其MATLAB程序将n阶实对称矩阵A规约成三对角形式的MATLAB主程序function T=house(A)[n,n]=size(A);for k=1:n-2s=norm(A(k+1:n,k),2);if (A(k+1,k)<0)s=-s;endr=sqrt(2*s*(A(k+1,k)+s));U(1:k)=zeros(1,k);U(k+1)=(A(k+1,k)+s)/r;U(k+2:n)=A(k+2:n,k)'/r;V(1:k)=zeros(1,k);V(k+1:n)=A(k+1:n,k+1:n)*U(k+1:n)';C=U(k+1:n)*V(k+1:n)';P(1:k)=zeros(1,k);P(k+1:n)=V(k+1:n)-C*U(k+1:n);A(k+2:n,k)=zeros(n-k-1,1);A(k,k+2:n)=zeros(1,n-k-1);A(k+1,k)=-s; A(k,k+1)=-s;A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-2*U(k+1:n)'*P(k+1:n)-2*P( k+1:n)'*U(k+1:n);endT=A;例5.5.4 用初等反射矩阵正交相似约化实对称矩阵A为三对角矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------------------=5261215121416134237299021237312611451721253233219612371564901435612A 解 输入MATLAB 程序>> A=[12 -56 3 -14 -90 -4;-56 71 23 61 -9 -21;3 23 53 12 -72 51;-14 61 12 73 23 21;-90 -9 -72 23 -34 -61;-41 -21 51 21 -61 -52];T=house(A)运行后屏幕显示结果T =12.0000 114.5513 0 0 0 0114.5513 -43.2395 -108.2763 0 0 00 -108.2763 49.7411 -22.7766 0 00 0 -22.7766 40.2476 -89.1355 00 0 0 -89.1355 44.9606 39.30900 0 0 0 39.3090 19.29025.6 QR 方法及其MATLAB 程序5.6.5 最末元位移QR 法计算实对称矩阵特征值及其MATLAB 程序用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序function tzg=qr4(A,t,max1)[n,n]=size(A); k=0;Ak=A;tzg=zeros(n); state=1;for i=1:n;while ((k<=max1)&(state==1)&(n>1))b1=abs(Ak(n,n-1)); b2=abs(Ak(n,n));b3=abs(Ak(n-1,n-1));b4=min(b2, b3); jd=10^(-t); jd1=jd*b4;if (b1>=jd1)sk=Ak(n,n); Bk=Ak-sk*eye(n); [Qk,Rk]=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp('请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,')disp(' Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解,At=Rk*Qk+sk*eye(n)迭代矩阵:')i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;elsedisp('请注意:i 表示求第i 个特征值,tzgk 是矩阵A 的特征值的近似值,k 是迭代次数,')disp(' 下面的矩阵B 是m 阶矩阵At 的(m-1)阶主子矩阵,即At 降一阶.')i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state==1; i=i+1;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1));tzgk=Akdisp('请注意:n 阶实对称矩阵A 的全部真特征值lamoda 和至少含t个有效数字的近似特征值tzg 如下:')lamoda=sort(eig(A))例5.6.5 用最末元位移QR 方法求下列实对称矩阵的全部近似特征值,并将计算结果与A 全部真特征值比较.其中,2 1 1 1 1 3 1 21 1 4- 21 2 2 5)1(⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=A 精度为=ε510-; ,52612151214161342372990212373126114517212532332196123715641901435612)2(⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------------------=A 精度为=ε410-.解 (1)首先保存用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序为M 文件,取名为qr4.m.在MATLAB 工作窗口输入程序>> A=[5 2 2 1;2 -4 1 1;2 1 3 1;1 1 1 2]; tzg=qr4(A,5,100) 运行后屏幕显示结果请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =1k =1sk =2Bk =3 2 2 12 -6 1 12 1 1 11 1 1 0Qk =-0.70710678118655 0.38807526285317 0.12674485010490 -0.57735026918963-0.47140452079103 -0.87963726246718 0.06337242505245 0-0.47140452079103 0.20697347352169 -0.63372425052448 0.57735026918963-0.23570226039552 0.18110178933148 0.76046910062937 0.57735026918963 Rk =-4.24264068711929 0.70710678118655 -2.59272486435067 -1.649915822768610 6.44204936336256 0.28458852609232 -0.284588526092320 0 0.44360697536713 -0.443606975367130 0 0 0.00000000000000At =6.27777777777778 -3.10388935193069 -0.10455916682125 0.00000000000000-3.10388935193069 -3.65930388219545 0.01147685957127 0.00000000000000-0.10455916682125 0.01147685957127 1.38152610441767 0.00000000000000 -0.00000000000000 0.00000000000000 0.00000000000000 2.00000000000000 请注意:i 表示求第i 个特征值,tzgk 是矩阵A 的特征值的近似值,k 是迭代次数,下面的矩阵B 是m 阶矩阵At 的(m-1)阶主子矩阵,即At 降一阶.i =1tzgk =2.00000000000000k =1sk =2B =6.27777777777778 -3.10388935193069 -0.10455916682125-3.10388935193069 -3.65930388219545 0.01147685957127-0.10455916682125 0.01147685957127 1.38152610441767请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =2k =2sk =1.38152610441767Bk =4.89625167336011 -3.10388935193069 -0.10455916682125-3.10388935193069 -5.04082998661312 0.01147685957127-0.10455916682125 0.01147685957127 0Qk =-0.84445320114929 -0.53537837009187 0.016394874396770.53532568873289 -0.84460953959679 -0.007818734217300.01803324849744 0.00217404228940 0.99983502413586Rk =-5.79813264571247 -0.07718952005739 0.094439180886190 5.91931326753920 0.046285251232420 0 -0.00180396892170At =6.23815929000691 3.16959512520840 -0.000032531419853.16959512520840 -3.61788172311421 -0.00000392190472-0.00003253141985 -0.00000392190472 1.37972243310730请注意:i表示求第i个特征值,tzgk是矩阵A的特征值的近似值,k是迭代次数,下面的矩阵B是m阶矩阵At的(m-1)阶主子矩阵,即At降一阶.i =2tzgk =1.37972243310730k =2sk =1.38152610441767B =6.23815929000691 3.169595125208403.16959512520840 -3.61788172311421请注意:下面的i表示求第i个特征值,k是迭代次数,sk是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =3k =3sk =-3.61788172311421Bk =9.85604101312112 3.169595125208403.16959512520840 0Qk =-0.95198403663348 -0.30614766697629-0.30614766697629 0.95198403663348Rk =-10.35315786173815 -3.017403961789690 -0.97036415284199At =7.16193047323385 0.297074721510000.29707472151000 -4.54165290634115请注意:下面的i表示求第i个特征值,k是迭代次数,sk是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =3k =4sk =-4.54165290634115Bk =11.70358337957500 0.297074721510000.29707472151000 0。
matlab中exp的反函数
matlab中exp的反函数
在Matlab中,exp函数可以求出自然常数e的幂次方。
但是,如果我们想要求出e的幂次方的反函数,该怎么办呢?
答案是使用log函数。
在Matlab中,log函数可以求出以e为底数的对数。
因此,如果我们想要求出e的幂次方的反函数,可以使用log函数,具体代码如下:
y = log(x)
其中,x为自然常数e的幂次方的值,y为其反函数的值。
需要注意的是,由于log函数的底数为e,因此x的值必须为e的幂次方才能使用log函数求出其反函数。
举个例子,如果我们想要求出e的2次方的反函数,即e^2的反函数,可以使用以下代码:
x = exp(2);
y = log(x);
执行完以上代码后,y的值为2,即e^2的反函数为2。
- 1 -。