How to linearize a nonlinear system and find the Jacobian matrix?

Hi, I am working with a highly nonlinear system and would like to simulate its behavior around the equilibrium point in Simulink. To achieve this, I employed a symbolic method to linearize the nonlinear system in MATLAB first and obtain its Jacobian matrix at the equilibrium point. Below is the code I have developed for this purpose.

clear;
clc;

syms q_0 q_1 q_2 q_3 q_dot0 q_dot1 q_dot2 q_dot3 q_ddot0 q_ddot1 q_ddot2 q_ddot3;

q_ddot0 = ((26487*sin(2*q_2))/2 - (214839*sin(2*q_1))/2 - (2943*sin(2*q_3))/2 + (26487*sin(2*q_1 - 2*q_2 + 2*q_3))/2 + (26487*sin(2*q_1 + 2*q_2 - 2*q_3))/2 - 405000*cos(2*q_1 - 2*q_2) + 45000*cos(2*q_1 - 2*q_3) - 225000*cos(2*q_2 - 2*q_3) + 18700*q_dot1^2*sin(q_1) + 3150*q_dot2^2*sin(q_2) + 50*q_dot3^2*sin(q_3) - 2250*q_dot1^2*sin(q_1 - 2*q_2 + 2*q_3) - 2250*q_dot1^2*sin(q_1 + 2*q_2 - 2*q_3) - 450*q_dot2^2*sin(2*q_1 + q_2 - 2*q_3) + 1350*q_dot3^2*sin(2*q_1 - 2*q_2 + q_3) + 1350*q_dot1^2*sin(q_1 - 2*q_2) - 150*q_dot1^2*sin(q_1 - 2*q_3) + 300*q_dot2^2*sin(q_2 - 2*q_3) + 5400*q_dot2^2*sin(2*q_1 - q_2) + 300*q_dot3^2*sin(2*q_1 - q_3) + 900*q_dot3^2*sin(2*q_2 - q_3) + 845000)/(50*(27*cos(2*q_2) - 219*cos(2*q_1) - 3*cos(2*q_3) + 27*cos(2*q_1 - 2*q_2 + 2*q_3) + 27*cos(2*q_1 + 2*q_2 - 2*q_3) - 243*cos(2*q_1 - 2*q_2) + 27*cos(2*q_1 - 2*q_3) - 171*cos(2*q_2 - 2*q_3) + 662));
q_ddot1 = -(3*(15000*cos(q_1 - 2*q_3) - 135000*cos(q_1 - 2*q_2) - (79461*sin(q_1 - 2*q_2))/2 + (8829*sin(q_1 - 2*q_3))/2 + 320000*cos(q_1) - 45000*cos(q_1 - 2*q_2 + 2*q_3) - 45000*cos(q_1 + 2*q_2 - 2*q_3) - 165789*sin(q_1) + (44145*sin(q_1 - 2*q_2 + 2*q_3))/2 + (44145*sin(q_1 + 2*q_2 - 2*q_3))/2 - 150*q_dot2^2*sin(q_1 - q_2 + 2*q_3) + 450*q_dot3^2*sin(q_1 + 2*q_2 - q_3) + 8550*q_dot2^2*sin(q_1 - q_2) + 850*q_dot3^2*sin(q_1 - q_3) + 3650*q_dot1^2*sin(2*q_1) - 450*q_dot1^2*sin(2*q_1 - 2*q_2 + 2*q_3) - 450*q_dot1^2*sin(2*q_1 + 2*q_2 - 2*q_3) + 4050*q_dot1^2*sin(2*q_1 - 2*q_2) - 450*q_dot1^2*sin(2*q_1 - 2*q_3) - 600*q_dot2^2*sin(q_1 + q_2 - 2*q_3) + 1800*q_dot3^2*sin(q_1 - 2*q_2 + q_3) + 1800*q_dot2^2*sin(q_1 + q_2) + 100*q_dot3^2*sin(q_1 + q_3)))/(50*(27*cos(2*q_2) - 219*cos(2*q_1) - 3*cos(2*q_3) + 27*cos(2*q_1 - 2*q_2 + 2*q_3) + 27*cos(2*q_1 + 2*q_2 - 2*q_3) - 243*cos(2*q_1 - 2*q_2) + 27*cos(2*q_1 - 2*q_3) - 171*cos(2*q_2 - 2*q_3) + 662));
q_ddot2 = (3*(30000*cos(q_2 - 2*q_3) + 5886*sin(q_2 - 2*q_3) + 270000*cos(2*q_1 - q_2) - 79461*sin(2*q_1 - q_2) - 225000*cos(q_2) - 45000*cos(2*q_1 + q_2 - 2*q_3) + (97119*sin(q_2))/2 + (26487*sin(2*q_1 + q_2 - 2*q_3))/2 - 150*q_dot1^2*sin(q_1 - q_2 + 2*q_3) + 1050*q_dot3^2*sin(2*q_1 + q_2 - q_3) - 1050*q_dot3^2*sin(2*q_1 - q_2 + q_3) + 13950*q_dot1^2*sin(q_1 - q_2) - 5900*q_dot3^2*sin(q_2 - q_3) + 450*q_dot2^2*sin(2*q_2) - 450*q_dot2^2*sin(2*q_1 - 2*q_2 + 2*q_3) + 450*q_dot2^2*sin(2*q_1 + 2*q_2 - 2*q_3) + 4050*q_dot2^2*sin(2*q_1 - 2*q_2) - 2850*q_dot2^2*sin(2*q_2 - 2*q_3) - 2100*q_dot1^2*sin(q_1 + q_2 - 2*q_3) - 1350*q_dot3^2*sin(q_2 - 2*q_1 + q_3) + 900*q_dot1^2*sin(q_1 + q_2) + 150*q_dot3^2*sin(q_2 + q_3)))/(50*(27*cos(2*q_2) - 219*cos(2*q_1) - 3*cos(2*q_3) + 27*cos(2*q_1 - 2*q_2 + 2*q_3) + 27*cos(2*q_1 + 2*q_2 - 2*q_3) - 243*cos(2*q_1 - 2*q_2) + 27*cos(2*q_1 - 2*q_3) - 171*cos(2*q_2 - 2*q_3) + 662));
q_ddot3 = (3*(60000*cos(2*q_1 - q_3) + 90000*cos(2*q_2 - q_3) - 17658*sin(2*q_1 - q_3) - 17658*sin(2*q_2 - q_3) - 25000*cos(q_3) - 135000*cos(2*q_1 - 2*q_2 + q_3) + (6867*sin(q_3))/2 + (79461*sin(2*q_1 - 2*q_2 + q_3))/2 - 450*q_dot1^2*sin(q_1 + 2*q_2 - q_3) - 1950*q_dot2^2*sin(2*q_1 + q_2 - q_3) + 1950*q_dot2^2*sin(2*q_1 - q_2 + q_3) + 2350*q_dot1^2*sin(q_1 - q_3) + 13100*q_dot2^2*sin(q_2 - q_3) - 50*q_dot3^2*sin(2*q_3) + 450*q_dot3^2*sin(2*q_1 - 2*q_2 + 2*q_3) - 450*q_dot3^2*sin(2*q_1 + 2*q_2 - 2*q_3) - 450*q_dot3^2*sin(2*q_1 - 2*q_3) + 2850*q_dot3^2*sin(2*q_2 - 2*q_3) - 6300*q_dot1^2*sin(q_1 - 2*q_2 + q_3) + 1350*q_dot2^2*sin(q_2 - 2*q_1 + q_3) + 200*q_dot1^2*sin(q_1 + q_3) - 150*q_dot2^2*sin(q_2 + q_3)))/(50*(27*cos(2*q_2) - 219*cos(2*q_1) - 3*cos(2*q_3) + 27*cos(2*q_1 - 2*q_2 + 2*q_3) + 27*cos(2*q_1 + 2*q_2 - 2*q_3) - 243*cos(2*q_1 - 2*q_2) + 27*cos(2*q_1 - 2*q_3) - 171*cos(2*q_2 - 2*q_3) + 662));
eqn     = [q_ddot0,q_ddot1,q_ddot2,q_ddot3];

J    = jacobian(eqn, [q_0 ,q_1 ,q_2 ,q_3 ,q_dot0, q_dot1 ,q_dot2 ,q_dot3]);

% Solve the equation to find equilibrium point
sol  = solve(eqn, [q_0 ,q_1 ,q_2 ,q_3 ,q_dot0, q_dot1 ,q_dot2 ,q_dot3]);

q_00    = sol.q_0;
q_10    = sol.q_1;
q_20    = sol.q_2;
q_30    = sol.q_3;
q_dot00 = sol.q_dot0;
q_dot10 = sol.q_dot1;
q_dot20 = sol.q_dot2;
q_dot30 = sol.q_dot3;

A = double(subs(J, [q_0, q_1, q_2, q_3, q_dot0, q_dot1, q_dot2, q_dot3], [q_00, q_10, q_20, q_30, q_dot00, q_dot10, q_dot20, q_dot30]));
disp('Linearized system matrix A:');
disp(A);

[Sam Chak]

Hi @TADDESE DIRIBA

Trigonometric functions, like sine and cosine, in dynamical systems are used to model the periodic behavior. Depending on how these functions are incorporated into the system equations, it is possible to have one or more equilibrium points or none at all.