eigen value of stochastic matrix

Question

rakesh kumar 2023-6-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1978609-eigen-value-of-stochastic-matrix

评论： rakesh kumar 2023-6-8

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 个评论
显示无隐藏无

Kautuk Raj 2023-6-8

Can you tell us what answer are you expecting?

rakesh kumar 2023-6-8

在 MATLAB Online 中打开

the lambdaf and uf values should be small.

tlength=1;
nel=2;
kmn=1;
w=0.95;
sigma=0.6;
E0=2.0*10^11;
el=E0*ones(1,nel);
[alpha,kl]=klterm4_spanos3(w,sigma,kmn,nel);
[kk,mm]=stiffm_new(nel,el);
kk=kk(3:end,3:end);
[ka,ma]=stiffm(nel,1);
m=ma(3:end,3:end);
for t=1:kmn
for j=1:kl
[kk2(:,:,j)]=stiffm_new(nel,E0*alpha(j,:,t));
kk1(:,:,j)=kk2(3:end,3:end,j);
end
end
k00=reshape(kk1,2*nel,[]);
k=[kk k00];
k=reshape(k,2*nel,2*nel,[]);
M=kl;
p=1;
A0=k(:,:,1);
A=k(:,:,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the dimensions and variables
n = 2*nel; % Size of vectors and matrices
%M = 2; % Number of matrices
%p = 2; % order of polynomials
P=factorial(M+p)/(factorial(M)*factorial(p));
m1 = 1; % m1th eigen value
maxIterations = 1; % Maximum number of iterations
tolerance = 1e-6; % Tolerance for convergence
% Initialize variables
u1 = ones(n*P, 1); % Initial guess for u, each column represents u1, u2, ..., uP
lambda1 = ones(P, 1); % Initial guess for lambda, each element represents lambda1, lambda2, ..., lambdaP
method=1;
c00= c_ijk(M,p,method);
e= e_ijk(M,p);
c_ijk_mat=my_cell2mat(c00,M,p);
e1=my_cell2mat_e(e,M,p);
% Define the matrices and tensors
%c0 = eye(P); % Example: identity matrix
c0 = c_ijk_mat(:,:,1);
%A0 = magic(n); % Example: identity matrix
%c = rand(P, P, M); % Example: random P x P x m tensor
c = c_ijk_mat(:,:,2:end);
%A = 10^-2*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
e= e1;
[newLambda, newU] = replaceLambdaAndU(m1,P,A0);
lambda=newLambda;
u=reshape(newU,[],1);
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
I want first elemnets of lambdaf vetor to be smallest value of eigen value of A0 and all other values of lambdaf to be smaller values than lambdaf ,of the order of 10^-2*lambdaf

请先登录，再进行评论。

请先登录，再回答此问题。

eigen value of stochastic matrix

2 个评论
显示无隐藏无

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

eigen value of stochastic matrix

2 个评论 显示 无隐藏 无

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无