Mass spring system equation help

I am good at Matlab programming but over here I am stuck in the maths of the problem, I am dealing with the differential equation of spring mass system mx’’+cx’+kx=0 where x’’=dx2/dt2 and x’=dx/dt. I have the values of mass and I also have the array of time and x i.e x is given for a particular value of time so I can find x’’ and x ‘ easily. I am stuck at what method to apply to find the value of c and k. I can program any method but have searched several books but didn’t get how to find c and k. If I get to know the method I can program it easily.

Answer 1

[Star Strider]

I suggest you write an objective function containing your differential equation that would integrate the differential equation, then do curve-fitting of x(t) with nlinfit or lsqcurvefit.

I have provided an objective function here that may work for you. You do not have to pass the parameters specifically to DifEq. It has access to them in the function space. This code assumes all variables and parameters are column vectors.

function X = SMD(B, t, m)     % ‘SMD’ for ‘Spring-Mass-Damper’
% B = parameter and initial conditions column vector (unless you want to
%       supply the initial conditions in this function rather than passing
%       them as parameters).

X0 = B(3:4);                  

% This gives you the option of passing the last two entries of the
%   parameter vector as the initial values of your ODE.  In this case, the
%   curve-fitting function will also fit the initial conditions as well as
%   Kd and Ks.  If you don't want the initial conditions to be parameters,
%   B becomes [2 x 1] and you define X0 as whatever you like in the
%   function.

[T,X] = ode45(@DifEq, t, X0);

      function xdot = DifEq(t, x)
            %   B(1) is the coefficient of viscous friction (‘damper’), Kd;
            %   B(2) is the spring constant, Ks;
            xdot = zeros(2,1);
            xdot(1) =  x(2);
            xdot(2) = -x(1)*B(2)/m  -x(2)*B(1)/m;
      end

X = xdot(:,2);     % Assumes ‘xdot’ is a [N x 2] matrix where N = length(t)

end

I've had success with this general approach in several situations. You may have to experiment with it a bit to fit your particular situation. See Passing Extra Parameters for details in passing the mass m to SMD, or you can define m in the function if it never changes. If you do, then remove it from the SMD argument list.

[Star Strider]

Thank you for accepting my answer!

[Jerry]

I am grateful to all for the help.

Answer 2

[Sumit Tandon]

I feel that if you can calculate the values of x' and x'', then you could take a couple of sets of x, x' and x'' and get a system of linear equation of the form Ax - B = 0.

After this, you could use the backslash operator (\) or MLDIVIDE in MATLAB to solve for c and k.

Any other ideas?

[Jerry]

i dont thnk this would work.

[Jerry]

I need the method to find the coefficients.

Answer 3

[Jacob]

In response to Sumit above:

x = 2xn matrix of x' and x, A = 1x2 matrix of c and k, and B = 1xn matrix of mx''. Solve for A.

Example:

x=[v1 v2 v3 v4;x1 x2 x3 x4];
B=[mx1'' mx2'' mx3'' mx4''];

A=x\B;

This should get you the values of c and k (c=A(1), k=A(2)).

[Jerry]

yes that can be done but i am looking for a constant value for my system that is a constant number for c and k

[Jacob]

That would return a single value for c and k. A is 1x2 ([c k]), and MATLAB would find the "best fit" of c and k using that method.

[Jerry]

what i didnt understand is how you have used x, meaning i have the equation mx’’+cx’+kx=0 where x’’=dx2/dt2 and x’=dx/dt. and i have the array for x and t so i can find dx/dt and stuff easily but OK LET ME TRY AND THEN I WILL LET YOU KNOW HERE.

[Greg Heath]

1. Divide the ODE by m and treat the normalized parameters c1 = c/m and k1 = k/m as unknowns

2. Convert the ODE and Initial Conditions to a set of Linear Finite Difference Equations

3. As suggested by previous repliers, rearrange the equations into linear matrix equation form with C = [ c1 ; k1 ] as the unknown.

4. Solve using backslash.

[Jerry]

yeah something like this mx'' + cx' + kx = 0

determine roots of the equation mr² + cr + k = 0 mr² + cr + k = 0 r² + (c/m)r + k/m = 0

completing the square r² + (c/m)r + k/m = 0 ( r + 0.5c/m)² - (0.5c/m)² + k/m = 0 ( r + 0.5c/m)² - 0.25c²/m² + k/m = 0 ( r + 0.5c/m)² = 0.25c²/m² - k/m ( r + 0.5c/m)² = (c² - 4mk) / ( 4m² ) ( r + c / 2m) = ±√[ (c² - 4mk) / ( 4m² ) ] r = -(c / 2m ) ± √[ (c² - 4mk) / ( 4m² ) ]

[Jerry]

Now i have got enough information to formulate a code i can now but i hew one little query, if i have a vector x =[0 5.5.-----} and a time vector of same length, is it possible that we can find the derivative of the x vector with respect to the time vector? I just need to get this done, i have the code for rest.

Answer 4

[Greg Heath]

There appears to be 2 straightforward approaches:

1. For c1=c/2m, k1=k/m and sufficiently small dt(i) = t(i)-t(i-1)

   a. Obtain an inhomogeneous system of linear equations for C = [c1 ; k1] by       converting the differential equation to a difference equation in x(i).

   b. Obtain the solution to A*C=B via C = A\B.

2. Use NLINFIT or LSQCURVEFIT to estimate c1, k1, A and B from the form of the exact solution

 x(i) = exp(-c1*t(i)) * ( A * cos( sqrt(c1^2-k1) * t(i)) 

        + B * sin( sqrt(c1^2-k1) * (t(i) )

However, since the equation is linear in A and B, a two step estimation of [c1 ; k1 ] and [ A B] might be useful.

Hope this helps.

Greg