Find State Space Representation of Linearized non-linear System to find State Space Representation in Canonical form MATLAB

I have the following non-linear system and accompanying 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. To find the solution indirectly I introduced a change of variables to reduce the non-linear ODE to a first order linear ODE.

I would like to find the state space (SS) representation of my linearized non-linear system and corresponding ODE. More specifically, I would like to find the SS matrices A,B,C,D. I did not explicitly perform the change of variables. Rather, MATLAB did it for me. The ultimate goal of finding the SS representation (matrices A,B,C,D) is to then be able to find the canonical form of the non-linear ODE. I'm not sure if this can be done with a transfer function estimate of my ODE45 solution or some other MATLAB tool.

Answer 1

[Sam Chak]

Hi @Jonathan Frutschy

Given that no original MATLAB code for the nonlinear system is provided, I believe it would be more efficient for me to derive the linear state-space matrices intuitively by hand. If my derivation is correct, you should obtain the following symbolic state matrix :

syms m1 m2 c1 c2 c3 k1 k2 k3

%% symbolic state matrix
As  = [zeros(2), eye(2);
       -(k1+k2)/m1, k2/m1, -(c1+c2)/m1, 0;
       k2/m2, -(k2+k3)/m2, 0, -(c2+c3)/m2]

As = 

%% symbolic input matrix
Bs  = [sym('0'); sym('0'); sym('1'); sym('0')]

Bs = 

%% State-space model
A   = matlabFunction(As);
B   = matlabFunction(Bs);
C   = eye(4);
D   = 0;

m1  = 1;
m2  = 1;
c1  = 1;
c2  = 1; 
c3  = 1; 
k1  = 1; 
k2  = 1; 
k3  = 1;
sys = ss(A(c1,c2,c3,k1,k2,k3,m1,m2), B(), C, D)

sys =
 
  A = 
       x1  x2  x3  x4
   x1   0   0   1   0
   x2   0   0   0   1
   x3  -2   1  -2   0
   x4   1  -2   0  -2
 
  B = 
       u1
   x1   0
   x2   0
   x3   1
   x4   0
 
  C = 
       x1  x2  x3  x4
   y1   1   0   0   0
   y2   0   1   0   0
   y3   0   0   1   0
   y4   0   0   0   1
 
  D = 
       u1
   y1   0
   y2   0
   y3   0
   y4   0
 
Continuous-time state-space model.

%% Step response
step(sys), grid on

[Ibrahim]

Hi @Sam Chak, I have a question that didn't find an answer to it, can you tell me how to contact you for further information. The Matlab community closed my question. I can pay you to solve the question, let me know how to contact you, Thanks.

The question: "https://www.mathworks.com/matlabcentral/answers/2107781-write-code-and-driver-matlab-programs-to-solve-the-previous-combined-differential-equations-to-simu"

[Jonathan Frutschy]

@Sam Chak Sam, this would work if the springs were linear, but they are non-linear. The force in each i^th spring is: F_i = k_i*x_i+gamma_i*(x_i)^3. That Is why I used a change of variables in my stiffeness matrix by factoring out x_i to make the stiffness matrix a function of x_i.

[Sam Chak]

Hi @Jonathan Frutschy

I believe there might be some confusion regarding your request to find the state-space matrices , , , . These matrices are typically associated with linear systems or linearized systems , . Nonlinear systems, on the other hand, do not have , , , matrices.

In the case of your mass-damper-spring system, it is nonlinear with respect to the states but linear with respect to the input signal F (often denoted as 'u' in standard terminology). Therefore, the general state-space form of an affine nonlinear system can be expressed as:

.

In this equation, represents a nonlinear function, while and behave linearly, similar to the input matrix and the output matrix , respectively.

To ensure clear understanding for everyone, would you consider revising the description of your problem above?

[Ibrahim]

Hi @Sam Chak, How I contact you?

[Sam Chak]

Hi @Ibrahim, I wanted to mention that it would be inappropriate for me to respond to your request within @Jonathan Frutschy's thread, especially since it is unrelated to his modeling problem for solving in the ode45 solver.

I didn't close your question. If you believe your question was closed without a valid reason, you could try reopening your question to provide additional context and to explain why the question is relevant to you and MATLAB Answers community.

Here is an example:

