Define recursive function over two indices

I have $x_i=(i-1)/50$ for $i=1,\ldots,51$. Given a function $f(x) = x^2$, I want to find $y_i^k$ given by

\begin{eqnarray} y_i^1 &=& f(x_i) \\ y_i^{k+1} &=& \max\left[\frac{y_{i-1}^{k+1}+y_{i+1}^{k}}{2},\;f(x_i)\right] \end{eqnarray}

This is what I did:

x    = 0 ; x = 1;
x[i_]   := x[i] = (i - 1)/50/; i <=51;
y    = y = 0;
f[x_]   :=x^2;
ff[xi_] :=Composition[f,x][xi];

Then I don't know how to program the relationship between $y_i^k$ and $y_i^{k+1}$.

• I find your question confusing. Perhaps it's just me. – Dr. belisarius Nov 19 '15 at 18:36
• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Nov 19 '15 at 18:36
• What is the value of $y_{1}^{1}$? I ask as it would make use of $y_{0}^{2}$, which is not defined. – Edmund Nov 19 '15 at 19:17

You came close; you can almost transcribe the equations and then let Mathematica do the recursion for you. This isn't necessarily the most efficient way to do things, but it makes up for that in simplicity, and given the size of your problem it's plenty fast enough. First, let's get rid of any stale definitions:

Clear[f, x, y];

Then let's define our $x_i$s, like so:

x[i_Integer] := = i/50.;

It will be important for the $y_i^k$s. Next, we want to define $f(x)$ and $y_i^1$.

f[x_] := x^2;

y[i_Integer, 1] /; 0 <= i < 50 = y[i, 1] = f[x[i]];

The second assignment, to y[i, 1], is how we do memoization in Mathematica, and it saves us from having to recompute things later. I'll be using it throughout the problem, since it's necessary for getting performance that isn't grindingly awful.

y[i_Integer, 1] /; 0 <= i <= 50 := y[i, 1] = f[x[i]];

Now let's impose boundary conditions for $y_i^k = 0$ for $k < 0$ or $50 < k$, which allows for much simpler code:

y[51, _] = y[0, _] = 0.0;

Then we can define our recursion for $y_i^k$:

y[i_Integer?Positive, k_Integer?Positive] := y[i, k] =
Max[0.5*(y[i - 1, k] + y[i + 1, k - 1]), f[i]]

Again, memoization is going to be crucial to making this work. We can now get our answer, which I'll display as a plot, with the initial values for comparison:

ListPlot[Table[{x[i], y[i, #]}, {i, 0, 50}] & /@ {1, 100},
Joined -> True,
AspectRatio -> Automatic,

PlotLegends -> {Subscript[y, i]^1, Subscript[y, i]^100}] If you want to play with larger maximum values, the recursion can become a problem, as Mathematica has a $RecursionLimit that prevents very deep chains of recursive calls. If you compute in stages, perhaps using Table or Do, this will be fine. Do[y[50, i], {i, 200, 1000, 100}]; As you can see, the result converges to$ y_i^k = x_i $for$ k \to \infty \$:

ListPlot[Table[{x[i], y[i, #]}, {i, 0, 50}] & /@ {1, 1000},
Joined -> True,
AspectRatio -> Automatic,

PlotLegends -> {Subscript[y, i]^1, Subscript[y, i]^1000}
] • Thank you for your answer I understand now, I wasn't considering y as a function of two variables ! – Mathelp Nov 19 '15 at 21:13