ODE to state space conversion
16 views (last 30 days)
Show older comments
aakash dewangan
on 8 Dec 2021
Edited: aakash dewangan
on 25 Oct 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
0 Comments
Accepted Answer
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!