Re: how use NDSolve with an ODE having parameters
- To: mathgroup at smc.vnet.net
- Subject: [mg39711] Re: how use NDSolve with an ODE having parameters
- From: "Carl K. Woll" <carlw at u.washington.edu>
- Date: Mon, 3 Mar 2003 04:26:05 -0500 (EST)
- Organization: University of Washington
- References: <b3rsp0$dmt$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Murray,
This topic has come up before, and I will give you my favorite approach.
Suppose your model is the solution to the differential equation given in
your post (more complicated ODEs can be easily accomodated):
f[a_?NumericQ,b_?NumericQ,t_?NumericQ]:=
y[t]/.NDSolve[{y'[s]==a y[s]+b,y[0]==1},y,{s,0,1}][[1]]
Note that I've restricted the input parameters (a,b) here to strictly
numeric quantities to prevent NDSolve from being fed symbolic inputs as you
mention in your post. I've also restricted the variable t to a strictly
numeric quantity for a reason I'll explain later. Note that the variable
used inside the NDSolve is s instead of t. If s has a value, then you will
need to either choose a different variable of integration, or include the
NDSolve inside a Block or Module to make sure that s doesn't already have a
value.
Now, to use NonlinearFit, we need to teach Mathematica how to take
derivatives of this function. First, what are the derivatives of f[a,b,t]
with respect to a and b? Simply:
fa[a_?NumericQ,b_?NumericQ,t_]:=
ya[t]/.
NDSolve[
{y'[s]==a y[s]+b,ya'[s]==a ya[s]+y[s],y[0]==1,ya[0]==0},
{y,ya},
{s,0,1}
][[1]]
fb[a_?NumericQ,b_?NumericQ,t_]:=
yb[t]/.
NDSolve[
{y'[s]==a y[s]+b,yb'[s]==a yb[s]+1,y[0]==1,yb[0]==0},
{y,yb},
{s,0,1}
][[1]]
What I've done in the above is to simply differentiate the original ODE (as
well as the initial condition) with respect to the parameters, and performed
an NDSolve on the system of ODEs consisting of the original and the
differentiated ODEs. Note: Strictly speaking, the definition of fb doesn't
require the original ODE, but for pedagogical reasons I included them
anyway. Also, here there is no pressing need to restrict the variable t to
numeric quantities and so I haven't done so.
Now that we've figured out how to differentiate the model with respect to
the parameters, we need to teach Mathematica this information. So:
Derivative[1,0,0][f]=fa;
Derivative[0,1,0][f]=fb;
Pretty simple. For example,
In[21]:=
D[f[a,b,t],a]
Out[21]=
fa[a, b, t]
Now it's time to use NonlinearRegress, which we do in the usual way. Load in
the package:
In[22]:=
<<Statistics`NonlinearFit`
Create some data:
In[24]:=
data=Table[{t,f[1,2,t]+10^-3 Random[]},{t,0,1,1/100}];
Note that I've included some noise here. Finally, it's time to run
NonlinearRegress. However, rather than running NonlinearRegress, I'll use
NonlinearFit here, as the information one gets from NonlinearRegress would
probably look pretty ugly in this post.
In[25]:=
NonlinearFit[data,f[a,b,t],t,{a,b}]
Out[25]=
f[0.999229, 2.00225, t]
Note here that if I hadn't restricted the input to f to be strictly numeric
in t as well as a and b, then we would have gotten an InterpolatingFunction
here instead, as in
InterpolatingFunction[{{0., 1.}}, <>][t]
so that the parameter values a and b would not have been reported (not a
problem with NonlinearRegress though).
I hope the steps above are clear, and that translating the above procedure
to more complicated ODEs that Mathematica can only solve numerically is
transparent. Basically, there are three additional steps:
1. Define the model with both the parameters and varaibles restricted to
numeric quantities.
2. Define the derivatives of the model with respect to the parameters.
3. Teach Mathematica about these derivatives.
Then, go ahead and use NonlinearRegress as usual.
There are other approaches, but most of them require that you either use the
FindMinimum method of NonlinearRegress, or that you abandon NonlinearRegress
(and all of it's statistical feedback) and use FindMinimum directly. I think
my method allows one to use the full power of NonlinearRegress in a simple
and very natural way. Feel free to email me with a specific example if you
wish.
Carl Woll
Physics Dept
U of Washington
"Murray Eisenberg" <murraye at attbi.com> wrote in message
news:b3rsp0$dmt$1 at smc.vnet.net...
> This is a simplification of a question asked by a colleage. He wants to
> use as the model function argument to NonlinearRegress (from
> Statistics`NonlinearFit1) a solution of an initial-value problem for a
> differential equation, where the differential equation depends on a
> parameter.
>
> The catch is that the differential equation cannot be solved explicitly,
> so he has to resort to solving the initial-value problem by means of
> NDSolve. Of course, NDSolve will not do anything if the differential
> equation involves symbolic parameters. Thus the IDEA of what he wants
> to do is to use the "resulting function" from something like
>
> y[t]/.First@NDSolve[{y'[t] == a y[t] + b, y[0] == 1.}, y[t], {t, 0.,
1.}]
>
> -- where two parameters a and b are involved -- as the model. Of course
> if NDSolve above is changed to DSolve, no difficulty. But in the ACTUAL
> problem at issue, with a much more complicated differential equation,
> DSolve does nothing.
>
> Is there some way to make this work?
>
> There are evidently two difficulties:
>
> (1) How to deal with NDSolve when the differential equation involves
> parameters (perhaps there's something regarding use of Hold that will
> help?); and
>
> (2) For each pair of particular values of the parameters, the result
> from NDSolve is an InterpolatingPolynomial object and NOT the sort of
> "expression in the variable" that seems to be required for the model
> argument to NonlinearRegress. How should the InterpolatingPolynomial
> object be massaged to allow it to be used as an ordinary expression in
> the variable?
>
> --
> Murray Eisenberg murray at math.umass.edu
> Mathematics & Statistics Dept.
> Lederle Graduate Research Tower phone 413 549-1020 (H)
> University of Massachusetts 413 545-2859 (W)
> 710 North Pleasant Street
> Amherst, MA 01375
>
>