Parameter Estimation for a set of Differential Equations

Parameter Estimation for a set of Differential Equations

Post by Kamal Ajitsari » Wed, 23 Apr 1997 04:00:00



Hi!

i would like to know how can you estimate the parameters for a set of
differential equations using levenberg method.

thanks in advance

kamal

 
 
 

Parameter Estimation for a set of Differential Equations

Post by Bradley Jon » Wed, 30 Apr 1997 04:00:00




> Hi!

> i would like to know how can you estimate the parameters for a set of
> differential equations using levenberg method.

> thanks in advance

> kamal

One of my examples in the Statistics Talk at the MATLAB Conference in October
just happens to be an example of what you want. :-)

At the risk of stealing my own thunder, here is a simple example of fitting
data from a damped harmonic oscillator using the differential equation as
the model.

This example uses 3 tools: inline functions, ode45, and nlinfit from the
Statistics Toolbox.

First create a global inline function:

global decayode
decayode = inline('[y(2);-(k*y(2)+w*y(1))]','t','y','flag','k','w')

Now you will need a prediction function for nlinfit to use for the iterative
Levenberg-Marquardt procedure.

function ypred = decayfit(beta,t)
global decayode

k = beta(1);
w = beta(2);
y0 = beta(3);
v0 = beta(4);

[tnew,ypred] = ode45(decayode,t,[y0 v0],[],k,w);
ypred = ypred(:,1);

Note: decayfit calls ode45 with our inline function as an argument.

Here is a script to generate some pseudo-data, fit them, and plot the results.

% Supply number of points.
n = 25;
% Supply parameter values.
a = .5;
b = 4;
y0 = 0.00001;
v0 = 1;

% Generate data.
t = sort(unifrnd(0,2*pi,n,1));
[t,ytrue] = ode45(decayode,t,[y0 v0],[],a,b);
y = ytrue(:,1) + normrnd(0,0.1);

% Fit data to differential equation.
[beta,r,J] = nlinfit(t,y,'decayfit',[a;b;y0;v0]);  % <-- Here is nlinfit.
figure;
plot(t,y,'+')
tnew = linspace(0,2*pi,1000);

% Plot original data with fitted function.
[tnew,ypred] = ode45(decayode,[t(1) 4*pi],[beta(3) beta(4)],[],beta(1),beta(2));
hold on,plot(tnew,ypred(:,1),'r-')

I hope you find this useful.

--


|    24 Prime Park Way                http://www.mathworks.com    |
|    Natick, MA 01760-1500                   ftp.mathworks.com    |
==== Tel: 508-653-1415 ==== Fax: 508-653-6971 =====================

 
 
 

Parameter Estimation for a set of Differential Equations

Post by Torsten Schuetz » Wed, 30 Apr 1997 04:00:00





> > Hi!

> > i would like to know how can you estimate the parameters for a set of
> > differential equations using levenberg method.

> > thanks in advance

> > kamal

> One of my examples in the Statistics Talk at the MATLAB Conference in October
> just happens to be an example of what you want. :-)

> At the risk of stealing my own thunder, here is a simple example of fitting
> data from a damped harmonic oscillator using the differential equation as
> the model.

> This example uses 3 tools: inline functions, ode45, and nlinfit from the
> Statistics Toolbox.

Since we have no derivative information for nlinfit the Jacobian has to
be approximated by finite differences, i.e. we have to solve the
differential equation for different, nearby parameter values. How do you
assure that the automatic step size control of ode45 doesn't disturb
your derivative information. More specific:
the forward difference [ f(x0 + \epsilon ) - f(x0)] / \epsilon gives
useful derivative information only if the function f (the dynamic
system) is computed with the same step size sequence. How do you assure
this in your above example? Are there any options to say ode45 ,,use the
same step sizes as in the previuos run''?
I think it is this point that makes parameter estimation in differential
equations nontrivial. Either you must turn off the automatic step size
control or you must solve all equations (for all disturbed parameter
values) at once.

> I hope you find this useful.

> --


> |    24 Prime Park Way                http://www.mathworks.com    |
> |    Natick, MA 01760-1500                   ftp.mathworks.com    |
> ==== Tel: 508-653-1415 ==== Fax: 508-653-6971 =====================

Would be interesting to know if there any such options in your ode
solvers.

Torsten
--
Torsten Schuetze                        Dresden University of Technology
Phone: (+49 351) 463 4084               Department of Mathematics
Fax:   (+49 351) 463 4268               Institute of Numerical
Mathematics

WWW: http://www.math.tu-dresden.de/~schuetze

 
 
 

Parameter Estimation for a set of Differential Equations

Post by Murali Yeddanapud » Wed, 30 Apr 1997 04:00:00



> Since we have no derivative information for nlinfit the Jacobian has to
> be approximated by finite differences, i.e. we have to solve the
> differential equation for different, nearby parameter values. How do you
> assure that the automatic step size control of ode45 doesn't disturb
> your derivative information. More specific:
> the forward difference [ f(x0 + \epsilon ) - f(x0)] / \epsilon gives
> useful derivative information only if the function f (the dynamic
> system) is computed with the same step size sequence. How do you assure
> this in your above example? Are there any options to say ode45 ,,use the
> same step sizes as in the previuos run''?
> I think it is this point that makes parameter estimation in differential
> equations nontrivial. Either you must turn off the automatic step size
> control or you must solve all equations (for all disturbed parameter
> values) at once.

Hi Torsten:

Indeed. If the cost function f(x) being minimized is computed using
the result from an ODE solver, then the finite difference approximation
to the Jacobian, computed from successive calls to the cost function
for f(x) and f(x+dx), can become quite inaccurate.

We use the following approach (possibly similar to what you suggest in
your last statement) in the Nonlinear Control Design Toolbox, that does
not restrict the steps taken by the ODE Solver when computing f(x+dx),
but avoids the above problem.

Instead of letting the optimization routine compute the Jacobian, we
compute it in a seperate function (gradfun) by solving for f(x) and
f(x+dx) *simultaneously*.  In gradfun, we pass an augmented system to
the ODE solver, consisting of the original system looking at the
parameter values in x and a perturbed version looking at x+dx.

Murali

--
Murali Yeddanapudi
The MathWorks, Inc