ODE to state space conversion

43 views (last 30 days)
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

Accepted Answer

aakash dewangan
aakash dewangan on 25 Oct 2023
Edited: aakash dewangan on 25 Oct 2023
'odeToVectorField' command does the job

More Answers (0)

Categories

Find more on Function Creation in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!