Every winter term in my introductory ordinary differential equations class, I ask each of my students to choose from the scientific literature a paper that puts ODE models at the center of the work, and to reproduce the author’s results. For the past 5 years we’ve used Maxima for that purpose, and throughout the course for both symbolic and numerical tasks.
A perennial favorite is “Temperature Models for Ware Hall” by J. K. Denny and C. A. Yackel, The College Mathematics Journal, Vol. 35, No. 3 (May, 2004), pp. 162-170 http://www.jstor.org/stable/4146891 .
We have a first order ODE that combines two Newton’s law of cooling processes: the room is in thermal contact with the hot outside air at 90F and with cooled air at 60F coming from the air conditioner.
The idea is to model an air conditioning system that is controlled by a thermostatic switch. In the model above the switch function 
There might be another way to model this (I’d be delighted to learn what that is!) but here we treat the thermostat with a global variable that can be read and changed within a numerical solver as part of the right-hand-side function.
The question in this post is: how to implement a global flag like that in Maxima? Below are four variants involving various combinations of the Maxima quote operator. The first two give the hoped-for behavior and were used to produce the figure above. The last two don’t seem to have access to the global flag and the cooling system runs without ceasing. I include a graph of that incorrect computed solution at the bottom of this post.
The Details
- In the first case, the switch function is included in the right-hand-side function with a single quote.
- In the second, the switch in included with quote-quote, but the order in which the right hand side and switch is important: the switch function must be defined AFTER the right hand side function in which it is called
- In the third, the switch in included with quote-quote, but the order from above is reversed and the global flag seems undetected.
- In the fourth, the switch is included without quote or quote-quote and the behavior is same as quote-quote in the second and third cases.
I don’t know why each works or doesn’t! I’m hoping you do and can tell me 🙂 Thanks in advance.
/*This works as intended with single-quote on ThermSwitch
ThermSwitch can be defined either before or after ff */
kill(all);
globalAC:1;
ThermSwitch(T):=block(
if T<67 then globalAC:0
else if T>73 then globalAC:1,
return(globalAC))$
ff:.05*(90-T) + .2*'ThermSwitch(T)*(60-T);
s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));/*This works as intended with double-quote on ThermSwitch
ThermSwitch defined AFTER ff */
kill(all);
globalAC:1;
ff:.05*(90-T) + .2*''ThermSwitch(T)*(60-T);
ThermSwitch(T):=block(
if T<67 then globalAC:0
else if T>73 then globalAC:1,
return(globalAC))$
s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));/*This does NOT work as intended with double-quote on ThermSwitch
ThermSwitch defined BEFORE ff */
kill(all);
globalAC:1;
ThermSwitch(T):=block(
if T<67 then globalAC:0
else if T>73 then globalAC:1,
return(globalAC))$
ff:.05*(90-T) + .2*''ThermSwitch(T)*(60-T);
s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));/*This does NOT work as intended with unquoted ThermSwitch before ff
it does work as intended with unquoted ThermSwitch after ff */
kill(all);
globalAC:1;
ThermSwitch(T):=block(
if T<67 then globalAC:0
else if T>73 then globalAC:1,
return(globalAC))$
ff:.05*(90-T) + .2*ThermSwitch(T)*(60-T);
s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));
The MaximaList 
Well, first of all, you can eliminate Quote-Quote from consideration. Quote-Quote is NOT a quoting operator — it substitutes in a value early, almost the OPPOSITE of quotation. The confusion is compounded by some documentation about it being incorrect.
Without looking at your problem in detail (or trying out my solution), I think what you want is ff: ‘(.05*(90-T) + .2*ThermSwitch(T)*(60-T)). This prevents the contents of the parentheses from being evaluated.
You can also simplify your code by writing:
ThermSwitch(T):= globalAC: (if T73 then 1 else globalAC)$
It also seems like a bad practice a priori to pass a function which depends on modifying a global variable to the diffeq solver. It apparently works find in this case, but a more sophisticated solver might not evaluate the points in temporal order. But I’m afraid I don’t have any positive suggestions o how to do this better.
LikeLike
Thank you Stavros!
Can I please trouble a bit further about this?
I think I see that we want to prevent the function from being evaluated initially so that the rule of the function remains and not just a value. I’m imagining that eventually the solver rk() does need to evaluate the function. rk() maybe does something special for that, but in general what maxima mechanism is used to evaluate a quoted expression at a value of the argument?
grateful for your time and consideration.
LikeLike
To evaluate an expression, use `ev`:
> `param: 23`
> `formula: ‘(param*2)` → `param*2` —param is not evaluated
> `ev(formula)` → `46`
You can also use subst:
> `subst(23,’param,formula)` → `46` —note the quoting, otherwise it’ll try to substitute 23 for 23!
Internally, `rk()` actually uses a function called `coerce-float-fun`, which tries to be clever about compiling the expression for efficient floating-point evaluation.
LikeLike
Both of your scripts work if you change the ” to ‘. Using Mr. Macrackis’ quoting structure, ff: ‘(.05*(90-T) + .2*ThermSwitch(T)*(60-T)), also works in either of your scripts.
LikeLike
Thank you!
LikeLike
How do I format comments? Apparently Markdown doesn’t work here.
LikeLike
It appears <pre> and <\pre> work, as well as <code> and <\code>
So for example
And here is some code
subst(23,'param,form);
LikeLike
I’m also wondering if the WordPress LaTeX interpreter works in comments:
LikeLiked by 1 person
What was the markup for that,
x^2?
And how do you edit and delete comments…?
LikeLike
Hmm. I was trying to show my markup. Does ampersand work? <math>?
This must be documented somewhere?
LikeLike
LikeLike
Umm, I can’t see what markup you used to get that. I tried <math> x^2 </math> doesn’t work for me: x^2.
LikeLike
I don’t know how to quote this literally, so:
$, followed by the word latex, x^2, followed by a closing $
I wrote
a little post about this a few years ago
LikeLike