Bug: too subtle for me
- To: mathgroup at smc.vnet.net
- Subject: [mg2556] Bug: too subtle for me
- From: harrison at helios.physics.utoronto.ca (David Harrison)
- Date: Mon, 20 Nov 1995 01:14:34 -0500
- Organization: University of Toronto - Dept. of Physics
Dear Mathgroup:
I have a nasty bug that I've been unable to penetrate despite a day
and a half's worth of effort. Often someone fresh can look and see
what the problem is, and I'm hoping that will happen with this posting.
Warning: triggering the bug takes a few minutes of cpu (depending on
your machine). I have been unable to strip down the procedures without
making the bug go away.
Appended to this posting are two routines: VarySpaceGridExplicit[] and
Explicit[]. Explicit[] solves the 1-dimensional heat equation using
an explicit algorithm, and VarySpaceGridExplicit[] investigates the
accuracy of Explicit[] for varying space grids. Sorry for their length,
but as mentioned stripping them down makes the bug disappear.
The following commands work:
VarySpaceGridExplicit[0.05], VarySpaceGridExplicit[0.10],
VarySpaceGridExplicit[0.15], VarySpaceGridExplicit[0.20] and
VarySpaceGridExplicit[0.25].
However, the following command:
Table[VarySpaceGridExplicit[time], {time, 0.05, 0.25, 0.05} ]
does not, giving nonsense results for time == 0.15.
Implementation details:
Mathematica version: 2.2
Platform: Hewlett-Packard 9000/750 running HP-UX 9.01
Any insight will be greatly appreciated.
---
David Harrison | "Music is a hidden practice
Dept. of Physics, Univ. of Toronto | of the soul, that does not
Inet: harrison at faraday.physics.utoronto.ca | know that it is doing
Tel: 416-978-2977 Fax: 416-978-5848 | mathematics." -- Leibniz
--------------------------------------------------------------------------
VarySpaceGridExplicit[time_] := Module[
{ dt, n, max, result, uexact },
uexact[t_, 0.5] := 0.228 /; t == 0.05;
uexact[t_, 0.5] := 0.526 /; t == 0.10;
uexact[t_, 0.5] := 0.710 /; t == 0.15;
uexact[t_, 0.5] := 0.823 /; t == 0.20;
uexact[t_, 0.5] := 0.892 /; t == 0.25;
dt = 0.0001;
max = Round[ time/dt ];
(*
* Explicit[] will return only two solns,
* for t = 0 and t = 'time'. Rather than
* dig out the soln corresponding to x =
* 0.5, exploit the fact that it will be
* the minimum value.
*)
result = Table[ {1/n,
Min[ Last[
Explicit[ 0,1,1,dt,1/n,max, time]
]] - uexact[time, 0.5]
}, {n,50,10,-8}
];
result
] /; (time == 0.05 || time == 0.10 ||
time == 0.15 || time == 0.20 ||
time == 0.25)
(*
* Note that the notation used below for variables follows Golub
* & Ortega, "Scientific Computing and Differential Equations",
* Chapter 8, fairly closely.
*
* We check the input arguments to make sure that they make sense:
* + 1 / dx must be an integer.
* + the time step must divide evenly into the sample times.
*)
Explicit[ inittemp_,
alpha_?NumberQ,
beta_?NumberQ,
dt_?NumberQ,
dx_?NumberQ,
max_Integer,
samplet_?NumberQ ] /;
(IntegerQ[Rationalize[1/dx]] &&
IntegerQ[ Rationalize[ samplet / dt]]) :=
Module[
{
u, (* the current solution *)
mu, (* dt / dx^2 *)
n, (* the number of internal points in the grid *)
iter, (* iteration counter *)
samples (* the list of sampled solutions *)
},
(*
* Set n, the number of internal space grid points. Note
* that although most texts write of n running from 0
* through n + 1, the structure of vectors in MMA means
* we will have to consider n running from 1 through n + 2.
*)
n = Rationalize[ 1/dx ];
n -= 1;
mu = dt / (dx^2);
(* Set initial temperatures *)
If[NumberQ[inittemp],
u = Table[inittemp, {i,0,n + 1}],
(* else *)
u = Table[N[inittemp[i dx]], {i,0,n + 1}]
];
(* Now the boundary conditions *)
u[[1]] = alpha;
u[[n+2]] = beta;
iter = 1;
samples = {u};
While[iter <= max,
iter++;
(*
* The following code is highly functional. Coding
* this procedurally so it looks very much like Eqn.
* 8.2.28 of Golub & Ortega makes Explicit[] slower
* by a factor of about 3.
*)
u = u +
Insert[
Apply[ mu * (#1 + -2 * #2 + #3) & ,
Partition[ u, 3, 1], 1 ],
0, {{1}, {-1}}
];
(* BUG in Mod[]. Test should be "== 0" *)
If[Mod[iter*dt, samplet] < 0.0000000001,
AppendTo[samples, u]
];
];
samples
]
--
David Harrison | "The senses do not lie, only
Dept. of Physics, Univ. of Toronto | they do not tell the truth."
Inet: harrison at faraday.physics.utoronto.ca | -- Mach
Tel: 416-978-2977 Fax: 416-978-5848 |