How can i fit ODE parameters to given frequency and amplitude data pairs with "lsqcurvefit"?

Hello there,

i have a problem to implementing a ODE to the lsqcurvefit function. My ODE: a simple spring-damper-mass oszillator. And i have amplitude (mampl) and frequency (frqSet) data to which the parameters should be adjusted.

I have written this so far but i am not getting anywhere:

xd=frqSet;

yd=mampl;

x0=[1e-4 1e+0 1e-10 1e-9];

opts = optimoptions('lsqcurvefit','MaxFunctionEvaluations',10000,'FunctionTolerance',1e-8,'StepTolerance',1e-8,'MaxIterations',10000,'OptimalityTolerance',1e-8);

lx = [1e-4 1e+0 1e-10 1e-9];

ux = [1e-2 1e+2 1e-7 1e-6];

x = lsqcurvefit(@fitfun,x0,xd,yd,lx,ux,opts);

function fitfun(x,xd)

ts = [2*pi/1.001e+5 2*pi/1.002e+5];

u0 = 1e-7;

[ta,ua] = ode45(@(u,t) f(u,t),ts,u0);

function dudt= f(u,t)

dudt = x(4).*sin(xd.*t)./(x(2))-2.*x(1).*sqrt(x(2)./x(3)).*x(3).*u(1)./x(2)-u(2).*x(3)./x(2);

end

end

Thanks for help in advance.

[Matt J]

And i have amplitude (mampl) and frequency (frqSet) data to which the parameters should be adjusted.

What are the unknown parameters and what do mampl and frqSet correspond to in terms of the mathematical notation in the ODE,

that you have shown us?

[Paul-Adam]

The parameters what i am looking for are: resp. . With focus on . mampl are the Amplitudes of for different Frequencies (frqSet).

[Bjorn Gustavsson]

In classical mechanics at my school we used to call this mampl-as-function-of frequency "gain". The OP has a curve g(f) for this type of driven damped mass-spring system.

[Matt J]

mampl are the Amplitudes of for different Frequencies (frqSet).

Different frequencies ?

[Paul-Adam]

Yes. frqSet are an array of excitation frequencies.

Answer 1

[Bjorn Gustavsson]

There are two issues you will have to solve. The first is to convert your second order ODE to two coupled first-order ODEs. The second is to wrap everything into an "error-function" for the fit to your data such that you can optimise the parameters.

The first step will be to re-write the equation-of-motion function to something like this:

function dudtdu2dt2= f(u,t)
  % First derivative/Velocity
  dudt = u(2);
  % Second derivative/acceleration   note the changes of indices    v          v 
  du2dt2 = x(4).*sin(xd.*t)./(x(2))-2.*x(1).*sqrt(x(2)./x(3)).*m.*u(2)./x(2)-u(1).*x(3)./x(2);
  dudtdu2dt2 = [dudt;du2dt2]; % Return the column vector of [v;a]
end

For the second step have a look at the solution to a similar enough problem here: monod-kinetics-and-curve-fitting. The ODEs that are solved there are completely different but the technique/programming pattern is identical.

HTH

[Paul-Adam]

Thank you for your quick answer, but after implementing the two coupled first order ODE's. I got an error because the incoming function for the odesolver must produce a column vector, which I can't create, because du2dt2 produces a series vector with 4 values instead of one. And I think there is a general problem because I am trying to fit in the frequency domain instead of the time domain.

[Torsten]

And I think there is a general problem because I am trying to fit in the frequency domain instead of the time domain.

Yes, that's the main problem. And the arguments to ODE45 are reversed:

Use

[ta,ua] = ode45(@f,ts,u0);
function dudt= f(t,u)

instead of

[ta,ua] = ode45(@(u,t) f(u,t),ts,u0);
function dudt= f(u,t)

[Bjorn Gustavsson]

@Paul-Adam, that your input-data as frequency and amplitude makes your problem completely different than the setup you've implemented. There are a good couple of other problems as well.

