Define recursive function over two indices

Question

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[1]    = 0 ; x[51] = 1;
x[i_]   := x[i] = (i - 1)/50/; i <=51;
y[1]    = y[51] = 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}$.

Answer 1

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}
 ]