• "I have encountered a (particular) issue (or error) while attempting to solve this problem in MATLAB/Simulink. Here is my initial attempt, which was inspired by this example from the MathWorks documentation (provide the URL link). Could someone please point out what the issue might be and guide me regarding which MATLAB functions I should use to complete (this specific) task?"

However, it is important to note that using an instructional tone, such as "Write the MATLAB code..." or "Build the Simulink model..." is generally discouraged in this forum. Questions phrased in this manner may be closed by higher-ranked users, as they do not align with the guidelines set forth by MATLAB Answers.

[Jonathan Frutschy]

@Sam Chak This clear things up. I have revised the question to emphasize that the non-linear system must be linearized (via a change of variables in this case) to find a SS representation. I guess I'm just confused on why your A matrix doesn't have any gamma terms in it?

[Sam Chak]

Hi @Jonathan Frutschy

There are two key observations that can be made from the system diagram and the equations of motion image. These observations allow me to "intuitively linearize" the nonlinear system without requiring your MATLAB code. However, if you do provide the code for the nonlinear system, I can offer guidance using MATLAB tools.

1. The Equilibrium point

When there is no force acting () and the system is at rest (), the equilibrium positions of masses and are and , respectively.

2. The stiffness matrix, K

Furthermore, if the nonlinear system is linearized around the equilibrium positions ( & ), the influence of the parameter γ in the stiffness matrix diminishes. In other words, when the motion operates within a range of conditions centered around the equilibrium positions, the nonlinear effect of γ becomes negligible. This is why the parameter γ does not appear in the linear state matrix mentioned in my Answer above.

[Jonathan Frutschy]

@Sam Chak Thank you for clearing that up. See below for my complete implementation.

The functions carray,marray, and karray were custom written and return the K,C, and M matrices in the ODE shown in the original post for any number of masses (ie. resSize). The function farray() is currently set up to apply the force input to every 100th mass. Right now my biggest issue is that my mass displacements are very small, with some being as small as (10^-40 meters). MATLAB sometimes will also not be able to solve the ODE and simply crash.

Looking into stiff solvers like ode15s, it states that providing the Jacobian matrix in the odeset() option will improve performance. Given the code below, do you know how I can find the Jacobian matrix and specify it in odeset()? The only other thing I can think of is contuing to play with my k,c,m,and gamma values and the amplitude of the force input until I get some more reasonable displacement numbers.

resSize = 600; % number of masses
m = 1.*ones(resSize,1); % [kg]
k = 2.*ones(resSize+1,1); % [N/m]
c = 1.*ones(resSize+1,1); %[kg/s]
gamma = 1.*ones(resSize+1,1); % [N/m^3]

alpha = 0.2;
beta = 0.15;
dt = 1; % [s]
tspan = [0:dt:20]; % [s]. The units of time are dictated by the units of odefun().
X_0 = zeros(resSize*2,1);
f_0 = 0.007; % [N] input signal amplitude 
f = @(t) f_0.*sin(alpha.*t).*cos(beta.*t);
F = farray(f,resSize);
K = karray(k,gamma,resSize);
M = marray(m);
C = carray(c,resSize);
% Solving with ODE45
fun = @(t,X)odefun(X,K,C,M,F(t),resSize); % Defines the ODE function and it's input matrices
[t_ode,X_answer] = ode45(fun,tspan,X_0); % Calls the ODE45 function, sets the time and axial conditions. Outputs time and displacements.

[Torsten]

Why do you still solve the linearized system if you now switched to "ode45" instead of "ss" ?

[Sam Chak]

Hi @Jonathan Frutschy

The nonlinear system is unavailable in your code. So I tried solving a simple system with 600 masses with ode45 and it works well. I believe there must be another reason that your MATLAB crashed.

%% State-space method (Unforced)
n       = 600;      % number of masses
A       = [zeros(n),                               eye(n);
          -diag(1*ones(1, n))\diag(1*ones(1, n)), -diag(1*ones(1, n))\diag(2*ones(1, n))];
[t, x]  = ode45(@(t, x) A*x, [0 20], [1; zeros(2*n-1, 1)]);
plot(t, x), grid on, xlabel t, ylabel('\bf{x}\rm{(t)}')

[Paul]

@Jonathan Frutschy

