Re: subtract a list of interpolating functions from another
- To: mathgroup at smc.vnet.net
- Subject: [mg48550] Re: subtract a list of interpolating functions from another
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Sat, 5 Jun 2004 07:18:51 -0400 (EDT)
- Organization: The University of Western Australia
- References: <c9pfod$t27$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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