Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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];

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;

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
share|improve this question
Recursion or iteration. That is the question. :) – Dr. belisarius Nov 1 '13 at 19:10
@belisarius Iteration. That is the answer. [Hoping it's okay to answer a rhetorical question.] – Daniel Lichtblau Nov 1 '13 at 20:29
@DanielLichtblau Recursion or iteration. That is the question. :) – R. M. Nov 1 '13 at 20:37
Apparently it's already devolved into reiteration. – bill s Nov 1 '13 at 21:25
To understand recursion, you must first understand recursion. – Simon Woods Nov 1 '13 at 22:57
up vote 11 down vote accepted

Use Nest, to kill the recursion, as follows:

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}}  *)
share|improve this answer
This works nicely, thanks! – Mr Alpha Nov 2 '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 3 '13 at 1:55

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.