ODE to state space conversion

44 visualizaciones (últimos 30 días)
aakash dewangan
aakash dewangan el 8 de Dic. de 2021
Editada: aakash dewangan el 25 de Oct. de 2023
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

Respuesta aceptada

aakash dewangan
aakash dewangan el 25 de Oct. de 2023
Editada: aakash dewangan el 25 de Oct. de 2023
'odeToVectorField' command does the job

Más respuestas (0)

Categorías

Más información sobre Function Creation en Help Center y File Exchange.

Productos


Versión

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by