Message Boards Message Boards

0
|
12512 Views
|
4 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Fitting data with NMinimize to a ODE system model solved with NDSolve

I'm trying to define a general procedure to fit experimental data like {{t1,x11,x21...}, {t2,x12,x22...}, {t3,x13,x23...}...} where ti is the time and xj are different observable magnitudes of the system.

I want to fit the data to a model given by a set of differential equations combining for that the NMinimize and NDSolve functions.

I have considered the following simplified problem.

1) My experimental data are just given by one observable magnitude as a function of time:
Xdat={{0.,0.},{2.5,0.274535},{5.,0.402493},{7.5,0.443142},{10.,0.434248},{12.5,0.39945},{15.,0.353193},{17.5,0.304003},{20.,0.256645},{22.5,0.213543},{25.,0.175702},{27.5,0.143293},{30.,0.116034},{32.5,0.0934173},{35.,0.0748496},{37.5,0.0597335},{40.,0.0475107},{42.5,0.0376818},{45.,0.0298143},{47.5,0.0235409},{50.,0.0185546},{52.5,0.0146022},{55.,0.0114765},{57.5,0.00900965},{60.,0.00706603},{62.5,0.00553691},{65.,0.00433544},{67.5,0.00339245},{70.,0.00265304},{72.5,0.00207375},{75.,0.00162023},{77.5,0.0012654},{80.,0.000987946},{82.5,0.000771099},{85.,0.000601692},{87.5,0.000469395},{90.,0.000366112},{92.5,0.000285506},{95.,0.000222611},{97.5,0.000173548},{100.,0.000135282}}
ListPlot[Xdata, Joined -> True, Frame -> True, PlotRange -> All]


2) The next step I define the model given by a set of ODE's
model[p_List] := {x1'[t] == -p[[1]] x1[t],
  x2'[t] == p[[1]] x1[t] - p[[2]] x2[t], x3'[t] == p[[2]] x2[t],
  x1[0] == 1, x2[0] == 0, x3[0] == 0}

where the parameters p[[1]] and p[[2]] are the parameters to be optimized with NMinimize. I have decided to define the model in this way to allow the change of the name of the parameters.

3) The next step was the definition of the objective function where in this simple case reduces to the minimization of the sum of residuals, but in other problems involves several observable/state variables:
 fun[data_, ode_] :=
   Module[{nsol, x2teo, phi, yval, xval, tmax, param, nsolall},
    {xval, yval} = Transpose[data];
    tmax = Max[xval];
    nsolall = NDSolve[ode, {x1[t], x2[t], x3[t]}, {t, 0, tmax}] ;
    nsol = x2[t] /. nsolall[[1]];
    x2teo = Evaluate[nsol] /. t -> xval;
    phi = Total[(x2teo - yval)^2]
    ];

I need a general function like fun[•] because I would like to try different models with the same data set without the requirement of writing a particular objective function for each model.

4) The function works properly with numerical arguments giving the value of the objective function at some point of the parameter space.
In[9]:= fun[Xdata, model[{0.01, 0.05}]]
Out[9]= 0.998479
ContourPlot[fun[Xdata, model[{x, y}]], {x, 0, 0.5}, {y, 0, 0.5}]



The function shows a well defined minima.

5) The last step is the use of NMinimize to obtain the best fit parameters. I have try to use the function as:
NMinimize[fun[Xdata, model[{p1, p2}]], {p1, p2}]
The function returns several errors on the evaluations basically saying that the function fun[•] is not numerical.

I'm not sure if the problem is to understand how the NMinimize function works with symbols but I cannot fix this errors.
The use of Mathematica functions like FindFit[•] or NonlinearModelFit[•] are not feasible because are valid for problems where only there is one observable magnitude. The problem I'm trying to solve involve several state/observable variables and is not as simple as the precedent example.

I'll appreciate any comment and suggestion about the procedure or code that want to develop.
Regards,
Javier
4 Replies
It might be possible to do what you want by using Repeated Patterns.

http://reference.wolfram.com/mathematica/tutorial/RepeatedPatterns.html
POSTED BY: Frank Kampas
Hello again, Frank
I have tested your proposal defining the function
fun1[p1_?NumberQ, p2_?NumberQ] := fun[Xdata, model[{p1, p2}]]
and then, when I use NMinimize gives the right answer...
In[9]:= NMinimize[fun1[p1, p2], {p1, p2}]
Out[9]= {4.95555*10^-15, {p1 -> 0.15, p2 -> 0.1}}
Thnak you for your suggestion. That works.
But, I'm trying to generalize the use of NMinimize for a List of parameters. When we define fun1[p1,p2], this function has two arguments, or in general I need 'n' arguments for a model of 'n' parameters. I'm wondering if there is a way to define a function for NMinimize where the only argument is a list of parameters.
Something like:
fun2[p_List]:=fun[Xdata,model[{p[[1]],p[[2]],p[[3]]...}]
and then minimizing like:
NMinimize[fun2[p],p]
but until I have tried, that doesn't work.
Thanks, again because I can explore other possibilities with your suggestion
Regards
Javier
I defined a function
fun1[p1_?NumberQ, p2_?NumberQ] := fun[XData, model[{p1, p2}]]
using your data.  When I try to evaluate it, I get error messages such as
Set::shape: "Lists {xval$8593,yval$8593} and Transpose[XData] are not the same shape."
That may give some clue as to the error messages you are seeing when you try to use NMinimize.
It might be a good idea to start by putting the data in the function rather than making it an argument to the function call.
POSTED BY: Frank Kampas
Hi Frank
Sorry, but I have noticed that the name given to the data set  was 'Xdat'  instead  'Xdata'  used in other places of the code. Its a mistake.
Thanks
Javier
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract