Solving a forced mass-spring-damper system with Runge Kutta method in matlab

Question

I am trying to solve a forced mass-spring-damper system in matlab by using the Runge-Kutta method. Currently the code uses constant values for system input but instead I would like to vectors as input. For examples, I would like to replace my force amplitude F0 with a vector value.

Should I be using for loops or what is the simplest way to do it?

function O = MSDSRK(m,b,k,F0,w,x0,v0) 

% ----- Input argument -----
% m: mass for particle 
% b: damping coefficient 
% k: spring constant 
% F0: amplitude of external force 
% w: angular freuency of external force 
% x0: initial condition for the position x(0)
% v0: initial condition for the velocity v(0)

dt=0.1;

options=odeset('InitialStep',dt,'MaxStep',dt);

td=[0:dt:50];

% Solve differential equation with Runge-Kutta solver 

[t,x]=ode45(@(t,X)MSD(t,X,m,b,k,F0,w),td,[x0;v0],options);

% Extract only particle position trajectory

O=[t x(:,1)];

end


function dX=MSD(t,X,m,b,k,F0,w)

% With two 1st order diffeential equations,
% obtain the derivative of each variables
% (position and velocity of the particle)

dX(1,1)=X(2,1); 

dX(2,1)=(1/m)*(F0*sin(w*t)-b*X(2,1)-k*X(1,1)); 
end

I would use loop in this case, since the output of your function is already a vector. If you really want to vectorize it, make you function to produce output in the form of matrix, one of dimensions of which corresponds to different values of F0. — freude

Rody Oldenhuis Rody Oldenhuis · Accepted Answer · 2014-02-05T08:06:16

A for-loop will work for sure, but it's not what I'd advise in this situation. That approach has the drawback of making you think in a certain way, namely, that you are solving a simple 1D system, just multiple times.

The for-loop approach doesn't teach you very much more than what you already knew how to do. You'll get a passing grade I'm sure. But if you want to excell, not only here but also in future classes (and more importantly, in your job later), you have to do more. You can only learn how to tackle more complicated problems by thinking about this simple problem in a different way, namely, as a multi-dimensional mass/spring/damper system, of which the components all happen to be uncoupled and all happen to have the same initial values.

For such systems, vectorization and matrix manipulation is really the way to go. That is the generalization of 1D/2D/3D systems, to arbitrary-D systems. And matrix manipulation and vecorization, that is just what MATLAB is best at. That is indeed what you can learn here.

Here's how to do it for your case. I made it a bit smaller in scope, but the principles remain the same:

function O = MSDSRK

    % Assign some random scalar values to all parameters, but 
    % a *vector* to the force amplitudes 
    [m,b,k,F0,w,x0,v0] = deal(...
        1,0.2,3,  rand(13,1),  2,0,0);

    % We need to simulate as many mass-spring-damper systems as there are
    % force amplites, so replicate the initial values
    x0 = repmat(x0, numel(F0),1);
    v0 = repmat(v0, numel(F0),1);

    % The rest is the same as before
    dt = 0.1;
    options = odeset('InitialStep', dt,'MaxStep', dt);

    td = 0:dt:50;

    [t,x] = ode45(@(t,X) MSD(t,X,m,b,k,F0,w), td, [x0;v0], options);

    % The output is now: 
    %
    % x(:,1)  position of mass in system with forcing term F0(1)
    % x(:,2)  position of mass in system with forcing term F0(2)
    % ...
    % x(:,14) speed of mass in system  with forcing term F0(1)
    % x(:,15) speed of mass in system  with forcing term F0(2)
    % ...    
    O = [t x(:,1:numel(F0))];

    % And, to make the differences more clear:
    plot(t, x(:,1:numel(F0)))

end


function dX = MSD(t,X,m,b,k,F0,w)

    % Split input up in position and velocity components
    x = X(1:numel(X)/2);
    v = X(numel(X)/2+1:end);

    % Compute derivative
    dX = [
        v
        (F0*sin(w*t) - b*v - k*x)/m
    ];
end

Solving a forced mass-spring-damper system with Runge Kutta method in matlab

1 Answers