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]} instead of solving for {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