## Parameter Estimation for a set of Differential Equations

### Parameter Estimation for a set of Differential Equations

Hi!

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

kamal

### Parameter Estimation for a set of Differential Equations

> Hi!

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

> 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

> > Hi!

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

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

> 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

Hi!

Does anybody has experience in doing parameter estimation of
differential equations?
The system is described by only one equation of y'(t)=A*y(t). In my case
data of y(t) is known and I has to fit 3 or three parameters.

gernot