Quotes, global variables, and ODE solvers.

Computed solution: room temperature in a model of a thermostatically switched air conditioning system

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.

\frac{dT}{dt}=0.2(90-T)+0.2S(T)(60-T)

The idea is to model an air conditioning system that is controlled by a thermostatic switch. In the model above the switch function S(T)=1 when the cooling sytem is running, becomes zero if the system is running and the temperature reaches 67F, becomes 1 if the system is not running and the temperature reaches 73F. The behavior of the room temperature is shown in the figure at the top of this piece.

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)));
Incorrect Computed solution: room temperature in a model of a thermostatically switched air conditioning system — the thermostatic switch in the model isn’t functioning correctly.

13 thoughts on “Quotes, global variables, and ODE solvers.”

  1. 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.

    Like

    1. 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.

      Like

      1. 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.

        Like

  2. 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.

    Like

    1. It appears <pre> and <\pre> work, as well as <code> and <\code>

      So for example

      Here is some preformatted text
      

      And here is some code

      subst(23,'param,form);

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s