# How to deal with too much recursion

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.;
heatContent = temp*joulesToHeatWater;
outgoingFlux = emissivity *boltzmanConstant*temp^4;
flux = (incomingFlux - outgoingFlux)*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

• 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. :) – rm -rf 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

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, flux, outgoingFlux, temp}, 5000] // AbsoluteTiming

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

• 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