Solving Linear differential Equation

I have an array

X=[ 0.1000    0.0844    0.0434   -0.0090   -0.0559   -0.0831   -0.0832   -0.0574   -0.0152    0.029]

Which shows the position and time array with respect to that position is

T= [0    0.5000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.5000]

And I have the mass spring equation mx’’ + c x’ + kx = 0 where x’’ is the double derivative of x which I have found by using dx=diff(x.2) and dt2=diff(t,2) and x’ is found by dx=diff(x) and dt=(diff). Problem is I have implemented the code to find the value of c and k in the equation using A=x\b formula

I have impleneted the code by using

xx=dx./dt;   xx2=dx2./dt2;

The values obtained using the formula A=x\b are Nan and Nan both for c and k because my dt2=diff(t,2) comes out to be zero and I have even pad zeros to make the size equal for xx and xx2 but what can I do to make the size equal apart from padding zeros since I thnk padding zeros is causing a lot of issue :( is there a way I can like interpolate and get the sizes equal for the diff since diff is reducing the size by n-1, correct and what can be done about dt2 is it fine or it should be dt2=dt^.2 since its coming out to be all zero.

Below is my code.

    x=[ 0.1000    0.0844    0.0434   -0.0090   -0.0559   -0.0831   -0.0832   -0.0574   -0.0152    0.029]';
    t= [0    0.5000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.5000]';  
    dx=diff(x);
    dt=diff(t);
    dx2=diff(x,2);
    dt2=diff(t,2); % this comes out zero
    xx=dx./dt;

    xx2=dx2./dt2;
    % padding zeros to make size equal

    xx2=padarray(xx2,size(x)-size(xx2),'post');
    xx=padarray(xx,size(x)-size(xx),'post');
    mass=100;
    gh=horzcat(xx,x);
    A=gh\(m*xx2)

[Jerry]

Thanks Walter for editing.

Answer 1

[Star Strider]

I got an excellent fit with the method I suggested earlier. I just discovered you listed mass m = 100, so I used it to estimate these parameters (output copied from the Command Window):

  Coefficient of viscous friction, Kd =   1.058501468043754E+01
  Spring Constant, Ks =   1.305291301729188E+02
  Initial X(1) =  -7.899752805718143E-03
  Initial X(2) =   1.001228985786043E-01

The fit I got to your X vector with these parameters is:

          X            Estimate
     100.0000e-003   100.1229e-003
      84.4000e-003    84.3690e-003
      43.4000e-003    43.3323e-003
      -9.0000e-003    -8.9916e-003
     -55.9000e-003   -55.8367e-003
     -83.1000e-003   -82.9966e-003
     -83.2000e-003   -83.0860e-003
     -57.4000e-003   -57.4732e-003
     -15.2000e-003   -15.4055e-003
      29.0000e-003    29.2575e-003

Did the objective function that integrated your differential equation to estimate the parameters with nlinfit or lsqcurvefit work in your application? (Sometimes you have to try a few times with different initial parameter guesses to get a good fit. I used the objective function I sent you earlier with good results in my own work, as well as in estimating the parameters from the data you posted.)

If you want the both X(1) and X(2), it's fairly easy to get them once you have your estimated parameters. You have to add an extra tilde (~) argument to the end of the argument list of the objective function, and add a nargin ‘if’ block. The function will return both X vectors when you want them, and will still work with nlinfit or lsqcurvefit when you only want X(2).

I had to transpose your T and X vectors to column vectors to work with my functions, since I use the column-major convention for everything. This was my call to lsqcurvefit:

          OptsSMD = optimset('MaxFunEvals',5000, 'MaxIter',5000, 'TolFun',1E-010, 'TolX',1E-010, 'Algorithm','Levenberg-Marquardt');
          Lobnd = [];
          Hibnd = [];
          B0a = rand(4,1) * 100;
          [Ba,Rnma,Rsda] = lsqcurvefit(@SpringMassDamper,B0a,T,X,Lobnd,Hibnd,OptsSMD, m);
          SimSMDa = SpringMassDamper(Ba,T, m)     % Calculate fitted data for plot

This is the function I used to estimate these values:

function Xv = SpringMassDamper(B, t, m)
% 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 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

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

end

% ====================  END: SpringMassDamper.m  ====================

[Jerry]

I used the function above but i dont get values as you have got. I used the code

x=[ 0.1000 0.0844 0.0434 -0.0090 -0.0559 -0.0831 -0.0832 -0.0574 -0.0152 0.029]'; % Which shows the position and time array with respect to that position is

t= [0    0.5000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.5000]';  

m=100;
OptsSMD = optimset('MaxFunEvals',5000, 'MaxIter',5000, 'TolFun',1E-010, 'TolX',1E-010, 'Algorithm','Levenberg-Marquardt');
          Lobnd = [];
          Hibnd = [];
          B0a = rand(4,1) * 100;
          [Ba,Rnma,Rsda] = lsqcurvefit(@SpringMassDamper,B0a,t,x,Lobnd,Hibnd,OptsSMD, m);
          SimSMDa = SpringMassDamper(x,t, m)     % Calculate fitted data for plot

I have also taken the column vector but how ave to got the values of Kd and Ks, i didnt get any such values but dont you think the equation was mx’’ + c x’ + kx = 0 where x and t is given and there is no purpose for x'' and x'. I completely agree with your code and i am grateful but i arnt able to get the values of Kd and Ks as you did, please see my above code is it correct and how can i obtain kd and Ks as you have obtained.

[Jerry]

these are my outputs

Optimization terminated: the relative change in the sum-of-squares of the functions is less than options.TolFun

Ba =

     10.5850
    130.5291
     -0.0079
      0.1001

Rnma =

    1.6234e-007

Rsda =

    1.0e-003 *

      0.1229
     -0.0310
     -0.0677
      0.0084
      0.0633
      0.1034
      0.1140
     -0.0732
     -0.2055
      0.2575

SimSMDa =

     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090
     -0.0090

[Star Strider]

First, I apologise for the earlier confusion.

I see the problem you're having. The first argument to SpringMassDamper is the parameter vector, in this case Ba, so your call to it using:

SimSMDa = SpringMassDamper(x,t, m)     % Calculate fitted data for plot

isn't going to work. Use mine instead:

SimSMDa = SpringMassDamper(Ba, t, m)   % Calculate fitted data for plot

In my code, Kd is Ba(1) and Ks is Ba(2). (The initial conditions are Ba(3:4)). I used this line:

fprintf(1,'\n\tCoefficient of viscous friction, Kd = %23.15E\n\tSpring Constant, Ks = %23.15E\n\tInitial X(1) = %23.15E\n\tInitial X(2) = %23.15E\n\n', Ba)

to produce the output I posted. (I'm sending it along for your convenience.) I suggest you copy it and paste it as the command directly after your call to lsqcurvefit.

You may have to run the regression a few times to get a good fit to your data. That's one of the characteristics of nonlinear parameter estimation. Also, if it seems to be taking longer than it should to converge, interrupt it with a ‘Control-C’ and start it again. Yours isn't a difficult problem, but the curve-fitting algorithm can occasionally get lost in the wilderness on even straightforward problems.

Here are a couple other code snippets you may find useful:

          tc = [0:0.1:4.5]';                        % Higher-resolution time vector for plot
          SimSMDac = SpringMassDamper(Ba, tc, m);   % Evaluate function with ‘tc’ time vector

          figure(8)
          plot(t, x, 'pb')
          hold on
          plot(t, SimSMDa, '+r', 'MarkerSize',5)
          plot(tc, SimSMDac, '--r')
          hold off
          xlabel('Time (\mus)')
          ylabel('Displacement (parsecs)')
          grid

That should work. I'll keep this open for a while to be sure you don't have problems with my code.

[Jerry]

yes i have got the above results similar to your but may i know why have you taken X0=B(3:4), and what if i change it to a scalar value it doesnt give the right answer? why (3:4)?

I must say your method is awesome

[Jerry]

Values are comong absolutely correct but there is one probelm that the answer obtained are

Coefficient of viscous friction, Kd =   1.058501468043754E+01
  Spring Constant, Ks =   1.305291301729188E+02

thats with +ve power and my actual values through which i have to compute error is

   Kd =   1.058501468043754E-01
   Ks =   1.305291301729188E-02

Actual values are these, what can be done to attain these, i understood why you have taken B(3,4) after signle steping each and every line, i understood that why :) Thanks Star for the help you have given me, i wish i can be of any worth to you.

[Star Strider]

There may be many solutions to the differential equation, all of which are proportional, because of the nature of the eigenvalues of the equation. If you want to see what solutions may exist for the values of ( 0 < Kd < 1 ) and ( 0 < Ks < 1 ), then change the B0a statement to this:

B0a = rand(4,1) / 100;

or even:

B0a = rand(4,1) / 1000;

You may have to run the fit several times until it converges.

It is my pleasure to help you.

[Jerry]

thanks :)

[Star Strider]

My guess is that it converged without a problem on the values with the orders-of-magnitude you want.

Again, my pleasure!