Matt Miller and I have recently encountered the same bug with numeric solutions of initial value problems in Maple 7. This bug has been reported to Maple Support; the incident number is 347182. I am posting this message on MUG as an alert to other Maple users. The bug might appear to be rather specific, but the two of us discovered the same bug within the same 24 hour period. (To be honest, Matt showed me the problem he was having, I diagnosed the problem, then encountered the same situation in the process of updating one of my worksheets later in the day.)
The importance of this error is that it occurs in situations where users are often lazy about checking Maple’s output. An unsuspecting user is not likely to discover this problem without much effort.
The problem seems to arise when the procedure returned by dsolve,numeric is evaluated in
the same command in which it is defined. Here is a simple example to illustrate the
problem:
> restart; > interface( version ); Maple Worksheet Interface, Maple 7.00, SUN SPARC SOLARIS, May 28\ 2001 Build ID 96223 > ode := diff( y(t), t ) = y(t); d ode := -- y(t) = y(t) dt > ic1 := y(0) = 1; > ic2 := y(0) = 2; ic1 := y(0) = 1 ic2 := y(0) = 2 > dsolve( { ode, ic1 }, numeric )(1); # exact value should be exp(1) [t = 1., y(t) = 2.71828131105043] > dsolve( { ode, ic2 }, numeric )(1); # exact value should be 2*exp(1) [t = 1., y(t) = 2.71828131105043] > sol1 := dsolve( { ode, ic1 }, numeric ): > sol2 := dsolve( { ode, ic2 }, numeric ): > sol1(1); [t = 1., y(t) = 2.71828131105041937] > sol2(1); [t = 1., y(t) = 5.43656262948751578]
Another workaround, suggested by Brian Fox in Maple’s Technical Support, is to include
output=listprocedure as an additional argument to dsolve:
> dsolve ({ode,ic1}, y(t), numeric, output=listprocedure)(1); [t(1) = 1, y(t)(1) = 2.71828131105041937] > dsolve ({ode,ic2}, y(t), numeric, output=listprocedure)(1); [t(1) = 1, y(t)(1) = 5.43656262948751490]
Note, however, the unfortunate appearance of the left-hand sides of the equations. The use of this output would require modifications to nonstandard notation in subsequent eval or subs commands.
I have not attempted to uncover the full extent or precise source of this problem. Is the bug-generating syntax used in other Maple commands or packages? I suspect it has something to do with the remember table, but don’t have the time to track this down. (Hint, hint!)
You’ve been warned! (:-))
It is corrected with Maple 8 (U. Klein)
I haven’t looked at your specific problem very closely, but I’ll bet that it is related to this remember table bug that I discovered last week in Maple 6. Under certain circumstances, a returned local procedure has a persistent remember table.
> restart; > kernelopts(version); Maple 6.01, SUN SPARC SOLARIS, June 9 2000 Build ID 79514 > interface(verboseproc=3); > P:= proc(ic) > local f; > f:= proc() 0 end; > f(ic[1]):= f(ic[1])+ic[2]; > f > end proc: > P([0,1])(0); 1 > P([0,1])(0); 2
Unlike your example, in this case the bug still occurs if the returned procedure is assigned (why the difference?).
> g:= P([2,2]): > g(2); 2 > h:= P([2,2]): > h(2); 4 > g:= P([2,2]): > g(2); 6
Clear the table.
> g:= subsop(4= NULL, op(g)); g := proc() 0 end proc > g:= P([2,2]); g := f > g(2); 2 > print(h); proc() 0 end proc # (2) = 2
h has acquired a remember table.
> g:= subsop(4= NULL, op(g)); g := proc() 0 end proc > print(h); proc() 0 end proc
h has lost its table.
> addressof(eval(g)); 3339372 > addressof(eval(h)); 3339372
h and g are really the same thing.
This latter fact is true even if there is no remember table.
> restart; > P:= proc() proc() 0 end end: > g:= P(); > h:= P(); > addressof(eval(g)); 3521784 > addressof(eval(h)); 3521784
The bug can be corrected like this:
> P:= proc(ic) > local f; > f:= proc() 0 end; > f:= subsop(4= NULL, eval(f)); #Clear table > f(ic[1]):= f(ic[1])+ic[2]; > f > end: > P([1,1])(1); 1 > P([1,1])(1); 1
But this attempted correction will not work. (Why?)
> P:= proc(ic) > local f; > f:= proc() option remember, system; 0 end; > gc(); #Attempt to remove tables from "system" procedures > f(ic[1]):= f(ic[1])+ic[2]; > f > end: > P([1,1])(1); 1 > P([1,1])(1); 2
I believe that this subs trick is the "official" way of avoiding these problems. Subs is capable of making a true copy of a procedure (including its remember table – although that’s not relevant here). (Thanks to Robert Israel for teaching me that a long time ago.)
> P:= proc(ic) > local f; > f:= subs(_= 0, proc() _ end); > f(ic[1]):= f(ic[1])+ic[2]; > f > end: > P([0,1])(0); 1 > P([0,1])(0); 1 > g:= P([0,1]): > h:= P([0,1]): > addressof(eval(g)); 3509456 > addressof(eval(h)); 3509584
(continuing my previous article on the same subject...)
I have fixed the bug. Yes, it was related to the ideas I gave in the article above, but it was a little more subtle than described there because it only affected the remember table of the procedure named ‘unknown‘.
Correction of the Bug about Unknown Remember Tables in dsolve/numeric
> restart; > kernelopts(version); Maple 7.00, SUN SPARC SOLARIS, May 28 2001 Build ID 96223 > ode:= diff(y(t), t) = y(t): > dsolve({ode, y(0)=1}, numeric, method= rkf45)(1); [t = 1., y(t) = 2.71828131105043]
The problem is that when a procedure is returned, but not assigned to a name, it is assigned to the name "unknown". Unknown can have a remember table. These remember tables suffer from the same globality properties as the remember tables that I showed in my previous article on this subject.
> T:= op(4, eval(unknown)); T := table(["complex" = true, [ 1..25 1-D Array ] "right_comp_soln_data" = [1, 25, [ Data Type: complex ], [ Storage: rectangular ] [ Order: C_order ] [ 1..25 x 1..1 2-D Array ] [ Data Type: complex ], 4, false, 37, 0, true, 15, 3] [ Storage: rectangular ] [ Order: C_order ] ]) > interface(rtablesize=25); > T["right_comp_soln_data"][3]; [0., .0138810615785643, .0277621231571285, .0416431847356928, .0555242463142571, .104703338600836, .153882430887414, .203061523173993, .252240615460571, .302611511380123, .352982407299675, .403353303219226, .453724199138778, .503996888692608, .554269578246439, .604542267800270, .654814957354100, .697963087684838, .741111218015575, .784259348346312, .827407478677050, .870555609007788, .913703739338525, .956851869669262, 1.]
We see that that is the array of t-values.
> convert(T["right_comp_soln_data"][4], listlist); [[1.], [1.01397785097315], [1.02815108215469], [1.04252242449701], [1.05709464738379], [1.11038122051575], [1.16635381252075], [1.22514781794812], [1.28690559131751], [1.35338861554478], [1.42330612136986], [1.49683553403133], [1.57416362907186], [1.65532416183425], [1.74066902853877], [1.83041395956687], [1.92478603124900], [2.00965481087161], [2.09826560191513], [2.19078339230133], [2.28738057535581], [2.38823708346178], [2.49354052961592], [2.60348698245198], [2.71828131105043]]
And that is the array of y-values. These rather complicated remember tables are only used in the new rkf45 and rosenbrock methods of dsolve/numeric.
To correct the bug, we unassign unknown at the start of dsolve/numeric.
> p:= Patch(`dsolve/numeric`):
The lengthy procedure/module Patch is at the end of this article.
> with(p); [&-, &r, &s, Recompile, c, oo, s] > &-1; 1 userinfo(1,{dsolve, `dsolve/numeric`},'entering'); > 0&r"`Start patch: Carl 12Sep2001`; \ > proc(p) \ > if p then \ > userinfo(1, {dsolve, `dsolve/numeric`}, `Unassigning protected unknown.`); \ > unprotect(':-unknown') \ > fi; > unassign(':-unknown'); \ > if p then protect(':-unknown') fi \ > end \ > (evalb(:-unknown::protected)); \ > `End patch`; \ > "; 1 `Start patch: Carl 12Sep2001`; proc(p) if p then userinfo(1, {dsolve, `dsolve/numeric`}, `Unassigning protected unknown.`); unprotect(':-unknown') fi; unassign(':-unknown'); if p then protect(':-unknown') fi end (evalb(:-unknown::protected)); `End patch`; userinfo(1,{dsolve, `dsolve/numeric`},'entering'); Note the trick above for including *permanent* comments in code. > Recompile();
The Patch/Recompile procedure does not currently install the patched code in a library. That is being worked on. For the time being, you may install it in a library yourself.
Test the patch:
> infolevel[all]:= 1: > dsolve({ode, y(0)=1}, numeric)(1); dsolve/numeric: entering [t = 1., y(t) = 2.71828131105043] > dsolve({ode, y(0)=2}, numeric)(1); dsolve/numeric: entering [t = 1., y(t) = 5.43656262948750] > f:= dsolve({ode, y(0)=1}, numeric); dsolve/numeric: entering > f(1); [t = 1., y(t) = 2.71828131105041937] > f:= dsolve({ode, y(0)=2}, numeric); dsolve/numeric: entering > f(1); [t = 1., y(t) = 5.43656262948751578]
Here is the procedure Patch. This is a work-in-progress that will soon be distributed as a patch-library-maintainence system.
A Patching System for Maple’s Internal Source Code
Author: Carl Devore <devore@math.udel.edu>
12 September 2001
Thanks to Joseph Riel for the idea that it might be possible to recompile the modified output of showstat. Thanks to Robert Israel for the idea that local procedures can be extracted from modules.
> restart; > Patch:= proc(ProcName::evaln, m::`module`) > option `Copyright (C) 2001, Carl James DeVore, Jr. All rights reserved.`; > global _PatchLock; > local M, IsModule, _TheProc, S, `&->`, p, LockCheck, findline, findlines, strip, ReplaceText; > # If this module is the Patcher, then this procedure below is the Snatcher. > # It will retrieve any export or local from any module. Note that in the > # case of exports, this procedure works in cases where ":-" does not. > `&->`:= proc(M::`module`, Item::evaln) > local items, p, stringer, item; > stringer:= proc(item) > local s, p, newp; > s:= sprintf("%A", item); > p:= 1; > do # Might need to remove some module name prefixes > newp:= searchtext(":-", s, p..-1); > if newp=0 then return s[p..-1] else p:= p+newp+1 fi > od > end proc; > items:= [op(3, eval(M))]; #The locals > item:= stringer(Item); > if not member(item, map(stringer, items), p) then > items:= [op(1, eval(M))]; # The exports > if not member(item, map(stringer, items), 'p') then error "%1 not found in module", Item fi > fi; > eval(items[p]) > end proc; > # Find index in S of a line number. > findline:= (n::nonnegint)-> `if`(n=0, 0, searchtext(cat("\n", " " $ 4-length(n), n), S)); > # Find indices into S that delimit a pair of line number. > findlines:= proc(m::nonnegint, n::posint) > local p; > p:= findline(m); > if p=0 and m<>0 then error "linenumber %1 is too large", m fi; > p+1..findline(n)-1 > end proc; > > # Strip replaces \n's with blanks. When a single line is displayed, this > # format is used. That way the arrows in search string markers and compile > # error markers will make sense. > strip:= proc(In::string) > # This proc can be done better in Maple 7 with the StringTools. > local Out; > Out:= In; > do try Out:= ReplaceText(Out, "\n", " ") catch: return Out end od > end proc; > # ReplaceText finds the nth occurence of s1 in s and replaces it with s2. > # n defaults to 1. > # This procedure is a companion to Maple's SearchText. > # The fourth argument is not currently used in this module. > ReplaceText:= proc(s::string, s1::string, s2::string, n::posint) > local i,p,l; > i:= 0; > l:= length(s); > to `if`(nargs=4, n, 1) do > if i=l then error "String not found: %1", s1 fi; > p:= SearchText(s1, s, i+1..-1); > if p=0 then error "String not found: %1", s1 else i:= i+p fi > od; > cat(s[1..i-1], s2, s[i+length(s1)..-1]) > end proc; > IsModule:= evalb(nargs>1); > if IsModule then > M:= m; > if not eval(M)::`module` then error "2nd argument must the name of a module" fi; > _TheProc:= eval(M)&->ProcName; > # At this point, it might me good to look for lexical vars, and return > # an error in that case. > else > M:= ``; > _TheProc:= ProcName > fi; > if not eval(_TheProc)::procedure then > error "First argument must be the NAME of a procedure" > fi; > > # The only data structure in this module is a copy of the procedure > # specifed by P that is put into a string > # and has had \n's (newlines) and statement numbers added by debugopts. > > S:= sprintf("%s", debugopts(procdump= _TheProc)); > # Note the above statement causes "_TheProc" to appear in an assignment > # statement in S. > p:= SearchText("|Calls Seconds Words|\nPROC |", S); > if p>0 and p<SearchText("\n 1", S) then > error "Turn off statement-level tracing before patching." > fi; > > # Since the parse statement operates at the global level, we only allow > # one patch-in-progress at a time. > # I could do some complicated stuff to get around this, but I can't think > # of any good reason why someone would want to do this. A patch is > # considered to be "in progress" if it has not been recompiled. > if assigned(:-_PatchLock) and :-_PatchLock then > WARNING("can only patch one proc at a time. Previous patch will be scrapped.") > fi; > :-_PatchLock:= true; > LockCheck:= proc() > if not :-_PatchLock then error "Already Recompiled. Reload module to patch again." fi > end proc; > module() > export `&-`, `&s`, `&r`, Recompile, oo, c, s; > local findcurrentlineno, i, Relink; > > s:= ``; > c:= ``; > # Find the last line number. Set to oo. > for i do if S[-i]="\n" then try oo:= sscanf(S[-i+1..-i+5], "%d")[]; if oo::posint then break fi catch: end fi od; > findcurrentlineno:= proc(S::string) > local i,n; > for i to length(S) do > if S[i]="\n" then try n:= sscanf(S[i+1..i+5], "%d")[]; if n::posint then return n-1 fi catch: next end fi > od; > oo > end proc; > # The procedure Recompile: > # 1. Replaces all showstat formatting codes and statement numbers by blanks. > # The length is not changed to facilitate returning to the exact point > # that caused on error. The original S is not changed, a copy is made. > # 2. Compiles string S to a procedure using parse. Sometimes debugopts > # creates syntax that parse won't accept. > # The only example I know is "catch NULL:". I attempt to figure out if > # the compilation error is caused by this, and then change these to > # simply "catch :" > # 3. If compilation fails, then the user is returned to the editor with > # an "^" pointing to the position that parse is complaining about. > # 4. If the procedure came from module, then it is relinked to the module. > # > # Note that parse statements operate at the global level. That means that > # the name of the procedure is assigned by the parse statement. > > Recompile:= proc() > local In, Out, B, N, Protected, LE; > LockCheck(); > In:= S; > Out:= ""; > do > B,In:= sscanf(In, "%[^\n]%[^\r]")[]; > Out:= cat(Out, B, " "); > if In="\n" then break fi; > B,N,In:= sscanf(In[2..-1], "%[ ]%[0-9]%[^\r]")[]; > Out:= cat(Out, " " $ length(B)+length(N)) > od; > Protected:= evalb(_TheProc::protected); > if Protected then unprotect(_TheProc) fi; > try > parse(cat(Out, ";"), statement) > catch: > do try Out:= ReplaceText(Out, "catch NULL:", "catch :") catch: break end od; > try > parse(cat(Out, ";"), statement) > catch: > LE:= [lastexception]; > if nops(LE)=4 and LE[-1]::posint then > # Display the compilation error > &-findcurrentlineno(S[LE[-1]..-1]); > printf("\n%s", cat(" " $ LE[-1]-findline(c)-2, "^")); > error LE[2][1..-5], LE[3..-2] > fi; > error LE[2..-1] > end try > finally > if Protected then protect(_TheProc) fi > end try; > if IsModule then > # The _TheProc in the next statement is the global _TheProc. Note > # how this trick allows the assignment to > # any export of any module as long as the name of module is global. > try parse(sprintf("%A:-%A:= _TheProc;", M, ProcName), statement) > catch: Relink(Procname, eval(:-_TheProc), eval(M)) > end > fi; > # If the compilation and relink were successful, then the current patch > # is not "in progress" > :-_PatchLock:= false; > [][] > end proc; > # Relink a procedure local to a module. > Relink:= proc(N::name, P::procedure, M::`module`) > # Waiting for Helmut Kahovec to write this procedure. > error "Relinking local procedures to their parent modules not yet implemented." > end; > > # Search command. Sets the current line and the current search string. > `&s`:= proc(N::{nonnegint,string}, Search::string) > local n,p,pos; > LockCheck(); > if nargs=0 then #repeat current search, advancing one line > if not s::string then > error "Use &s\"string\"; or s:= \"string\"; to set the current search." > fi; > n:= `if`(c::nonnegint, c, 0) > elif nargs=1 then > if not N::string then error "String expected: %1", N fi; > s:= N; > n:= `if`(c::nonnegint, c, 0) > else > if N::string then error "Line number expected for left argument." fi; > n:= N; > s:= Search > fi; > if n=oo then n:= 0 fi; #Recycle the search > #If current line is 0, check it; otherwise advance 1 line before searching in the no-arg call. > p:= `if`(n<1, 1, findline(`if`(nargs=0, n+1, n))); > pos:= SearchText(s, S, p..-1); > if pos=0 then error "Not found." fi; > &-findcurrentlineno(S[p+pos-1..-1]); > printf("\n%s", cat(" " $ p+pos-findline(c)-2, "^" $ length(s))) > end proc; > # Line display command. If used with a single argument, it also sets the > # current line. > `&-`:= proc(m::nonnegint, n::posint) > local s; > LockCheck(); > s:= S[findlines(m, `if`(nargs=1, m+1, min(10^5, n+1)))]; > if nargs=1 then s:= strip(s) fi; > printf("%s", s); > c:= `if`(nargs=1, m, ``); > [][] > end proc; > # Replace command. Only works on the current line. > `&r`:= proc(N::{string, identical(0), identical(1)}, Replace::string) > local p,q,replace,n,new; > LockCheck(); > if nargs=1 and not s::string then > error "Use &s\"string\"; or s:= \"string\"; to set current search." > fi; > if not c::nonnegint then error "Use &-n; or c:=n; to set line number to n before replacing." fi; > if nargs=1 then replace:= N; n:= `` else n:= N; replace:= Replace fi; > p,q:= op(findlines(c, c+1)); > if n="" then new:= replace #replace whole line > elif n=`` then new:= ReplaceText(S[p+6..q], s, replace) #replace search string > elif n=0 then new:= cat(replace, S[p+6..q]) #prepend to line > elif n=1 then new:= cat(S[p+6..q], replace) #append to line > else error "Left argument can only be \"\", 0, or 1, but received %1", n > fi; > S:= cat(S[1..p+5], new, `if`(c=oo, "", S[q+1..-1])); > &-c > end proc > end module > end proc:
A bit of help about the various editting commands:
Let’s test it on a procedure local to a module. Note, thus, that we are looking at code that is ordinarily difficult to look at.
> p:= Patch(`Map/Internal`, eval(LinearAlgebra)): > with(p); [&-, &r, &s, Recompile, c, oo, s]
List the entire proc. Lines start at 0 (the procedure header), and go to oo.
> 0&-oo; _TheProc := proc(MV, f) local x, i, j, rows, cols; 1 if type(MV,{'Vector', 'Matrix'}) then 2 if member(rtable_options(MV,storage),{sparse, sparse[upper], sparse[lower]}) then 3 for x in rtable_elems(MV) do 4 try 5 MV[lhs(x)] := f(rhs(x),args[3 .. -1]) catch NULL: 6 NULL end try end do elif type(MV,'Vector') then 7 for i to op(1,MV) do 8 try 9 MV[i] := f(MV[i],args[3 .. -1]) catch NULL: 10 NULL end try end do else 11 rows, cols := op(1,MV); 12 for i to rows do 13 for j to cols do 14 try 15 MV[i,j] := f(MV[i,j],args[3 .. -1]) catch NULL: 16 NULL end try end do end do end if; 17 MV else 18 f(MV,args[3 .. -1]) end if end proc
Note that the line numbers are actually the statement numbers as reported by showstat. These numbers are fixed throughtout a single editting session.
c refers to the current line number.
> c; There is no current line.
Search for 1st occurence of "MV"
> &s"MV"; _TheProc := proc(MV, f) local x, i, j, rows, cols; ^^ Check current line and search. > c; 0 > s; "MV"
Set current line to 3 and display.
> &-3; 3 for x in rtable_elems(MV) do
Replace entire line.
> ""&r" `junk`;"; 3 `junk`;
Search for "cols" starting on the this line.
> &s"cols"; 11 rows, cols := op(1,MV); ^^^^
Put something at the beginning of the line.
> 0&r"junk2; "; 11 junk2; rows, cols := op(1,MV);
Search for "MV", starting at line 7.
> 7&s"MV"; 7 for i to op(1,MV) do
^^ Put something at the end of the line.
> 1&r"junk 3"; 7 for i to op(1,MV) dojunk 3
Set current line to 0.
> &-0; _TheProc := proc(MV, f) local x, i, j, rows, cols;
Replace current search string in current line.
> &r"junk 4"; _TheProc := proc(junk 4, f) local x, i, j, rows, cols;
Search again for the current search string.
> &s(); 1 if type(MV,{'Vector', 'Matrix'}) then ^^ > Recompile(); _TheProc := proc(junk 4, f) local x, i, j, rows, cols; ^ Error, (in Recompile) incorrect syntax in parse: unexpected number
Note that the position of the compilation error is indicated.
Scrap this editting session. Let’s write our own module to edit.
> mymodule:= module() > export myproc; > myproc:= proc(P) try print(`Hello, world.`) catch: end end > end; mymodule := module _m649292 () export myproc; end module > p:= Patch(myproc, mymodule): Warning, can only patch one proc at a time. Previous patch will be scrapped. > with(p); Warning, these names have been rebound: &-, &r, &s, Recompile, c, oo, s [&-, &r, &s, Recompile, c, oo, s] > 0&-oo; _TheProc := proc(P) 1 try 2 print(`Hello, world.`) catch NULL: 3 NULL end try end proc
Note the "catch NULL:" This is incorrect syntax generated by showstat. The Recompile
handles it gracefully.
> &s"Hello"; 2 print(`Hello, world.`) catch NULL: ^^^^^ > &r"Goodbye"; 2 print(`Goodbye, world.`) catch NULL: > Recompile(); > mymodule:-myproc(); Goodbye, world.
That proves that the module was actually changed. Also note that the procedure was automatically reinserted into its parent module. The full version of Patch will also reinsert local procedures into their parent modules.
Test the protection feature:
> P:= proc() "junk" end; > protect(P); > evalb(P::protected); true > p:= Patch(P): > with(p); Warning, these names have been rebound: &-, &r, &s, Recompile, c, oo, s [&-, &r, &s, Recompile, c, oo, s] > &s"junk"; 1 "junk" end proc ^^^^ > &r"new junk"; 1 "new junk" end proc > Recompile(); > P(); "new junk" > evalb(P::protected); true