Though not showing the code for odefun (or the other functions), I don't understand how it's used. If i understand correctly, and I might not, when the anonymous function fun is defined, aren't K and C just numeric matrices? If so, then the values of K and C are fixed for the entire simulation. Is that what's intended? Or are K and C both functions, in which case they should be K(X) and C(x) as inputs to odefun? Though I guess K(x) and C(x) can be evalutated inside odefun. Is F, as returned from f_array, an anonymous function?

Can you post the entire code, maybe just for a few masses?

[Jonathan Frutschy]

@Torsten I've been using ode45() the entire time. I was thinking of using ode15s() as another option.

[Jonathan Frutschy]

@Thales @Paul @Sam Chak To help everyone out, here are the missing functions:

function output = odefun(X,K,C,M,F,resSize) % odefun, with input matrices
 X_variable = X(1:resSize); % position
 X_variable_dot = X(resSize+1:2*resSize); % 1st derivitive
 X_variable_dot_dot = M\(F-K(X)*X_variable-C*X_variable_dot); % 2nd derivitive, from the equation of motion
 output=[X_variable_dot;X_variable_dot_dot]; % places results as a column vector
end

function F = farray (f,N)
     % for 600 masses
     F = @(t) [f(t); zeros(175,1); f(t); zeros(100,1); f(t); zeros(74,1); f(t); zeros(175,1); f(t); zeros(71,1)];
end

function K = karray(k, gamma, N)
    K = @(x) make_diagonal(x, k, gamma, N);
end

function out = make_diagonal(x, k, gamma, N)
    x = x(:);
    x = [x(1:N); 0];

    ck = circshift(k, -1);
    cg = circshift(gamma, -1);
    cx = circshift(x, -1);
    ccx = circshift(x, 1);

    d1 = -3 .* ck .* cg .* cx .^ 2 - ck;
    d2 = (k .* gamma + ck .* cg) .* x .^ 2 + k + ck;
    d3 = -3 .* k .* ccx .^ 2 - k;

    out = full(spdiags([d1 d2 d3], -1:1, N, N));
end

function C = carray(c,N)
    C = zeros(N,N);
    for i =1:N
         C(i,i) = c(i) + c(i+1);
    end
end

function M = marray(m)
    M = diag(m);
end

My goal now is to find the Jacobian matrix for any number of masses. @Torsten provided this answer for a two mass network, where z1,z2 are the state functions of displacement and velocity of mass 1 respectively, and z3,z4 are the state functions of displacement and velocity of mass 2 respectively.

syms z1 z2 z3 z4 m1 m2 c1 c2 c3 k1 k2 k3 gamma1 gamma2 gamma3 F
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;

From here I'm not sure how to find the Jacobian. Below is my (wrong) attempt:

% initial guess
z0_s = [z1,z2,z3,z4]; % Symbolic form
z0_d =[0,0,0,0]; % Double form
% Jacobian
J=jacobian([z1, z2, z3, z4],z0_s); % Symbolic Jacobian matrix
J=subs(J,[z1, z2, z3, z4],z0_d); % Double Jacobian matrix

How do I find the Jacobian for my two mass system? Then, how do I find the Jacobian for any number of masses by leveraging my karray(), marray(), carray(), and farray() functions?

[Sam Chak]

Hi @Jonathan Frutschy

I don't think that @Torsten would provide equations like "z1 = z2 ..." directly. It's possible that you may have interpreted them incorrectly in your MATLAB code. However, after making slight adjustments to the code, it produces the essentially same result as the linear state matrix I derived intuitively in my answer above (although the sequence is different).

syms z1 z2 z3 z4 m1 m2 c1 c2 c3 k1 k2 k3 gamma1 gamma2 gamma3 F

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

z   = [z1, z2, z3, z4];     % state vector
z0  = [0, 0, 0, 0];         % equilibrium point (NOT initial guess)

% Jacobian
J   = jacobian([Eq1; Eq2; Eq3; Eq4], z)  % Symbolic Jacobian matrix

J = 

J   = subs(J, [z1, z2, z3, z4], z0)       % Double Jacobian matrix

J = 

[Sam Chak]

@Jonathan Frutschy, I tested below it appears that there might be an error in your code for creating the stiffness matrix when dealing with a two-mass network. I seldom use the commands circshift(), spdiags(), and full(). Would you be able to rectify the 'make_diagonal()' section of the code?

