# ODE 45 to solve coupled ODE

4 views (last 30 days)
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).
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);
aakash dewangan on 30 May 2021
That was the expression of ftotal, when i printed it..

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']
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)