Re: subtract a list of interpolating functions from another

• To: mathgroup at smc.vnet.net
• Subject: [mg48604] Re: subtract a list of interpolating functions from another
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Mon, 7 Jun 2004 05:33:36 -0400 (EDT)
• Organization: The University of Western Australia
• References: <c9pfod\$t27\$1@smc.vnet.net> <c9sdmq\$c78\$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```In article <c9sdmq\$c78\$1 at smc.vnet.net>,
sean_incali at yahoo.com (sean kim) wrote:

> I can't seem to get your suggestion to work. It brings back lotta
> plotting error. so I looked at,
>
> Evaluate[(Subtract@@@stepde)/.ndsolstep]
> Evaluate[(Subtract@@@nostepdes)/.ndsol]
>
> first one brings back,
>
> {0.05 InterpolatingFunction[{{0.,400.}},<>][t]
> InterpolatingFunction[{{0.,400.}},<>][t]+ b'[t],
> 0.1 a0[t] InterpolatingFunction[{{0., 400.}}, <>][t] -
> 0.05\InterpolatingFunction[{{0., 400.}}, <>][t]
> InterpolatingFunction[{{0., 400.}}, <>][t] + x'[t],
> -0.1 a0[t] InterpolatingFunction[{{0., 400.}}, <>][t] + 0.05
> \InterpolatingFunction[{{0., 400.}}, <>][t]
> InterpolatingFunction[{{0., 400.}}, <>][t] + y'[t]}
>
> It doesn't look like they are evaluated...

As you can see, b'[t] has not been evaluated. This happens because, for
some reason, you've now solved for

{b[t], x[t], y[t]}

{b, x, y}

The difference _is_ important. When you solve for {b,x,y}, derivatives
such as b'[t] are easily evaluated. For example, with

sol1 = NDSolve[{h'[x] == h[x], h[0] == 1}, h[x], {x,0,5}]
and
sol2 = NDSolve[{h'[x] == h[x], h[0] == 1}, h, {x,0,5}]

compare
h'[1] /. sol1
to
h'[1] /. sol2

In sol2 solving for h means that Mathematica knows how to evaluate
derivatives of h.

> What is Subtract@@@stepde supposed to do up there?

stepde is a system of equations. You want to turn these into a set of
functions that should evaluate to zero. For example,

Subtract @@@ {a == b, c == d}

> below is the code that was ran.
>
>
> k1 = 0.1;
> k2 = 0.05;
> a0[t_]:= 0.5/; t< 0 ;
> a0[t_]:= 0.5/;0 <= t<=20 ;
> a0[t_]:= 0.5/; 20<=t<= 60 ;
> a0[t_]:= 0.5/; 60<=t<= 400;
>
> a=0.5;
>
> stepde = {b'[t]== -k2 b[t] y[t], x'[t]== -k1 a0[t] x[t] + k2 b[t]
> y[t], y'[t]== k1 a0[t] x[t] - k2 b[t] y[t]};
> nostepde = {b'[t]== -k2 b[t] y[t], x'[t]== -k1 a x[t] + k2 b[t] y[t],
> y'[t]== k1 a  x[t] - k2 b[t] y[t]};
>
> ndsolstep =  NDSolve[ Join[stepde, { b[0] == 1, x[0] == 1, y[0] ==
> 0}], {b[t], x[t], y[t]}, {t, 0,0,20, 60, 400,400},
> Method->ExplicitRungeKutta][[1]] ;
> ndsol = NDSolve[Join[ nostepde,  {b[0] == 1, x[0] == 1, y[0] == 0}],
> {b[t], x[t], y[t]}, {t, 0,400}][[1]] ;

Solve for {b,x,y} here.

Cheers,
Paul

> SetOptions[Plot,PlotRange -> All];
>
> Evaluate[(Subtract@@@stepde)/.ndsolstep]
> Evaluate[(Subtract@@@nostepdes)/.ndsol]
>
> Plot[Evaluate[(Subtract@@@stepde)/.ndsolstep],{t,0,400}]
> Plot[Evaluate[(Subtract@@@nostepdes)/.ndsol],{t,0,400}]
>
> --- Paul Abbott <paul at physics.uwa.edu.au> wrote:
> > In article <c9pfod\$t27\$1 at smc.vnet.net>,
> >  sean_incali at yahoo.com (sean kim) wrote:
> >
> > > still stepwise ode accuracy related question,
> > but...
> > >
> > > consider two lists of three interpolating
> > functions.
> > >
> > > like below.
> > >
> > > k1 = 1/10; k2 = 1/20;
> > >
> > > a0[t_] :=  0.5 /; t < 0 ;
> > > a0[t_] := 0.5 /; 0 <= t <=20 ;
> > > a0[t_] :=  0.5 /; 20 <= t <=60 ;
> > > a0[t_] := 0.5 /; 60 <= t <=400;
> > >
> > > a = 0.5;
> > >
> > > ndsolstep = NDSolve[{ b'[t] == -k2 b[t] y[t],
> > x'[t] == -k1 a0[t]  x[t]
> > > + k2 b[t] y[t], y'[t] == k1 a0[t] x[t] - k2 b[t]
> > y[t],  b[0] == 1,
> > > x[0] == 1, y[0] == 0}, {b, x, y}, {t, 0, 0, 20,
> > 60, 400, 400}, Method
> > > -> ExplicitRungeKutta][[1]]
> > >
> > > ndsol = NDSolve[{ b'[t] == -k2 b[t] y[t], x'[t] ==
> > -k1 a  x[t] + k2
> > > b[t] y[t], y'[t] == k1 a x[t] - k2 b[t] y[t],
> > b[0] == 1, x[0] == 1,
> > > y[0] == 0}, {b, x, y}, {t, 0, 400}][[1]]
> > >
> > >
> > > will give two lists of
> > >
> > > {b -> InterpolatingFunction[{{0., 400.}}, <>],
> > > x -> InterpolatingFunction[{{0., 400.}}, <>],
> > > y -> InterpolatingFunction[{{0., 400.}}, <>]}
> > >
> > > one from normal system and another from stepwise
> > defined( which has
> > > Rob Knapp's fix in it) they should be same if not
> > very close.
> >
> > And they are.
> >
> > > I thought maybe I would take a value of
> > interpolating function at time
> > > poiints and subtract to see the differences. (to
> > check how close they
> > > are)
> >
> > Instead, why not just substitute the solutions back
> > into the
> > differential equations (Mathematica knows how to
> > compute derivatives of
> > InterpolatingFunctions) and see how well they are
> > satisfied:
> >
> >   des = {b'[t] == -k2 b[t] y[t], x'[t] == -k1 a0[t]
> > x[t] + k2 b[t] y[t],
> >    y'[t] == k1 a0[t] x[t] - k2 b[t] y[t]};
> >
> >   SetOptions[Plot, PlotRange -> All];
> >
> >   Plot[Evaluate[(Subtract @@@ des) /. ndsolstep],
> > {t, 0, 400}];
> >
> >   Plot[Evaluate[(Subtract @@@ des) /. ndsol], {t, 0,
> > 400}];
> >
> > Cheers,
> > Paul
> >
> > --
> > Paul Abbott                                   Phone:
> > +61 8 9380 2734
> > School of Physics, M013                         Fax:
> > +61 8 9380 1014
> > The University of Western Australia      (CRICOS
> > Provider No 00126G)
> > 35 Stirling Highway
> > Crawley WA 6009
> > mailto:paul at physics.uwa.edu.au
> > AUSTRALIA
> > http://physics.uwa.edu.au/~paul
> >
>

--
Paul Abbott                                   Phone: +61 8 9380 2734
School of Physics, M013                         Fax: +61 8 9380 1014
The University of Western Australia      (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009                      mailto:paul at physics.uwa.edu.au
AUSTRALIA                            http://physics.uwa.edu.au/~paul

```

• Prev by Date: Re: what actually is in the WRI "functions" database?
• Next by Date: Hurst exponent (R/S)
• Previous by thread: Re: subtract a list of interpolating functions from another
• Next by thread: A module to write a module