Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1995
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1995

[Date Index] [Thread Index] [Author Index]

Search the Archive

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       |


  • Prev by Date: Re: How to write a loop
  • Next by Date: Q: Help processing large arrays?
  • Previous by thread: Re: Problem: Plotting list of {InterpolatingFunction[]}
  • Next by thread: Q: Help processing large arrays?