ODE to state space conversion

Question

aakash dewangan 2021-12-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1605960-ode-to-state-space-conversion

编辑： aakash dewangan 2023-10-25

Hello,

I have written the program given below. In this program, I have 3 ODEs. I am converting these ODEs into statespace form using in-build function of MATLAB. When I run it for different values/cases, it changes the substitution "S". Can you please tell me how does it decide the substitution?

This is very important to know in my program to contunue the work.

Thank you :)

clc; clear all; close all 
syms p1(t) p2(t) p3(t) rho L m v T k G             
rho = 1.3;  T = 45000;  L = 60;  m = 1;  v = 400*1000/3600;  k = 10;  G = 0.1
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);    
% Mass matrix terms
AA = rho*L/2 + m*(sin(pi*v*t/L))^2;
BB = m*sin(2*pi*v*t/L)*sin(pi*v*t/L);
CC = m*sin(3*pi*v*t/L)*sin(pi*v*t/L);
DD = rho*L/2 + m*(sin(2*pi*v*t/L))^2;
EE = m*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
FF = rho*L/2 + m*(sin(3*pi*v*t/L))^2;
% Stiffness matrix terms 
GG = T*(pi/L)^2*(L/2) + k*(sin(pi*v*t/L))^2;
HH = k*sin(2*pi*v*t/L)*sin(pi*v*t/L);
II = k*sin(pi*v*t/L)*sin(3*pi*v*t/L);
JJ = T*(2*pi/L)^2*(L/2) + k*(sin(2*pi*v*t/L))^2;
KK = k*sin(2*pi*v*t/L)*sin(3*pi*v*t/L);
LL = T*(3*pi/L)^2*(L/2) + k*(sin(3*pi*v*t/L))^2;
% RHS
MM = k*G*sin(pi*v*t/L);
NN = k*G*sin(2*pi*v*t/L);
OO = k*G*sin(3*pi*v*t/L);
% Equation (coupled system of ODE to solve for p)
Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) +  GG*p1 + HH*p2 + II*p3 == MM; % Equation 1
Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) +  HH*p1 + JJ*p2 + KK*p3 == NN; % Equation 2
Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) +  II*p1 + KK*p2 + LL*p3 == OO; % Equation 3
%% ODE to state space conversion
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);             % converts ODE in state space form
ftotal = matlabFunction(V, 'Vars',{'t','Y'});        % Using readymade MATLAB function to solve using ODE 45
interval = [0 L/v];                                  % Time Interval to run the program 
y0 = [0 0 0 0 0 0];                                  % Initial conditions
pSol = ode45(@(t,Y)ftotal(t,Y),interval,y0);         % Using ODE 45 to solve stste space form of ODE
tValues = linspace(interval(1),interval(2),180);     % deviding time interval    
p2Values = deval(pSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
p1Values = deval(pSol,tValues,3); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
p3Values = deval(pSol,tValues,5); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions