ODE to state space conversion

조회 수: 43 (최근 30일)
aakash dewangan
aakash dewangan 2021년 12월 8일
편집: 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

채택된 답변

aakash dewangan
aakash dewangan 2023년 10월 25일
편집: aakash dewangan 2023년 10월 25일
'odeToVectorField' command does the job

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Function Creation에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by