Application Center - Maplesoft

App Preview:

Eigenvalue/Eigenvector Sensitivity

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Image 

Eigenvalue/Eigenvector Sensitivity 

Hakan Tiftikci
Turkey
hakan.tiftikci@yahoo.com.tr 

Introduction 

In this worksheet, a method to compute sensitivities Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mo(of eigenvalues/eigenvectors of a matrix Typesetting:-mrow(Typesetting:-mi(with respect to parameters of the matrix Typesetting:-mrow(Typesetting:-mi( is presented. Method is based on the Cayley-Hamilton theorem, which states that, the matrix Typesetting:-mrow(Typesetting:-mi(itself satisfies its characteristic polynomial, i.e. 

 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

 

where 

 

Typesetting:-mrow(Typesetting:-mi( 

 

is the characteristic polynomial of matrix Typesetting:-mrow(Typesetting:-mi(. Factored form of characteristic polynomial is given by 

 

Typesetting:-mrow(Typesetting:-mi( 

 

which have a connection to algebraic form by symmetric functions of eigenvalues by 

 

Typesetting:-mrow(Typesetting:-mo( 

 

Roots Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( are the eigenvalues of matrix Typesetting:-mrow(Typesetting:-mi(. 

 

Cayley-Hamilton theorem can be expressed using the factored form of characteristic polynomials, viz. 

 

Typesetting:-mrow(Typesetting:-mi( 

 

Taking partial derivative of charactristic polynomial Typesetting:-mrow(Typesetting:-mi(with respect to a parameter Typesetting:-mrow(Typesetting:-mi( that the matrix Typesetting:-mrow(Typesetting:-mi( depends on it, 

 

Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mo(
Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mo(
 

 

or if written more concisely 

 

Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mo( 

 

Denoting 

 

Typesetting:-mrow(Typesetting:-munderover(Typesetting:-mo( 

 

the last equation becomes 

 

Typesetting:-mrow(Typesetting:-mi( 

 

expanding and isolating eigenvalue derivatives on one side, following final matrix equation is obtained 

 

Typesetting:-mrow(Typesetting:-munderover(Typesetting:-mo( 

 

which is Typesetting:-mrow(Typesetting:-msup(Typesetting:-mi( equations in Typesetting:-mrow(Typesetting:-mi( unknown eigenvalue derivatives Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(Selecting Typesetting:-mrow(Typesetting:-mi(appropriate (well-formed/nonsingular) equations from a total of Typesetting:-mrow(Typesetting:-msup(Typesetting:-mi(, eigenvalue derivatives Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(can be obtained. For the derivatives Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mo(of eigenvectors Typesetting:-mrow(Typesetting:-mmultiscripts(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(, the defining equation 

 

Typesetting:-mrow(Typesetting:-mi( 

 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi( 

 

is differentiated with respect to Typesetting:-mrow(Typesetting:-mi( to yield 

 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

 

Initialization 

Implementation and Application 

In the first subsection, method is first illustrated step-by-step for a simple system and then in second part is put to procedural form.  

 

Application to a simple system (Spring-Mass-Damper) 

When spring-mass-damper system Typesetting:-mrow(Typesetting:-mi(is converted to state-space form by Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( matrix Typesetting:-mrow(Typesetting:-mi(in linear representation Typesetting:-mrow(Typesetting:-mi(is given by following 

> A := <<0, -omega^2>|<1, -2*omega*zeta>>;
 

Typesetting:-mrow(Typesetting:-mi( (3.1.1)
 

Characteristic polynomial is given by 

> CharacteristicPolynomial(A, 's');
 

`+`(`*`(`^`(s, 2)), `*`(2, `*`(omega, `*`(zeta, `*`(s)))), `*`(`^`(omega, 2))) (3.1.2)
 

Eigenvalues are 

> Lambda := Eigenvalues(A);
 

Typesetting:-mrow(Typesetting:-mi( (3.1.3)
 

Matrix satisfies Cayley-Hamilton theorem,  

 

Typesetting:-mrow(Typesetting:-msup(Typesetting:-mi( 

 

as shown by following computation 

> J2 := Matrix(2,2,shape=identity);
 

Typesetting:-mrow(Typesetting:-mi( (3.1.4)
 

> testCayleyHamilton := A.A+MatrixScalarMultiply(A, 2*omega*zeta)+omega^2*J2;
 

Typesetting:-mrow(Typesetting:-mi( (3.1.5)
 

Sensitivities of eigenvalues (by direct computation) are 

> Mu1 := map(diff, Lambda, omega);
Mu2 := map(diff, Lambda, zeta);
 

 

Typesetting:-mrow(Typesetting:-mi(
Vector[column](%id = 149605472) (3.1.6)
 

Derivatives of the matrix with respect to its parameters 

> A1 := map(diff,A,omega);
A2 := map(diff,A,zeta);
 

 

Typesetting:-mrow(Typesetting:-mi(
Matrix(%id = 149626672) (3.1.7)
 

Derivative with respect to first parameter Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

Denote derivatives of Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( by Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(, then left and right hand side of  E.1 and the matrix equation are 

> LHS := mu1*(A-Lambda[2]*J2)+mu2*(A-Lambda[1]*J2):
RHS := A1.(A-Lambda[2]*J2)+(A-Lambda[1]*J2).A1:
EQU := LHS-RHS;
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
(3.1.1.1)
 

solving by Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn(elements of matrix equation 

> sol_11_22 := solve({EQU[1, 1], EQU[2, 2]}, {mu1, mu2});
 

{mu1 = `+`(`-`(zeta), `*`(`^`(`+`(`*`(`^`(zeta, 2)), `-`(1)), `/`(1, 2)))), mu2 = `+`(`-`(zeta), `-`(`*`(`^`(`+`(`*`(`^`(zeta, 2)), `-`(1)), `/`(1, 2)))))} (3.1.1.1.1)
 

compare last result with directly computed results 

> eval(<mu1,mu2>,sol_11_22)-Mu1;
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (3.1.1.1.2)
 

solving by Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn( 

> sol_11_22 := solve({EQU[1, 1], EQU[2, 2]}, {mu1, mu2});
 

{mu1 = `+`(`-`(zeta), `*`(`^`(`+`(`*`(`^`(zeta, 2)), `-`(1)), `/`(1, 2)))), mu2 = `+`(`-`(zeta), `-`(`*`(`^`(`+`(`*`(`^`(zeta, 2)), `-`(1)), `/`(1, 2)))))} (3.1.1.2.1)
 

compare last result with directly computed results 

> eval(<mu1,mu2>,sol_11_22)-Mu1;
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (3.1.1.2.2)
 

Performing computations for the Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn(and Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn(elements,  equivalence of solution obtained from different element selections of matrix equation (E.1) are demostrated (just) for the case. 

 

Derivative with respect to second parameter Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

Denote derivatives of Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( by Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(, then left and right hand side of  E.1 and the matrix equation are 

> LHS := mu1*(A-Lambda[2]*J2)+mu2*(A-Lambda[1]*J2):
RHS := A2.(A-Lambda[2]*J2)+(A-Lambda[1]*J2).A2:
EQU := LHS-RHS;
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
(3.1.2.1)
 

solving by Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn( and Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn( elements of matrix equation 

> sol_11_12 := solve({EQU[1, 1], EQU[1, 2]}, {mu1, mu2});
 

{mu2 = `+`(`-`(`/`(`*`(`+`(`*`(zeta, `*`(`^`(`+`(`*`(`^`(zeta, 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta, 2)), `-`(1)), `*`(omega)), `*`(`+`(`*`(`^`(zeta, 2)), `-`(1)))))), mu1 = `+`(`-`(`/`(`*`(omega,... (3.1.2.2)
 

compare last result with directly computed results 

> eval(<mu1,mu2>,sol_11_12)-Mu2;
simplify(%);
 

 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mo(
Vector[column](%id = 152601912) (3.1.2.3)
 

Sensitivity Procedure  

In this section, the method whose step are explained in above example is defined as a Maple procedure. 

> EigenSensitivity := proc(A,q,L)
local LHS,RHS,N,i,Aq,AP,MATEQ,sol;
description "computes the partial derivatives of eigenvalues 'L' of matrix 'A' with respect to parameter 'q'";

# if not given, compute the eigenvalues
#if evalb(L=NULL) then
# L := Eigenvalues(A);
#end if;

N := RowDimension(A);  # system size

# derivative of matrix with respect to parameter "q"
Aq := map(diff,A,q);

# define auxiliary product function
AP := proc(a,b)
local i,result,J;
description "auxiliary product function for defining matrix equation";
J := Matrix(N,N,shape=identity);
result := Matrix(N,N,shape=identity);
for i from a to b do
result := result . (A-L[i]*J);
end do;
result;
end proc;

# construct the left and right side of matrix equation (E.1)
LHS := Matrix(N,N,0);
RHS := Matrix(N,N,0);
for i from 1 to N do
LHS := LHS + AP(1,i-1) . Aq . AP(i+1,N);
RHS := RHS + MatrixScalarMultiply( AP(1,i-1) . AP(i+1,N), (mu||i));

end do;
# form matrix equation
MATEQ := LHS-RHS;

# solve matrix equation for eigenvalue derivatives
sol := solve( {seq(MATEQ[i,i],i=1..N)}, {mu||(1..N)});
# as default, linear system is selected from diagonals of matrix equation

# return (simplified) result as vector
#map(simplify,eval(<mu||(1..N)>,sol));
eval(<mu||(1..N)>,sol);
end proc:
 

> #untrace(EigenSensitivity);
 

 

Sensitivity Procedure Test 

Test the procedure on Spring-Mass-Damper 

> EigenSensitivity(A,omega,Lambda);
EigenSensitivity(A,zeta,Lambda);
simplify(%%-Mu1);
simplify(%%-Mu2);
 

 

 

 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mo(
Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mo(
Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn(
Vector[column](%id = 148249848) (3.3.1.1)
 

Test the procedure on general 2x2 matrix 

> A2 := <<a11,a21>|<a12,a22>>;
 

Typesetting:-mrow(Typesetting:-mi( (3.3.2.1)
 

> Lambda2 := Eigenvalues(A2);
 

Typesetting:-mrow(Typesetting:-mi( (3.3.2.2)
 

derivatives with respect to Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

compute sensitivities directly 

> Mu_11_direct := map(diff, Lambda2, a11);
 

Typesetting:-mrow(Typesetting:-mi( (3.3.2.1.1)
 

compute sensitivities by the method 

> Mu_11_method := EigenSensitivity(A2,a11,Lambda2);
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
(3.3.2.1.2)
 

> simplify(Mu_11_direct-Mu_11_method);
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (3.3.2.1.3)
 

derivatives with respect to Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

compute sensitivities directly 

> Mu_12_direct := map(diff, Lambda2, a12);
 

Typesetting:-mrow(Typesetting:-mi( (3.3.2.2.1)
 

compute sensitivities by the method 

> Mu_12_method := EigenSensitivity(A2,a12,Lambda2);
 

Typesetting:-mrow(Typesetting:-mi( (3.3.2.2.2)
 

> simplify(Mu_12_direct-Mu_12_method);
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (3.3.2.2.3)
 

Test the procedure on simplified 3x3 matrix 

For this case eigenvalues used in obtaining their derivatives are kept as symbols until last step in which results are compared to those of direct computation. Otherwise computational resources are quickly consumed. For the same reason, in this case, matrix is simplified by fixing some elements and reducing number of parameters. 

> A3 := Matrix(3,3, (i,j) -> a||i||j, shape=symmetric):
A3[2,1]:=1:
A3[3,1]:=2:
A3[3,2]:=2*a22:
A3;
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mi( (3.3.3.1)
 

> Lambda3 := Eigenvalues(A3):
 

derivatives with respect to Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( 

compute sensitivities directly 

> Mu_22_direct := map(diff, Lambda3, a22):
 

compute sensitivities by the method 

> Mu_22_method := EigenSensitivity(A3,a22,<L3_1,L3_2,L3_3>):
 

> mustBeZero := subs([L3_1=Lambda3[1],L3_2=Lambda3[2],L3_3=Lambda3[3]],Mu_22_method)-Mu_22_direct:
 

> simplify(mustBeZero);
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (3.3.3.1.1)
 

Eigenvalue Jacobian and a Related Application 

Jacobian of eigenvalue vector is obtained by repeated application of the method ("EigenSensitivity"). A more efficient version may be obtained by identifying the repeated multiplications and/or common factors.

Eigenvalue Jacobian can be used to determine required direction/rate in parameter space for attaining a direction/rate in eigenvalue space, as formulated by following expressions
 

 

Typesetting:-mrow(Typesetting:-mi( 

Typesetting:-mrow(Typesetting:-mi( 

or 

 

Typesetting:-mrow(Typesetting:-mi( 

Typesetting:-mrow(Typesetting:-mi( 

 

if it is desired to direct/move current eigenvalue vector Typesetting:-mrow(Typesetting:-mi( towards a target eigenvalue vector Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(then the choice  

 

Typesetting:-mrow(Typesetting:-mi( 

 

serves to this purpose. 

> EigenJacobian := proc(A,Q,L)
local result,i,Lq;
description "computes the Jacobian of matrix 'A' with respect to parameters 'Q'";
result := Matrix( RowDimension(A), Dimension(Q));
for i from 1 to Dimension(Q) do
Lq := EigenSensitivity(A,Q[i],L);
result[1..RowDimension(A),i] := Lq;
end do;
result;
end proc:
 

 

Application to Spring-Mass-Damper 

Above formulation is applied to move the parameters Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi(from an arbitrary starting point to a desired target point. First Jacobian is computed as 

> J := EigenJacobian(A,<omega,zeta>,Lambda);
 

Typesetting:-mrow(Typesetting:-mi( (3.4.1.1)
 

> Jinv := MatrixInverse(J);
 

Typesetting:-mrow(Typesetting:-mi( (3.4.1)
 

define gain (positive definite) and target eigenvalues (yet as symbols) 

> K := <<k11,k21>|<k12,k22>>;
Ltarget := <L1,L2>;
 

 

Typesetting:-mrow(Typesetting:-mi(
Vector[column](%id = 161345644) (3.4.2)
 

express differential vector in parameter space that corresponds to Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( direction in eigenvalue space. 

> dp := Jinv . K . (Ltarget-Lambda);
 

Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mi(
(3.4.3)
 

study-1 

Eigenvalues are moved to Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mo( 

> SL := [L1=-3,L2=-4,k11=1,k12=0,k21=0,k22=1,zeta=zeta(t),omega=omega(t)];
 

[L1 = -3, L2 = -4, k11 = 1, k12 = 0, k21 = 0, k22 = 1, zeta = zeta(t), omega = omega(t)] (3.4.2.1)
 

> DEQ := diff(omega(t),t)=subs(SL,dp[1]),
diff(zeta(t),t)=subs(SL,dp[2]);
 

diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`-`(3), `-`(`*`(`+`(`-`(zeta(t)), `*...
(3.4.2.2)
 

solve differential system numerically (note complex option) 

> dsol := dsolve({DEQ,zeta(0)=1/10,omega(0)=1},{zeta(t),omega(t)},'numeric',complex=true);
 

proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); if... (3.4.2.3)
 

Plot locus in parameter space 

> zetaplot := plots[odeplot](dsol,[Re(zeta(t)),Im(zeta(t))],t=0..5,style=POINT,color=red):
omegaplot := plots[odeplot](dsol,[Re(omega(t)),Im(omega(t))],t=0..5,style=POINT,color=blue):
 

> plots[display]({zetaplot,omegaplot},axes=boxed,labels=["Re","Im"],title="parameters in complex plane");
 

Plot_2d
 

plot history of eigenvalues (compare to target values) 

> solutionEigen1 := subs(omega=omega(t),zeta=zeta(t),Lambda[1]);
solutionEigen2 := subs(omega=omega(t),zeta=zeta(t),Lambda[2]);
 

 

`*`(`+`(`-`(zeta(t)), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(omega(t)))
`*`(`+`(`-`(zeta(t)), `-`(`*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2))))), `*`(omega(t))) (3.4.2.4)
 

> plots[odeplot](dsol,[t,Re(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L1)"]);
plots[odeplot](dsol,[t,Im(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L1)"]);
plots[odeplot](dsol,[Re(solutionEigen1),Im(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L1)","Im(L1)"]);
 

 

 

Plot_2d
Plot_2d
Plot_2d
 

> plots[odeplot](dsol,[t,Re(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L2)"]);
plots[odeplot](dsol,[t,Im(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L2)"]);
plots[odeplot](dsol,[Re(solutionEigen2),Im(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L2)","Im(L2)"]);
 

 

 

Plot_2d
Plot_2d
Plot_2d
 

study-2 

eigenvalues are moved to Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mo( 

> SL := [L1=-3+2*I,L2=-4,k11=1,k12=0,k21=0,k22=1,zeta=zeta(t),omega=omega(t)];
 

[L1 = `+`(`-`(3), `*`(2, `*`(I))), L2 = -4, k11 = 1, k12 = 0, k21 = 0, k22 = 1, zeta = zeta(t), omega = omega(t)] (3.4.3.1)
 

> DEQ := diff(omega(t),t)=subs(SL,dp[1]),
diff(zeta(t),t)=subs(SL,dp[2]);
 

diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
diff(omega(t), t) = `+`(`-`(`/`(`*`(`/`(1, 2), `*`(`+`(`*`(zeta(t), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(`^`(zeta(t), 2)), `-`(1)), `*`(`+`(`+`(`-`(3), `*`(2, `*`(I))), `-`(`*`...
(3.4.3.2)
 

> dsol := dsolve({DEQ,zeta(0)=1/10,omega(0)=1},{zeta(t),omega(t)},'numeric',complex=true);
 

proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; `:=`(_EnvDSNumericSaveDigits, Digits); `:=`(Digits, 14); if... (3.4.3.3)
 

> zetaplot := plots[odeplot](dsol,[Re(zeta(t)),Im(zeta(t))],t=0..1.2,numpoints=64,style=POINT,color=red):
omegaplot := plots[odeplot](dsol,[Re(omega(t)),Im(omega(t))],t=0..1.2,numpoints=64,style=POINT,color=blue):
 

> plots[display]({zetaplot,omegaplot},axes=boxed,labels=["Re","Im"],title="parameters in complex plane");
 

Plot_2d
 

> solutionEigen1 := subs(omega=omega(t),zeta=zeta(t),Lambda[1]);
solutionEigen2 := subs(omega=omega(t),zeta=zeta(t),Lambda[2]);
 

 

`*`(`+`(`-`(zeta(t)), `*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2)))), `*`(omega(t)))
`*`(`+`(`-`(zeta(t)), `-`(`*`(`^`(`+`(`*`(`^`(zeta(t), 2)), `-`(1)), `/`(1, 2))))), `*`(omega(t))) (3.4.3.4)
 

> plots[odeplot](dsol,[t,Re(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L1)"]);
plots[odeplot](dsol,[t,Im(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L1)"]);
plots[odeplot](dsol,[Re(solutionEigen1),Im(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L1)","Im(L1)"]);
 

 

 

Plot_2d
Plot_2d
Plot_2d
 

> plots[odeplot](dsol,[t,Re(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L2)"]);
plots[odeplot](dsol,[t,Im(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L2)"]);
plots[odeplot](dsol,[Re(solutionEigen2),Im(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L2)","Im(L2)"]);
 

 

 

Plot_2d
Plot_2d
Plot_2d
 

Note that eigenvalue locus is straight line as imposed while developing differential relations. 

References 

[1] Cayley-Hamilton@Wikipedia 

[2] Cayley-Hamilton@PlanetMath 

Conclusions 

A general method to compute sensitivities of eigenvalues/eigenvectors is described. Method is applied to simple Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn(system and simplified Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn(system. Even though, the method described does not necessarily rely on symbolic computation, CAS environment provides robust way of implementing and verifying results. 

 

The described method, since involves many matrix multiplications,  is computationally intensive, and for the symbolic case, eigenvalues (if available) are generally so complex that, resulting matrix equation bears all the complexity of eigenvalues involved, causing deficiency in computer resources even though equations are linear in eigenvalue derivatives to be determined. For such cases eigenvalues may be left as symbols until their values are needed (as illustrated in Typesetting:-mrow(Typesetting:-mn(test case) 

 

Additionaly,  following properties of method are observed: 

  • All the eigenvalues must be available. Hence, pure symbolic solutions are limited by degree 5 (Galois) for the general case.
 

  • Formulation has the structure of explicit ODE Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(. Thus, for continuous change of parameters, corresponding eigenvalues can be traced by integration of Typesetting:-mrow(Typesetting:-mi(.
 

  • Derivatives of the matrices with respect to parameters must be available. In pure numeric implementation, this requirement poses another problem. Either numerical (finite difference,..), Automatic-Differentiation (AD), or CAS generated code for differentiation may be used.
 

 


Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.
 

 

Image