ODE 45 to solve coupled ODE

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);

[Torsten]

Show us ftotal.

[aakash dewangan]

Torsten, ftotal expression in command window looks like this

@(t,Y)[Y(2);(sin(pi.*t.*(1.2e1./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(1).*6.579736267392906e-3-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(1).*4.0e2-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(1).*6.579736267392906e-3-Y(1).*5.921762640653616e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.2e1./5.0)).*Y(3).*3.999983550659332e2-sin(pi.*t.*(1.2e1./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(5).*3.999851955933984e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4);Y(4);(sin(pi.*t.*(6.0./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(3).*4.0e2-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(3).*1.644934066848227e-3-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(3).*1.644934066848227e-3-Y(3).*1.480440660163404e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.2e1./5.0)).*Y(1).*3.999934202637326e2-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(5).*3.999851955933984e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4);Y(6);(sin(pi.*t.*(1.8e1./5.0)).*8.0e1-sin(pi.*t.*(6.0./5.0)).^2.*Y(5).*1.480440660163404e-2-sin(pi.*t.*(1.2e1./5.0)).^2.*Y(5).*1.480440660163404e-2-sin(pi.*t.*(1.8e1./5.0)).^2.*Y(5).*4.0e2-Y(5).*1.332396594147063e3-sin(pi.*t.*(6.0./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(3).*3.999983550659332e2-sin(pi.*t.*(1.2e1./5.0)).*sin(pi.*t.*(1.8e1./5.0)).*Y(1).*3.999934202637326e2)./(sin(pi.*t.*(6.0./5.0)).^2+sin(pi.*t.*(1.2e1./5.0)).^2+sin(pi.*t.*(1.8e1./5.0)).^2+9.0e4)]

[Jan]

Which program has created a function including "1.8e1 / 5.0" and (6.0 / 5.0)?!

[aakash dewangan]

That was the expression of ftotal, when i printed it..

Answer 1

[Torsten]

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)