Solving Linear differential Equation

Question

Jerry 2012-8-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/45850-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)

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jerry 2012-8-12

Thanks Walter for editing.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2012-8-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/45850-solving-linear-differential-equation#answer_56065

编辑：Star Strider 2012-8-13

在 MATLAB Online 中打开

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  ====================

8 个评论
显示 6更早的评论隐藏 6更早的评论

Jerry 2012-8-13

在 MATLAB Online 中打开

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.

Star Strider 2012-8-13

编辑：Star Strider 2012-8-13

在 MATLAB Online 中打开

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 2012-8-16

thanks :)

Star Strider 2012-8-16

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

Again, my pleasure!

请先登录，再进行评论。

Solving Linear differential Equation

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Solving Linear differential Equation

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论