In your fitting-function you try to integrate the equation of motion. But you feed ode45 with a time-array with 2 components (that look as if you want them to be frequencies?). In order to get the amplitude you'll need to integrate the ode for at least a couple of period-times and then estimate the amplitude of the oscillations. In order to be able to compare those amplitude-estimates with the data you have, you will have to integrate the ode for each driving frequency and then collect all the amplitudes into one array.

In your ode you have a number or free parameters, D, k, m and F0, one problem here is that m only occurs as a denominator of k and F0 - that means that you get an identical amplitude-frequency-curve for m and 2*m if you simply scale the spring-constant and driving force similarly - this makes good physical sense, however it typically makes optimisation-functions perform poorly - because now the function will not have a localized minima but a curve in your 4-D parameterspace with that minima. This mucks up the convergence of the fitting function - it is not designed to pick a point along such constant-value curve. So the advice here is that you settle with estimating k_over_m and F0_over_m instead.

These steps can be implemented into the fitting-function, it will become a bit bigger, but you'll have to make bigger fitting-functions sooner or later anyways.

Another way to go about this specific problem is of course to calculate the analytical solution to the equation of motion and use that to calculate the amplitude-response - this is much more of manual labour and might not be what you intended to learn/practise.

[Paul-Adam]

I have allready swapped the two input arguments. ;) But there are still the problem with the column output vector.

[Bjorn Gustavsson]

Your first step is to understand what you get from the ode45-call and how that relates to the data you need to compare that output with.

To do that you should integrate the equations of motion for a long enough time to see the oscillations - then it will become clear what the two columns corresponds to - for example if you have a system of odes for dx/dt and dv_x/dt and integrate those two together then it stands to reason that the result might be x(t) and v_x(t).

You will still have to take the steps I outlined above.

[Paul-Adam]

I have already understood the actual problem. But I could not make a connection from my time based ODE to my frequency response data while using lsqcurvefit. In principle, I would have to output time ranges for several frequencies with the ODE solver and tap the amplitude for the frequency range in the steady state range. Then compare this with my frequency response. Maybe it would make more sense to work with a transfer function and use tfest?

[Bjorn Gustavsson]

Exactly! With your ODE-integration you calculate the time-variation of u for one driving frequency. From that oscillating curve of u(t) you then estimate the amplitude at that frequency. Then you repeat this procedure for each frequency of your data. That way your modified fitfun builds up the numerical curve for mampl for one set of parameters (D, k, m, and F0) once your fitfun does all of that then you can let lscurvefit do its magic.

If you have access to the symbolic toolbox you could also simply integrate the equation-of-motion, calculate the amplitude-as-a-function-of-frequency and turn that solution into a function with matlabFunction. From there it should be straightforward to fit 3 of your parameter.

You could also move this into the frequency-domain with a harmonic-expansion ansatz - that way you'd get away from onset-transients and get a closed-form expression for the frequency-variation of the amplitude.

[Paul-Adam]

I don't know how to generate an symbolic equation with relation between amplitude and frequency, but do you mean by harmonic expansion the hamonic balance method? Like

[Bjorn Gustavsson]

With harmonic expansion I mean that you can make the ansatz:

This happily ignores the onset but will give you the right long-term behaviour. The velocity and acceleration you get:

These you can then plug into your equation of motion - and notice that the complex exponentials cancels out - what remains is an algebraic equation for - the amplitude of the oscillation. You will find that it is now a complex variable, but that only means that the real and imaginary parts explain how much the oscillation is phase-shifted relative to the driving force. This equation you can solve for and calculate its absolute value for. This should be the analytical expression for the relation between the amplitude and frequency. We can do this because your equation of motion is linear, that gives us only complex exponentials that can be factored out. If you had a non-linear equation-of-motion this would not always be possible, for example if we had a term with a factor then this harmonic expansion would lead to one term with - this would couple the equations for different frequencies and be much more challenging.