Solving a system of second order differential equations (coupled mass-spring system)

19 views (last 30 days)
Hello, I am trying to solve the system defined by the following equations EqX1 and EqX2:
syms X1(t) X2(t) K f m1 m2 Y t
EqX1 = diff(X1,2) == K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f + K*(X2-X1)/m2;
[ODE,Vars] = odeToVectorField(EqX1, EqX2);
Sys = matlabFunction(ODE, 'Vars', {t, Y, K, f, m1, m2});
With arbitrary constants m1, m2, f, and K.
m1 = 1; % Mass 1
m2 = 2; % Mass 2
K = 100; % Spring constant
f = 5; % Forcing term on mass 2
All of my initial conditions are zero, except for the velocity of mass 2 (v2_i).
v2_i = 3;
ic = [0 v2_i 0 0];
tspan = [0, 1];
[T,Y] = ode45(@(t,Y) Sys, tspan, ic);
I feel like this should work, but when I try to run it, I get the following errors:
Error using odearguments (line 95)
@(T,Y)SYS returns a vector of length 1, but the length of initial conditions vector is 4. The vector returned by @(T,Y)SYS and the initial conditions vector must have the same
number of elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in CoupledMassSim (line 29)
[T,Y] = ode45(@(t,Y) Sys, tspan, ic);
Having looked at similar problems on this forum, I'm not sure why @(T,Y)SYS is returning a vector of length 1. Here is the Sys function for reference:
Sys =
function_handle with value:
@(t,Y,K,f,m1,m2)[Y(2);f+(K.*(Y(1)-Y(3)))./m2;Y(4);-(K.*(Y(1)-Y(3)))./m1]

Accepted Answer

Stephan
Stephan on 16 Aug 2019
syms X1(t) X2(t) K f m1 m2 Y t
EqX1 = diff(X1,2) == K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f + K*(X2-X1)/m2;
[ODE,Vars] = odeToVectorField(EqX1, EqX2)
Sys = matlabFunction(ODE, 'Vars', {t, Y, K, f, m1, m2})
m1 = 1; % Mass 1
m2 = 2; % Mass 2
K = 100; % Spring constant
f = 5; % Forcing term on mass 2
v2_i = 3;
ic = [0 v2_i 0 0]';
tspan = [0, 1];
[T,Y] = ode45(@(t,Y)Sys(t,Y,K,f,m1,m2), tspan, ic);
plot(T,Y)
  1 Comment
Benjamin Tonita
Benjamin Tonita on 16 Aug 2019
Thank you!
Also turns out I had my signs turned around. Equations should be:
EqX1 = diff(X1,2) == -1*K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f - K*(X2-X1)/m2;

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!