How to deal with too much recursion

Question

I have put together a very simple climate model based around four equations which define the state at time t based on time t-1, along with a initial state at time t = 0. The problem arises when I wish to evaluate, say, the temperature at year 2000 (t = 2000) or higher. I hit the recursion limit. I could keep increasing the $RecursionLimit all the time, but I wonder if there is a better way to deal with the problem?

ClearAll[timeStep, waterDepth, gramsPerM3ofWater, joulesPerCalorie, 
  joulesToHeatWater, solarConstant, albedo, boltzmanConstant, 
  emissivity, temp, heatContent, incomingFlux, outgoingFlux, flux];

(*Constants*)
timeStep = 31536000;(*One year in seconds*)
waterDepth = 4000.;
gramsPerM3ofWater = 1000000.;
joulesPerCalorie = 4.186;
joulesToHeatWater = waterDepth*gramsPerM3ofWater*joulesPerCalorie;
solarConstant = 1350.;
albedo = 0.3;
incomingFlux = solarConstant (1 - albedo)/4;
boltzmanConstant = 5.67*^-8;
emissivity = 1.;

(*Initial State*)
temp[0] = 0.;
heatContent[0] = temp[0]*joulesToHeatWater;
outgoingFlux[0] = emissivity *boltzmanConstant*temp[0]^4;
flux[0] = (incomingFlux - outgoingFlux[0])*timeStep;

(*Equations*)
heatContent[t_] := heatContent[t - 1] + flux[t - 1]
temp[t_] := heatContent[t]/joulesToHeat
outgoingFlux[t_] := boltzmanConstant*temp[t]^4
flux[t_] := (incomingFlux - outgoingFlux[t])*timeStep

@belisarius Iteration. That is the answer. [Hoping it's okay to answer a rhetorical question.] — Daniel Lichtblau, Nov 01 '13 at 20:29
@DanielLichtblau Recursion or iteration. That is the question. :) — rm -rf, Nov 01 '13 at 20:37
To understand recursion, you must first understand recursion. — Simon Woods, Nov 01 '13 at 22:57

score 11 · Accepted Answer · answered Nov 01 '13 at 20:27

11

Use Nest, to kill the recursion, as follows:

ClearAll[getNewValues];
getNewValues[{hc_, fl_, out_, temp_}] :=
   Module[{newhc, newfl, newout, newtemp},
      newhc = hc + fl;
      newtemp = newhc/joulesToHeatWater;
      newout = boltzmanConstant*newtemp^4;
      newfl = (incomingFlux - newout)*timeStep;
      {newhc, newfl, newout, newtemp}
   ];

Then, for example for the 5000 steps:

Nest[getNewValues, {heatContent[0], flux[0], outgoingFlux[0], temp[0]}, 5000] // AbsoluteTiming

(* {0.063185, {4.25409*10^12, 0.000237521, 236.25, 254.066}}  *)

answered Nov 01 '13 at 20:27

Leonid Shifrin

114,335
15
329
420

This works nicely, thanks! – Mr Alpha Nov 02 '13 at 17:48
@MrAlpha Was glad to help, and thanks for the accept. This is actually a quite nice application of Nest, I think, and one in which the use of Nest solves what seems to be a rather non-trivial problem in some other approaches. Definitely one of the examples to show the utility of Nest. – Leonid Shifrin Nov 03 '13 at 01:55

How to deal with too much recursion

1 Answers1

Linked