Re: Re: Vector-valued DE's
- To: mathgroup at smc.vnet.net
- Subject: [mg16609] Re: [mg16538] Re: Vector-valued DE's
- From: David Withoff <withoff>
- Date: Wed, 17 Mar 1999 23:55:11 -0500
- Sender: owner-wri-mathgroup at wolfram.com
> Kevin, thank you for your help.
> Unfortunately it does not help me too much.
> The problem is that my righthandside function (rhs[z[t]]) is not
> simply a function of t, but a function of all elements in z.
> Moreover, it is a computationally expensive function, which
> I want to evaluate only once for each iteration, and which
> I do not want to do symbolically (even if i could).
>
> I would like to do:
>
> NDSolve[ { z'[t]==rhs[z[t]] , z[0]=incon }, z , {t,0,1}]
>
> where z,rhs and incon are vector valued.
> They may contain in the order of 100 elements.
>
> This would mean however, that rhs would be re-evaluated
> for every element of z, after which the right element (rhs[vars][[i]])
> is chosen. That means that rhs is evaluated n times as often as
> necessary, which is unacceptable for me.
>
One way to set up NDSolve to solve a vector differential equation of this
form and to avoid re-evaluation of rhs is to use caching. The example
below uses the simple caching method describe under "caching" (or "Functions
That Remember Values They Have Found") in the Mathematica book.
In[1]:= zv = Table[z[i], {i, 5}];
In[2]:= zvt = Through[zv[t]]
Out[2]= {z[1][t], z[2][t], z[3][t], z[4][t], z[5][t]}
In[3]:= incon = Range[5] ;
In[4]:= rhs[z_] := Table[rhsv[i][z], {i, 5}]
In[5]:= rhsv[i_][p : {__?NumberQ}] := ExpensiveRHS[p][[i]]
In[6]:= ExpensiveRHS[p_] :=
(ExpensiveRHS[p] = p) /;
(Print["evaluating ExpensiveRHS at ", InputForm[p]]; True)
In[7]:= sol = NDSolve[{Thread[D[zvt, t] == rhs[zvt]],
Thread[Through[zv[0]] == incon]}, zvt, {t, 0, 1}]
evaluating ExpensiveRHS at {1., 2., 3., 4., 5.}
evaluating ExpensiveRHS at {1.000768221279597, 2.001536442559194,
> 3.002304663838792, 4.003072885118389, 5.003841106397986}
evaluating ExpensiveRHS at {1.000768811443531, 2.001537622887063,
> 3.002306434330595, 4.003075245774127, 5.003844057217659}
evaluating ExpensiveRHS at {1.001537623793816, 2.003075247587633,
[[[[[ invervening Print output omitted ]]]]]
evaluating ExpensiveRHS at {2.718285577657515, 5.436571155315031,
> 8.154856732972547, 10.87314231063006, 13.59142788828757}
Out[7]= {{z[1][t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
> z[2][t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
> z[3][t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
> z[4][t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
> z[5][t] -> InterpolatingFunction[{{0., 1.}}, <>][t]}}
In[8]:= Plot[Evaluate[zvt /. sol], {t, 0, 1}]
Out[8]= -Graphics-
The Print output is used here only to demonstrate that ExpensiveRHS
is evaluated only once for each new value of the vector argument.
You will likely want to drop this when it is no longer needed for
diagnostic purposes.
> It is a pity that the internals of (N)DSolve are "internal"
> indeed, because it probably would not be hard to solve on a
> slightly lower level.
> NDSolve should be able to work with list-equations easily.
Access to the internals would not be of much help here, but yes, it
would nevertheless be nice to make this easier. This suggestion has
been considered and facilities to make this easier will likely be added
for some future version of Mathematica.
> Dirk
>
> "Kevin J. McCann" wrote:
> >
> > I think the following is what you are after:
> >
> > Thread[{a,b,c}=={d,e,f}]
> >
> > {a==d,b==e,c==f}
> >
> > To use in DE's:
> >
> > Your vector function:
> >
> > z[t_]={a[t],b[t],c[t]}
> >
> > and
> >
> > rhs[t_]={d[t],e[t],f[t]}
> >
> > Now
> >
> > z''[t]==rhs[t]
> >
> > doesn't quite work, but
> >
> > Thread[z''[t]==rhs[t]]
> >
> > {a''[t]==d[t],b''[t]==e[t],c''[t]==f[t]}
> >
> > Similarly for the IC's, then just use Flatten to get a single list with DE +
> > IC.
> >
> > Hope this helps,
> >
> > Kevin
> >
Dave Withoff
Wolfram Research