Solving recursions

Question

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?

Answer 1

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]