ODE 45 to solve coupled ODE

1 view (last 30 days)
aakash dewangan
aakash dewangan on 28 May 2021
Commented: aakash dewangan on 30 May 2021
Hello,
I have written a code to solve a system of coupled second order differential equations using ODE 45.
I wrote this code by modifying the available code. I want to know that what are the initial conditions being applied to which variable (p1, p2, p3). Please let me know what are the initial conditions for p1, Dp1, p2, Dp2, p3, Dp3. And also which column of "y" in "pSol" corresponds to p1, p2 and p3 (in other words how can i get p1, p2, p3 from pSol).
Thanks in advance,
clc; clear all; close all
syms p1(t) p2(t) p3(t)
rho L m v T k G
% Parameters
rho = 9000; T = 15000; L = 100; m = 5; v = 120; k = 2000; G = 0.2
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);
Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) + GG*p1 + HH*p2 + II*p3 == MM;
Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) + HH*p1 + JJ*p2 + KK*p3 == NN;
Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) + II*p1 + KK*p2 + LL*p3 == OO;
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y'});
interval = [0 L/v];
y0 = [1 2; 3 4; 5 6]; %initial conditions
% v-k2
pSol = ode45( @(t,Y)ftotal(t,Y),interval,y0);
  4 Comments
Jan
Jan on 29 May 2021
Which program has created a function including "1.8e1 / 5.0" and (6.0 / 5.0)?!
aakash dewangan
aakash dewangan on 30 May 2021
That was the expression of ftotal, when i printed it..

Sign in to comment.

Answers (1)

Torsten
Torsten on 29 May 2021
Edited: Torsten on 29 May 2021
To be honest, I'd prefer to know what Matlab solves.
Order the variables as
[z1,z2,z3,z4,z5,z6] = [y1,y2,y3,y1',y2',y3']
Then your adapted mass matrix becomes
function MM = mass(t,z)
%Define M
MM = [eye(3),zeros(3);zeros(3),M];
end
and your right-hand side vector becomes
function dz = fun(t,z)
% Define K and F
dz = [zeros(3),eye(3);-K,zeros(3)]*z + [F;zeros(3,1)]
end
and the call to ode15s
options = odeset('Mass',@mass);
[T,Z] = ode15s(@fun,tspan,z0,options)

Community Treasure Hunt

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

Start Hunting!