Re: Re: NonlinearRegress and numerical functions...
- To: mathgroup at smc.vnet.net
- Subject: [mg20212] Re: [mg20197] Re: [mg20161] NonlinearRegress and numerical functions...
- From: "Carl K.Woll" <carlw at fermi.phys.washington.edu>
- Date: Wed, 6 Oct 1999 21:06:30 -0400
- Organization: Department of Physics
- References: <199910040107.VAA17261@smc.vnet.net.> <199910050804.EAA22807@smc.vnet.net.>
- Sender: owner-wri-mathgroup at wolfram.com
Daniel,
There is a way to get NonlinearFit to work even when using NDSolve in your
model. One way, which is essentially identical to what you did, is to use
Method->FindMinimum and give NonlinearFit two starting values. However,
suppose you wanted to include gradient information. One way to do this is to
use the Gradient option of NDSolve, but I prefer a different method. Teach
Mathematica how to take derivatives of your model! In your example, we had
Clear[f]
f[a_?NumberQ, b_?NumberQ, t_?NumberQ] :=
(y[t1] /. First[NDSolve[{y'[t1]==-a*y[t1], y[0]==b}, y[t1],
{t1,0,100}, AccuracyGoal->150, MaxSteps->2000]]) /. t1->t
Now, we need to figure out what the derivative of f[a,b,t] with respect to a
is.
fa[a_?NumberQ, b_?NumberQ, t_?NumberQ] :=
ya[t] /. First[
NDSolve[{ya'[t1] == -a ya[t1] - y[t1], y'[t1] == -a y[t1], ya[0] == 0,
y[0] == b}, {y, ya}, {t1, 0, 100}, AccuracyGoal -> 150,
MaxSteps -> 2000]]
Here I just differentiated the ODE with respect to a to obtain a new system of
ODEs. Next, we need to teach Mathematica what this derivative is
Derivative[1, 0, 0][f][a_, b_, t_] := fa[a, b, t]
Now we can use NonlinearFit
In[25]:=
NonlinearFit[data, f[a, 6, t], t, {a, 1}]
Out[25]=
f[3., 6, t]
An extension of this idea should work for Lars' problem.
Carl Woll
Physics Dept
U of Washington
Daniel Lichtblau wrote:
> Lars Ragnarsson wrote:
> >
> > Hi
> >
> > I've got a big problem making Mathematica realize that my function is
> > numerical. The following problem illustrates it in a rather simplified
> > way:
> >
> > g[t1_] := y[t] /. DSolve[{y'[t] == -y[t] 3, y[0] == 6}, y[t], t][[1]] /.
> >
> > t -> t1
> >
> > data = {0, 10, 20, 100};
> >
> > fit = Transpose[{data, g[data]}];
> >
> > f[t1_, a_,b_] := (y[t] /.NDSolve[{y'[t] == -y[t] a, y[0] == b}, y[t],
> > {t, 0, 100}][[1]]) /.t -> t1
> >
> > NonlinearRegress[data, f[t, a, 6], {t}, {a, {2.9, 3.2}, 2, 4}]
> >
> > It should be rather straightforward to solve this, but NonlinearRegress
> >
> > seem to evaluate the function 'f' first and the results:
> >
> > NDSolve::"ndnum": "Encountered non-numerical value for a derivative at \
> >
> > \!\(t\) == \!\(2.680466916142905`*^-274\)."
> >
> > etc....
> >
> > I've tried to use Unevaluated and Hold but nothing works!
> >
> > Any help would be greatly appreciated!
> >
> > Regards
> >
> > Larske Ragnarsson
> >
> > _________________________________________________
> > Larske Ragnarsson
> > Phone: +46 (0)31 7721867
> > Fax: +46 (0)31 7723622
> > Chalmers University of Technology
> > Solid State Electronics Laboratory
> > Department of Microelectronics ED
> > S-412 96 Guteborg, SWEDEN
> > _________________________________________________
>
> I'll change the notation just a bit. To start we define our function,
> create a bunch of data points, and show, with NonlinearFit, that we can
> recover the function.
>
> g[a_, b_, t_] := (y[t1] /.
> First[DSolve[{y'[t1]==-a*y[t1], y[0]==b}, y[t1], t1]]) /. t1->t
> pnts = {0, 10, 20, 100};
> data = {#, g[3,6,t] /. t->#}& /@ pnts
> <<Statistics`
>
> In[5]:= NonlinearFit[data, g[a,6,t], t, {a,1}]
> 6
> Out[5]= -----
> 3. t
> E
>
> Now suppose instead that we have a function defined by (numerical)
> solutions to a differential equation. I assume the example given was for
> purposes of comparison to a known solution though of course in general
> we do not have this recourse. Let me point out some issues, at least one
> of which is substantial.
>
> First, your numeric differential equation solution is not very accurate;
> to obtain a result similar to the exact solution for x as large as 100
> you need to fiddle with the NDSolve options (maybe this numeric error
> was intentional?).
>
> Clear[f]
> f[a_?NumberQ, b_?NumberQ, t_?NumberQ] :=
> (y[t1] /. First[NDSolve[{y'[t1]==-a*y[t1], y[0]==b}, y[t1],
> {t1,0,100}, AccuracyGoal->150, MaxSteps->2000]]) /. t1->t
>
> Let us check the relative errors at our data points.
>
> In[9]:= Map[(f[3,6,#]-g[3,6,#])/g[3,6,#]&, pnts]
> Out[9]= {0., -0.0000162022, -0.0000385361, -0.000219838}
>
> Not bad.
>
> Now for a more serious issue. You want to do a nonlinear function fit to
> a model that is, in effect, defined by a program. While this is a
> reasonable thing to want, I am not certain it can be done directly in
> Mathematica (I confess to no great expertise with NonlinearFit/Regress,
> so I may be mistaken). Anyway, you can instead formulate this as a
> minimization problem and use FindMinimum to solve for the parameter
> value. We do this explicitly below. Again, some option fiddling is
> helpful in order to get a good result.
>
> In[15]:= objfunc = Apply[Plus, Map[(f[a, 6, #[[1]]] - #[[2]])^2 &,
> data]];
>
> In[16]:= FindMinimum[objfunc, {a,2,2.1}, AccuracyGoal->50,
> MaxIterations->200, WorkingPrecision->100]
> -53
> Out[16]= {3.79385 10 , {a -> 2.99999837977173152356769137355}}
>
> Given the disparities between f and g at the data points, this is likely
> about as accurate a solution as can be obtained; to do better one would
> have to get more accurate results from the model function f and perhaps
> also use/request higher accuracy/precision in FindMinimum.
>
> Daniel Lichtblau
> Wolfram Research
- References:
- NonlinearRegress and numerical functions...
- From: Lars Ragnarsson <loke@ic.chalmers.se>
- Re: NonlinearRegress and numerical functions...
- From: Daniel Lichtblau <danl@wolfram.com>
- NonlinearRegress and numerical functions...