# 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 $u$ and $H$, 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 $\delta t$ is given as a table

\[Delta]t=Table[i,{i,\[CapitalDelta]t,8,\[CapitalDelta]t}]


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

And $A_R$ is also given in a table via a value $VAR$

VAR=Table[i/4,{i,\[CapitalDelta]t,8,\[CapitalDelta]t}]
AR=1-VAR


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?

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;

f[5]

• Ok, this seems to fix the issue, I'll try to run it, and see what I get. Aug 27, 2015 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). Aug 27, 2015 at 21:07
• You can memoize: tutorial/FunctionsThatRememberValuesTheyHaveFound Aug 27, 2015 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) Aug 27, 2015 at 21:30
• Thanks for the clarification! I'll try to do this :) Aug 28, 2015 at 5:33