gamma = 1;
N   = 2;
m   = ones(1, N);
c   = ones(1, N+1);
k   = ones(1, N+1);

%% Test
M   = marray(m)

M = 2x2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C   = carray(c, N)

C = 2x2

2 0 0 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x   = [0; 0];
K   = make_diagonal(x, k, gamma, N)

K = 2x2

-1 -1 -1 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function M = marray(m)
    M = diag(m);
end

function C = carray(c,N)
    C = zeros(N);
    for i = 1:N
         C(i,i) = c(i) + c(i+1);
    end
end

function K = karray(k, gamma, N)
    K = @(x) make_diagonal(x, k, gamma, N);
end

function out = make_diagonal(x, k, gamma, N)
    x   = x(:);
    x   = [x(1:N); 0];
    ck  = circshift(k, -1);
    cg  = circshift(gamma, -1);
    cx  = circshift(x, -1);
    ccx = circshift(x, 1);
    d1  = -3.*ck.*cg.*cx.^2 - ck;
    d2  = (k.*gamma + ck.*cg).*x.^2 + k + ck;
    d3  = -3.*k.*ccx.^ 2 - k;
    out = full(spdiags([d1 d2 d3], -1:1, N, N));
end

[Torsten]

@Sam Chak

I rewrote @Jonathan Frutschy 's system to first-order form to be used with ode45:

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

instead of

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;

I think @Jonathan Frutschy didn't interprete the ' as time derivative.

[Sam Chak]

@Torsten, thank you for providing clarification. I'm not entirely certain, but it's possible that some individuals may mistake z' as the transpose of vector z in MATLAB code, particularly if they are unfamiliar with the mathematical notations used in system dynamics. In the specific case you mentioned, it seems that the OP also misunderstood both and as the displacements, when was described as the time derivative of in the system of 1st-order derivatives that you rewrote.

OP: My goal now is to find the Jacobian matrix for any number of masses.

@Jonathan Frutschy, I'm curious as to why you require the Jacobian matrix when the ode45 solver is capable of solving the problem? The Jacobian matrix is typically employed for efficiently solving highly stiff systems. However, your multiple mass-damper-spring system does not appear to be stiff.

By the way, would you be able to attempt modeling a 4-mass network in MATLAB and share the results here? This experience should assist you in finding a generalized algorithm to compute an n-DOF stiffness matrix .

[Sam Chak]

@Jonathan Frutschy, I think you can find the Jacobian Matrix like this:

N   = 4;

m   = ones(1, N);
M   = marray(m)

M = 4x4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

c   = ones(1, N+1);
C   = carray(c, N)

C = 4x4

2 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

k   = 2;
K   = full(spdiags([-k 2*k -k], -1:1, N, N))

K = 4x4

4 -2 0 0 -2 4 -2 0 0 -2 4 -2 0 0 -2 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Jacobian Matrix
J   = [zeros(N), eye(N);
       - M\K,    - M\C]

J = 8x8

0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 -4 2 0 0 -2 0 0 0 2 -4 2 0 0 -2 0 0 0 2 -4 2 0 0 -2 0 0 0 2 -4 0 0 0 -2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function M = marray(m)
    M = diag(m);
end

function C = carray(c,N)
    C = zeros(N);
    for i = 1:N
         C(i,i) = c(i) + c(i+1);
    end
end

[Jonathan Frutschy]

@Sam Chak make_diagonal() shouldn't be called to generate the K matrix. Only karray() should be called. make_diagonal() is a subfunction of karray() That is why you were getting the wrong K matrix for the two mass example. For clarification, state variables z1,z2,z3,z4 are unitless representations of the system dynamics in real time? That is, z1 = state space (SS) representation of mass 1 displacement, z2 = SS representation of mass 1 velocity, z3 = SS representation of mass 2 displacement, z4 = SS representation of mass 2 velocity.

@Torsten Yes. I thought you meant transpose by ', not the time derivative.

[Jonathan Frutschy]

