Find Eigenvalues of ODE45 Solution MATLAB

I have the following non-linear ODE:

I have the following ODE45 solution:

fun = @(t,X)odefun(X,K,C,M,F(t),resSize); 
[t_ode,X_answer] = ode45(fun,tspan,X_0); 

The input matrices are stiffness K(X), damping C, mass M, and force F. resSize is the total number of masses in the system.

I would like to find the system's eigenvalues using either the Jacobian matrix, transfer function, or any other viable method.

I have tried using:

[vector,lambda,condition_number] = polyeig(K(X_answer),C,M);

This is tricky since my K matrix is a function handle of X. In other words, K=@(X). X represents a displacement vector of each mass in the system (x_1(t),x_2(t),...x_resSize(t)), where resSize is the total number of masses. My X_answer matrix is a double with dimensions of t_ode by resSize, where each row is the displacement vector of each mass in double form. Is there some way to substitute X_answer into my function handle for K so I can use polyeig()? If not, how would I go about finding my system's transfer function or Jacobian matrix so that I can find it's eigenvalues?

[Sam Chak]

Is a state-dependent matrix, as illustrated below?

If the system is not in equilibrium, the ode45 solver will provide a solution array 'X_answer' that corresponds to the values returned in the time vector 't_ode'. In this case, the stiffness matrix should vary over time.

Would you like to compute the eigenvalues at each time step from tspan(1) to tspan(end)?

[Jonathan Frutschy]

Yes. the stiffness matrix varies over time. (see edited post for the EOM). Yes. I would like to calulate the eigenvalues at each timestep from tspan(1) to tspan(end). @Sam Chak

[Jonathan Frutschy]

@Sam Chak Just as a follow up, is this how you would compute the eigenvalues at each timestep?

for ii=1:length(t_ode)
% This is a QUADRATIC eigenvalue problem, so you will have 2*resSize
% eigenvalues at each timestep, not resSize eigevalues !!!
[vector{ii},lambda{ii}] = polyeig(K(X_answer(ii,1:end/2)),C,M); % F(t) is accounted for in X_answer (i.e. the system response)
end

[Sam Chak]

@Jonathan Frutschy, while the syntax "polyeig(K, C, M)" appears correct, I never compute eigenvalues at each time step as eigenvalues are meaningful only for Linear Time-Invariant (LTI) systems. However, ensure that you perform the linearization correctly as described in @Torsten's answer.

What do you want to analyze from the array of eigenvalues?

[Jonathan Frutschy]

@Sam Chak I would like to find the system’s eigenfrequencies at each time step (omega_natural = (lambda_natural)^1/2). That way I have some insight into what frequencies to avoid stimulating the system with to prevent resonance.

[Sam Chak]

Hi @Jonathan Frutschy

The input force consists of the sum of a slower wave and a faster wave.

I'm curious about the scientific basis or any journal paper that describes the method of determining the system's resonance by computing the state-dependent eigenvalues at each time step. Is this related to the design of a Tuned Mass Damper for a high-rise building?

[Jonathan Frutschy]

@Sam Chak This is intended for a classified research project. While I can't divulge the details, a tuned mass damper on a high rise building is a great application!

Answer 1

[Torsten]

Your system reads

z1'= z2
z2' = ( F-(c1+c2)*z2-( ((k1*gamma1+k2*gamma2)*z1^2+k1+k2)*z1+(-k2*gamma2*z1^2-k2)*z3 ) )/m1
z3' = z4
z4' = (-(c2+c3)*z4-( (-k2*gamma2*z3^2-k2)*z1+((k2*gamma2+k3*gamma3)*z3^2+k2+k3)*z3))/m2

Linearize the right hand side by computing the 4x4 Jacobian (with respect to z1 - z4) and compute its eigenvalues in the course of the simulation.