Solving recursions

This question is tied to my previous question.

I have a few recursion formulas, and I keep getting Recursion depth of 1024 exceeded. error. I need to get the values of uuu and HHH, for each step so that in the end I get the table filled with values. The recursion formulas are:

u3[y_] := -((C’ (u30)^2 (AR[y])^2)/(2*H30)) + 1/2 (((C’)^2 (u30)^4 (AR[y])^4)/(H30)^2 – (4 (u30)^2 (AR[y])^2)/H30 ((C’ \[Lambda] u2[y – 1] Abs[u2[y – 1]])/(2 d) \[Delta]t[y] – C’ u2[y – 1] – H2[y – 1]))^(1/2)

H3[y_] := H2[y-1]-C'(u3[y]-u2[y-1])-(C’\[Lambda] u2[y-1]Abs[u2[y-1]])/(2d) \[Delta]t[y]

u2[y_] := 1/2 ((H10 – H3[y – 1])/C’ + u10 + u3[y – 1] – \[Lambda]/(2 d) (u10 Abs[u10] + u3[y – 1] Abs[u3[y – 1]]) \[Delta]t[y])

H2[y_] := 1/2 (H10 + H3[y – 1] – C’ (u3[y – 1] – u10) – (\[Lambda] C’)/(2 d) (u10 Abs[u10] – u3[y – 1] Abs[u3[y – 1]]) \[Delta]t[y])

The constants that I have are these:

L = 1000;
d = 0.5;
\[CapitalDelta]t = 0.5;
H10 = 100;
H30 = 4.75;
\[Lambda] = 0.012;
c = 1000;
g = 9.81;

u30 = Sqrt[(2 g d (H10 – H30))/(\[Lambda] L)]
u10 = u30;
u20 = u30;
C’ = c/g

H20 = H30 + (\[Lambda]*(L/2)*u30^2)/(2 g d)

The δt\delta t is given as a table


And you can pull the values by putting /. \[Delta]t[1]->\[Delta]t[[1]]

And ARA_R is also given in a table via a value VARVAR


And you can pull the values from here.

Now the recursion formulas are linked. I have the initial conditions, and if I run the first two recursions, and put

u3[1] /. AR[1] -> 0.875 /. \[Delta]t[1] -> 0.5 /. H2[0] -> H20 /. u2[0] -> u20

I get the correct value.

I need to somehow link these recursions. If you put in the first recursion y=3, you get a value that has u2[3] in it. You should be able to get that from the recursion with u2[y], but that recursion depends on the u3[y-1] values (from the previous step).

Can this be done in Mathematica or should I look for something like Python for solving this?



1 Answer


I think your problem is with the specification of the base case. You’re basically trying to do the following, I think:

f[n_] := f[n-1] + 1

f[5] /. f[0] -> 0

Of course, this fails because f[5] tries to evaluate completely before ReplaceAll even sees it.

The solution is to specify a base case in the definition of the function:

f[n_] := f[n-1] + 1
f[0] = 0;




Ok, this seems to fix the issue, I’ll try to run it, and see what I get.
– dingo_d
Aug 27 ’15 at 20:26



Works, but calculating 16 points took ages (for u3). I got to 11 points. Is there a way to optimize the code? Make it run in steps. Like first solve the initial one, assign it to u3[1], then use this in another loop and so on? Because I feel like currently it’s doing all things at once, and that makes it stall. And this is for 16 points (very crude scale). I’d like to have finer division (nicer plot).
– dingo_d
Aug 27 ’15 at 21:07



You can memoize: tutorial/FunctionsThatRememberValuesTheyHaveFound
– Patrick Stevens
Aug 27 ’15 at 21:21



@dingo_d Just to clarify, cause it took me a second, Patrick means that you can paste that in to the help center in MMA to get to that tutorial; alternatively, this is also online here. (+1)
– MarcoB
Aug 27 ’15 at 21:30



Thanks for the clarification! I’ll try to do this 🙂
– dingo_d
Aug 28 ’15 at 5:33