how use NDSolve with an ODE having parameters

[Murray Eisenberg]

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 mur...@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

[Kevin J. McCann]

Murray,

I have sent you a nb that is based on an old Mathematica Tech Note. It
should show you how to do a fit. The nb does a least squares fit, but you
could modify it to use one of the statistical packages.

Kevin

[Selwyn Hollis]

Murray,

I can't get NonlinearRegress to work either, but using FindMinimum on
the least squares error seems to work. Here's the general idea:

model[x_, a_, b_] := (y /. First@NDSolve[{ your DE, y[0] == initval},
y, {t, 0, 2}])[x]

LSE[a_, b_] := Sum[(data[[i, 2]] - model[data[[i, 1]], a, b])^2, {i,
Length[data]}]

FindMinimum[LSE[a, b], {a, a1, a2}, {b, b1, b2}]

---
Selwyn Hollis

[Reza Malek-Madani]

Dear Murray: Is this along the lines of what you're looking for?

f[omega_, a_, b_] :=
NDSolve[{x'[t] == y[t], y'[t] == -Sin[x[t]] - 0.1*y[t] + Cos[omega
t], x[0] == a, y[0] == b}, {x, y}, {t, 0, 10}]
sol = f[6, 1, 1];
newx[t_] = First[x[t] /. sol[[1]]];
Print[NIntegrate[(newx[t])^2, {t, 0, 10}]];
g[omega_] := Block[{sol}, sol = f[omega, 1, 1]; newx[t_] = x[t] /. sol;
First[newx[3]]]
Plot[g[omega], {omega, 0.1, 3}]

This solves the usual pendulum equation with omega and initial conditions
as parameters, computes the L^2 norm of x for a specific set of
parameters, and plots the x value at t=3 as omega varies.

Rez

-------------------------------------------------------------------------
Reza Malek-Madani Director of Research
Research Office, MS 10m Phone: 410-293-2504 (FAX -2507)
589 McNair Road DSN: 281-2504
U.S. Naval Academy Nimitz Room 17 in ERC
Annapolis MD 21402-5031 Email: rese...@usna.edu

--------------------------------------------------------------------------

[Murray Eisenberg]

Thanks, I'll pass this on to my colleague. Actually, he was hoping to
make use of the various information that NonlinearRegress provides with
the option RegressionReport.

Selwyn Hollis wrote:
> Murray,
>
> I can't get NonlinearRegress to work either, but using FindMinimum on
> the least squares error seems to work. Here's the general idea:
>
> model[x_, a_, b_] := (y /. First@NDSolve[{ your DE, y[0] == initval},
> y, {t, 0, 2}])[x]
>
> LSE[a_, b_] := Sum[(data[[i, 2]] - model[data[[i, 1]], a, b])^2, {i,
> Length[data]}]
>
> FindMinimum[LSE[a, b], {a, a1, a2}, {b, b1, b2}]
>
> ---
> Selwyn Hollis
>
>
>

[Carl K. Woll]

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" <mur...@attbi.com> wrote in message
news:b3rsp0$dmt$1...@smc.vnet.net...