MATLAB: Use ODE45 to solve a system of two coupled second order ODEs

ode45second order ode

I have the following 2nd order differential equations I need to solve:
x1''=(F(t)-(c1+c2)*x1'+c2*x2'-(k1+k2)*x1+k2*x2)/m1
x2''=(c2*x1'-c2*x2'+k2*x1-k2*x2)/m2
All k, c, m and F(t) are known. Derivatives are wrt time. I tried to lay it out as if it was a single 2nd order ode.
function [ xp ] = ftotal( t,x )
xp=zeros(4,1)
xp(1)=x(2);  %x1 dot
xp(2)=(F-(c1+c2)*x(2)+c2*x(4)-(k1+k2)*x(1)+k2*x(4))/m1;
xp(3)=x(5);
xp(4)=(c2*x(2)-c2*x(5)+k2*x(1)-k2*x(4))/m2;
end
[t,x]=ode45('ftotal',[0,1],[0 0])
As expected, it gave "Index exceeds matrix dimensions". All my initial conditions are 0, I know there should be 4 (x1, x1', x2, x2'=0) but I'm not sure how to input them. I have little experience with ODE45. My goal is to obtain x1 and x1' vs t, and x2 and x2' vs t. Any help is appreciated.

Best Answer

I prefer to let the Symbolic Math Toolbox do these derivations:
syms c1 c2 k1 k2 m1 m2 t Ftfcn x1(t) x2(t) Y 
%x1''=(F(t)-(c1+c2)*x1'+c2*x2'-(k1+k2)*x1+k2*x2)/m1
Dx1 = diff(x1);
D2x1 = diff(x1,2);
Dx2 = diff(x2);
D2x2 = diff(x2,2);
Eq1 = D2x1 == (Ftfcn-(c1+c2)*Dx1+c2*Dx2-(k1+k2)*x1+k2*x2)/m1
% x2''=(c2*x1'-c2*x2'+k2*x1-k2*x2)/m2
Eq2 = D2x2 == (c2*Dx1-c2*Dx2+k2*x1-k2*x2)/m2
[VF,Subs] = odeToVectorField(Eq1, Eq2)
ftotal = matlabFunction(VF, 'Vars',{t,Y,Ftfcn,c1,c2,k1,k2,m1,m2})
This produces (lightly edited):
ftotal = @(t,Y,Ftfcn,c1,c2,k1,k2,m1,m2)[Y(2);-(c2.*Y(2)-c2.*Y(4)+k2.*Y(1)-k2.*Y(3))./m2;Y(4);(Ftfcn(t)-(c1+c2).*Y(4)-(k1+k2).*Y(3)+c2.*Y(2)+k2.*Y(1))./m1];
Substituting random values and a random function:
Ftfcn = @(t) sin(t);
c1 = rand;
c2 = rand;
k1 = rand;
k2 = rand;
m1 = rand;
m2 = rand;
ic = zeros(4,1);
tspan = [0,1];
[T,Y] = ode45(@(t,Y) ftotal(t,Y,Ftfcn,c1,c2,k1,k2,m1,m2), tspan, ic);
plot(T, Y)
grid
This runs.
Experiment to get the result you want.
Related Question