eigen value of stochastic matrix
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I wan to find eigen values and eigen vector of a stochastic matrix using Newton Raphson method.
% Define the dimensions and variables
n = 3; % Size of vectors and matrices
m = 2; % Number of matrices
P = 4; % Number of eigenvalues and tensors
maxIterations = 100; % Maximum number of iterations
tolerance = 1e-6; % Tolerance for convergence
% Initialize variables
u = ones(n*P, 1); % Initial guess for u, each column represents u1, u2, ..., uP
lambda = ones(P, 1); % Initial guess for lambda, each element represents lambda1, lambda2, ..., lambdaP
% Define the matrices and tensors
c0 = eye(P); % Example: identity matrix
A0 = magic(n); % Example: identity matrix
c = rand(P, P, m); % Example: random P x P x m tensor
A = rand(n, n, m); % Example: random n x n x m tensor
e = rand(P, P, P); % Example: random P x P x P tensor
cA=0;
for i=1:m
cA= cA+kron(c(:,:,i),A(:,:,i));
end
% Perform Newton-Raphson iteration
for iter = 1:maxIterations
% Compute the left-hand side of the equation
lhs = (kron(c0 , A0) + cA) * u;
% Compute the right-hand side of the equation
lame=0;
for p = 1:P
lame = lame+kron(lambda(p) * e(:, :, p),eye(n)) ;
end
rhs=lame*u
% Compute the residual and check for convergence
residual_1 = lhs - rhs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the variables
I=eye(n,n);
cA1=zeros(P,1);
for i=1:P
cA1(i)= u'*kron(e(:,:,i),I)*u;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
residual_2=cA1;
% Compute the residual and check for convergence
residual = [residual_1;residual_2];
if norm(residual(:)) < tolerance
disp('Converged!');
break;
end
jacobian = jac_1(A,lambda , u, n,m,P);
% Update u and lambda using the Newton-Raphson method
delta = jacobian \ residual(:);
u = u - delta(1:n*P);
lambda = lambda - delta(n*P+1:end);
% Display the current iteration and residual
disp(['Iteration: ' num2str(iter) ', Residual: ' num2str(norm(residual(:)))]);
end
uf=reshape(u,n,[])
lambdaf=lambda
but the answer I am getting is
lambdaf =
1.0e+58 *
2.7709
-1.9052
0.3483
-0.3253.
It would be very helpful if someone help me improve this code
2 Kommentare
Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!