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}$.

  • 1
    $\begingroup$ I find your question confusing. Perhaps it's just me. $\endgroup$ Nov 19, 2015 at 18:36
  • $\begingroup$ 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! $\endgroup$ Nov 19, 2015 at 18:36
  • $\begingroup$ What is the value of $y_{1}^{1}$? I ask as it would make use of $y_{0}^{2}$, which is not defined. $\endgroup$
    – Edmund
    Nov 19, 2015 at 19:17

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}

more ys

  • $\begingroup$ Thank you for your answer I understand now, I wasn't considering y as a function of two variables ! $\endgroup$
    – Mathelp
    Nov 19, 2015 at 21:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.