Re: Re: Fitting parameters of nonlinear diff equation system
- To: mathgroup at smc.vnet.net
- Subject: [mg81567] Re: [mg81303] Re: Fitting parameters of nonlinear diff equation system
- From: "Szabolcs Horvát" <szhorvat at gmail.com>
- Date: Fri, 28 Sep 2007 02:05:42 -0400 (EDT)
- References: <fcnl8h$s85$1@smc.vnet.net> <200709181002.GAA04068@smc.vnet.net>
On 26/09/2007, Paul Thomas <paul.pgt at gmail.com> wrote:
> Hi Szabolcs,
>
> A couple things, thanks for your reply!
>
> Replacing f with f[t] in fun[] does seem to return a function. I tested this
> by just running fun with a given a,b,c and with the f[t] syntax it produces
> an interpolating function (as NDSolve should)--it doesn't with just f as I
> think it doesn't know what "f" is in that case--the output of just using "f"
> is literally just "f" with no function associated.
A few examples that demonstrate how NDSolve returns the result:
In[1]:= NDSolve[{f'[t]==f[t],f[0]==1},f,{t,0,1}]
Out[1]= {{f->InterpolatingFunction[{{0.,1.}},<>]}}
In[2]:= NDSolve[{f'[t]==f[t],f[0]==1},f[t],{t,0,1}]
Out[2]= {{f[t]->InterpolatingFunction[{{0.,1.}},<>][t]}}
In[3]:= f/.NDSolve[{f'[t]==f[t],f[0]==1},f,{t,0,1}]
Out[3]= {InterpolatingFunction[{{0.,1.}},<>]}
In[4]:= f[t]/.NDSolve[{f'[t]==f[t],f[0]==1},f[t],{t,0,1}]
Out[4]= {InterpolatingFunction[{{0.,1.}},<>][t]}
Using f[t] inside NDSolve gives the value of the
InpterpolatingFunction at the point 't'.
Using f (without the argument t) returns the function itself, not the
value of the function at a particular point.
> On the /@ priority issue, I actually changed one other thing, which maybe is
> a problem too.
> Instead of using "/@" there I used "/."--that I think does require the
> parentheses. I kept getting an error using "/@" and "/." seemed to fix it.
Here are two examples demonstrating /@ (Map) and /. (ReplaceAll).
Map applies f to each element in the list.
In[1]:= f/@{1,2,3}
Out[1]= {f[1],f[2],f[3]}
ReplaceAll can be used to replace parts of an expression with some
other expression. In this case we change 't' in f[t] to {1,2,3}.
In[2]:= f[t]/.{t->{1,2,3}}
Out[2]= f[{1,2,3}]
Both approaches may be used in the application we are discussing,
because InterpolatingFunctions behave as if they had the attribute
Listable, i.e. they accept a List as their argument, and "apply
themselves" to each element of the list. (Actually I did not know
about this before, so thanks for the hint!)
In[3]:= SetAttributes[f,Listable]
In[4]:= f[{1,2,3}]
Out[4]= {f[1],f[2],f[3]}
> The error would say that essentially a non-numerical equation was resulting
> from the /@ and it would write it out--it did seem to be doing the wrong
> equation, but it looks like the right one with /.
> Now it seems to be doing exactly what you designed it to do.
Now I learned that it is a really bad idea to provide a not fully
tested, partial example. It just leads to misunderstandings. I am
sure that your approach is correct, just slightly different than mine
(e.g. you might have used variables such as tvals and fvals slightly
differently).
Here is a full example with the very simple differential equation
f'(t) = a*f(t):
We use the values of an exponential function for fitting:
In[1]:=
tvals = Range[0, 1, .1];
fvals = 2 Exp[3 tvals];
In[3]:= fun[a_,b_]:=First[f/.NDSolve[{f'[t]==a*f[t],f[0]==b},f,{t,0,1}]]
In[4]:= opt[a_?NumericQ,b_?NumericQ]:=Norm[fun[a,b]/@tvals-fvals]
(* The more compact fun[a,b][tvals] - fvals could be used instead! *)
In[5]:= FindMinimum[opt[a,b],{{a,1},{b,1}}]
During evaluation of In[5]:= FindMinimum::lstol: The line search
decreased the step size to within tolerance specified by AccuracyGoal
and PrecisionGoal but was unable to find a sufficient decrease in the
function. You may need more than MachinePrecision digits of working
precision to meet these tolerances. >>
Out[5]= {0.0000218821,{a->3.,b->2.}}
> Once the parameters are fit I feed them back into "fun" and check the work
> essentially of the process and plot the optimization over time. For example,
> once I've done a particular fit:
>
> g[t_]:=fun[1,2,3] {pretending 1,2,3 were our best a,b,c fits}
This is where the misunderstanding came from: In my example
fun[1,2,3] returned a function (e.g. f), in your example it returns
the value of the function at point 't' (e.g. f[t]).
> g[t] is now my best fit function and matches the findminimum minimum at the
> particular t points I have. The curves also look generally as they should.
>
> Do you think I have an error in this?
No, I think that everything is OK. Here is the previous simplified
example converted to what I think you did:
In[1]:=
tvals=Range[0,1,.1];
fvals=2Exp[3tvals];
In[3]:= fun[a_,b_]:=First[f[t]/.NDSolve[{f'[t]==a*f[t],f[0]==b},f[t],{t,0,1}]]
In[4]:= opt[a_?NumericQ,b_?NumericQ]:=Norm[(fun[a,b]/.{t->tvals})-fvals]
In[5]:= FindMinimum[opt[a,b],{{a,1},{b,1}}]
During evaluation of In[5]:= FindMinimum::lstol: The line search
decreased the step size to within tolerance specified by AccuracyGoal
and PrecisionGoal but was unable to find a sufficient decrease in the
function. You may need more than MachinePrecision digits of working
precision to meet these tolerances. >>
Out[5]= {0.0000218821,{a->3.,b->2.}}
> Thanks for all your help with it, like
> I said, it seems like a very powerful approach at least for the types of
> things that I've been trying to do--very flexible!
I'm glad you found it useful!
> Paul
Szabolcs
- References:
- Re: Fitting parameters of nonlinear diff equation system
- From: Szabolcs Horvát <szhorvat@gmail.com>
- Re: Fitting parameters of nonlinear diff equation system