@Sam Chak Reason behind finding the Jacobian: Yes, you are right, I do not need to specify the Jacobian matrix for ode45() since my system is not stiff. I want to find the Lyapunov exponents of my system, which are defined as the eigenvalues of the Jacobian matrix evauated at the systems initial conditions given by the vector x0 = [x1;x2;x3;x4], where x1 = mass 1 displacement, x2 = mass 1 velocity, x3 = mass 2 displacement, x4 = mass 2 velocity. In state space represenation, the vector x0 becomes z0 = [z1;z2;z3;z4]. I'm curious, why do you call z0 the equilibrium point? What do you mean by equilibrium? Isn't z0 just the initial state of the system and not neccesarly its equilibrium state?

[Sam Chak]

Hi @Jonathan Frutschy

OP: why do you call z0 the equilibrium point? What do you mean by equilibrium?

In dynamical systems , an equilibrium point, also known as a steady state or critical point, is a state where the system remains unchanged over time, when unperturbed. I called it an equilibrium point because it is a solution of the equations of motion where the rates of change of the state variables are all zero .

OP: Isn't z0 just the initial state of the system and not neccesarly its equilibrium state?

Initial state values refer to the values of the system's state variables at a specific starting point or time. They are the conditions that determine the state of the system at the beginning of its evolution. The initial values are typically chosen to be different from the equilibrium point so that we can study the transient behavior. But in your case, you chose the equilibrium point as initial values.

[Torsten]

@Sam Chak

Since F is a time-dependent function, an equilibrium point as such seems to be impossible. Maybe a "periodic equilibrium function" ?

[Jonathan Frutschy]

@Sam Chak Thank you for clearing this up.

[Jonathan Frutschy]

@Torsten @Sam Chak "Since F is a time-dependent function, an equilibrium point as such seems to be impossible. Maybe a "periodic equilibrium function? " Torsten brings up a good point.

Since my K matrix varies as a function of mass position (and consequently time), shouldn't I have a seperate Jacobian matrix for every timestep that ode45() outputs?

[Sam Chak]

The equilibrium point of the nonlinear system mentioned in my previous comment refers to the actual solution of when the system is unperturbed, specifically when . However, if the input signal u has a constant portion c, the equilibrium point may be found at solving .

As mentioned by @Torsten, the input signal u is a unbiased periodic time function, resulting in periodic oscillations observed in the simulation. Although the nonlinear system does not reach the equilibrium point in finite time, the equilibrium point still exists according to its definition. This can be likened to an undamped pendulum described by , where there are two equilibrium points, and . However, if the initial angle deviates from an equilibrium point, stable periodic oscillations will be observed.

Regarding the Jacobian matrix , which is dependent on the displacement vector , it can be computed at any given point . In this particular discussion, the Jacobian matrix mentioned is evaluated specifically at the equilibrium point. It's important to note that the forcing function is independent of the state variable x and can be treated similarly to the derivative of a constant function, .

syms x1 x2 x3 x4 m1 m2 c1 c2 c3 k1 k2 k3 gamma1 gamma2 gamma3 F

Eq1 = x3;
Eq2 = x4;
Eq3 = (F-(c1+c2)*x3-(((k1*gamma1+k2*gamma2)*x1^2+k1+k2)*x1+(-k2*gamma2*x1^2-k2)*x2))/m1;
Eq4 = ( -(c2+c3)*x4-((-k2*gamma2*x2^2-k2)*x1+((k2*gamma2+k3*gamma3)*x2^2+k2+k3)*x2))/m2;

% Jacobian matrix for a two mass network
J   = jacobian([Eq1; Eq2; Eq3; Eq4], [x1, x2, x3, x4])  % Symbolic Jacobian matrix

J = 

[Jonathan Frutschy]

%% Jacobian Matrix
J   = [zeros(N), eye(N);
       - M\K,    - M\C]

Ok, so this is the Jacobian matrix of the system for any mass number N evaluated at the system's equilibrium point X = [x1; x_dot1; x2; x_dot2; x_3; x_dot3; x_4; x_dot4; ... x_N; x_dotN]? @Sam Chak

Answer 2

[Sam Chak]

Hi @Jonathan Frutschy

I have conducted a simulation of the motion for 600 masses using the ode45 solver. Please refer to the comments in the code for further details, as I was unable to utilize your karray() function handle.

N       = 600;
m       = 1*ones(N, 1);
c       = 1*ones(N+1,1);
k       = 2;    % <-- this one must be a scalar, or else K will be incorrect
gamma   = 1;

fun     = @(t, X) odefun(t, X, N, marray(m), carray(c, N), make_diagonal(X, k, gamma, N));
tspan   = [0 60*pi];
X0      = zeros(2*N, 1);
[t, X]  = ode45(fun, tspan, X0); 
plot(t/pi, X), grid on, xlabel('t/\pi'), ylabel('\bf{x}\rm(t)'), title('System responses')

%% Multi-body Mass-Damper-Spring System
function dX = odefun(t, X, N, M, C, K)
    %% Definitions
    x       = X(1:N);               % position state vector
    dx      = X(N+1:2*N);           % velocity state vector
    
    %% Force vector
    alpha   = 0.2;                  % angular frequency of sin
    beta    = 0.15;                 % angular frequency of cos
    famp    = 0.007;                % amplitude of sinusoidal wave
    f       = famp*sin(alpha*t)*cos(beta*t);
    F       = [f; zeros(175,1); f; zeros(100,1); f; zeros(74,1); f; zeros(175,1); f; zeros(71,1)];
    
    %% Equations of Motion
    ddx     = M\(F - K*x - C*dx);
    
    %% State-space model
    dX      = [dx; ...
               ddx];
end

%% Mass matrix, M
function M = marray(m)
    M = diag(m);
end

%% Damping matrix, C
function C = carray(c,N)
    C = zeros(N);
    for i = 1:N
         C(i,i) = c(i) + c(i+1);
    end
end

%% function handle for Stiffness matrix (but unused)
function K = karray(k, gamma, N)
    K = @(x) make_diagonal(x, k, gamma, N);
end

%% The true Stiffness matrix, K (but k must be a scalar, which is a limitation)
function out = make_diagonal(x, k, gamma, N)
    x   = x(:);
    x   = [x(1:N); 0];
    ck  = circshift(k, -1);
    cg  = circshift(gamma, -1);
    cx  = circshift(x, -1);
    ccx = circshift(x, 1);
    d1  = -3.*ck.*cg.*cx.^2 - ck;
    d2  = (k.*gamma + ck.*cg).*x.^2 + k + ck;
    d3  = -3.*k.*ccx.^ 2 - k;
    out = full(spdiags([d1 d2 d3], -1:1, N, N));
end

[Jonathan Frutschy]

@Sam Chak This is great. I'm not sure if I was clear on how karray() works. make_diagonal() is never intended to be called. You only need to call karray() once. The only reason that my system is nonlinear is because the K matrix varies as a function of the mass displacements and velocities. The output of karray() is a function handle of @(X), that way it can be passed to odefun(), which is a function handle of @(X,t). Lets go back to the four mass example for finding the Jacobian. See below:

N = 4;
k = 2;
gamma = 1;

K_X = karray(k,gamma,N); % Function handle @(X) version of K for use in odefun()
K = K_X([zeros(N,1);zeros(N,1)]) % Find K matrix at x1(t), x2(t) = 0, x1_dot(t) = 0, and x2_dot(t) = 0

K = 4x4

4 -2 0 0 -2 4 -2 0 0 -2 4 -2 0 0 -2 4

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% function handle for Stiffness matrix (but unused)
function K = karray(k, gamma, N)
    K = @(x) make_diagonal(x, k, gamma, N);
end
%% The true Stiffness matrix, K (but k must be a scalar, which is a limitation)
function out = make_diagonal(x, k, gamma, N)
    x   = x(:);
    x   = [x(1:N); 0];
    ck  = circshift(k, -1);
    cg  = circshift(gamma, -1);
    cx  = circshift(x, -1);
    ccx = circshift(x, 1);
    d1  = -3.*ck.*cg.*cx.^2 - ck;
    d2  = (k.*gamma + ck.*cg).*x.^2 + k + ck;
    d3  = -3.*k.*ccx.^ 2 - k;
    out = full(spdiags([d1 d2 d3], -1:1, N, N));
end

As you can see the output for the K matrix is the same as what you got. No need to call make_diagonal(), just call karray() to get the K matrix. Does that make sense, or do you still think that my karray() function is wrong?

[Jonathan Frutschy]

%% Jacobian Matrix
J   = [zeros(N), eye(N);
       - M\K,    - M\C]

@Sam Chak How is the forcing function F(t) accounted for in this Jacobian matrix?