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