Solving Linear differential Equation
2 次查看(过去 30 天)
显示 更早的评论
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 个)
Star Strider
2012-8-13
编辑:Star Strider
2012-8-13
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 个评论
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!
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!