# Defining a recursive integral sequence

I have a recursive integral sequence as follows: $$y_0(x)=1+r,$$ $$y_1(x)=-\frac{1}{\Gamma(a)}\int_0^x (x-t)^{a-1} y_0(t) dt,$$ $$\vdots$$ $$y_n(x)=-\frac{1}{\Gamma(a)}\int_0^x (x-t)^{a-1} y_{n-1}(t) dt,\quad 0<a<1.$$ I wrote the following code but the output answer isn't correct. I know this problem is caused by the change of variables($t‎\leftrightarrow‎ x$ ) in $y_n(x)$. How to modify this problem?

Y[0, x_] = 1 + r;
Y[k_, x_] :=
Y[k, x] =FullSimplify[-(1/Gamma[a])*Integrate[(x - t)^(a - 1)*(Y[k - 1, t]), {t, 0, x}]]
For[i = 0, i <= 10, i++, Print["Y", i, "[x]=", Simplify[Y[i, x]]]]

• 1. $y_0(x)=1+r$ or $y_0(x)=1+x$? 2. Does $y_n(x)$ depend on $a$ or not? Jan 3 '14 at 22:00

Y[x_, n_] := Y[x, n] = -Simplify[
Integrate[(x - t)^(a - 1) (Y[x, n - 1] /. x -> t), {t, 0, x},
Assumptions -> 0 < a < 1, GenerateConditions -> False]/Gamma[a]]


If we set Y[x_, 0] = 1+r then this yields the result $$y_n(x) = \frac{(-1)^n (1+r) x^{na}}{\Gamma(1+na)}.$$ Is this the answer you are looking for? In the case where you intended $y_0(x) = 1 + x$ instead of $1+r$, the result would be $$y_n(x) = \frac{(-1)^n (1+na+x)x^{na}}{\Gamma(2+na)}.$$ You can also do the whole thing as a single NestList command:

NestList[-Simplify[
Integrate[(x - t)^(a - 1) (# /. x -> t), {t, 0, x},
Assumptions -> 0 < a < 1, GenerateConditions -> False]/
Gamma[a]] &, 1 + r, 10]


Once we know the form of $y_n(x)$, we can just create an explicit definition so that it is not necessary to do the computation:

Y[x_, n_, a_, r_] := (-1)^n (1 + r) x^(n a)/Gamma[1 + n a]


Then to get a plot for $a = 1/2$, $r = 2$, $n = 1, 2, \ldots, 10$, we just do something like

Plot[Evaluate[Table[Y[x, n, 1/2, 2], {n, 1, 10}]], {x, 0, 3}]


Or if you want to be fancy and make an interactive plot,

Manipulate[Plot[Y[x, n, 1/2, 2], {x, 0, 3}, PlotRange -> {-6, 6}], {n, 1, 10, 1}]


And for a two-column set of formulas for $a = 1/2$,

TableForm[Table[TraditionalForm /@ {Y[x, n, a, r], Y[x, n, 1/2, r]}, {n, 1, 10}]]


If you want TeX output, instead use

TeXForm[Table[{Y[x, n, a, r], Y[x, n, 1/2, r]}, {n, 1, 10}]]

• @herppup, Thanks, $Y[x_, 0] = 1+r$. How to plot all together in a coordinate axis when $a=1/2, r=2$? Jan 4 '14 at 8:56
• How to print Y[x, n,a] and Y[x, n,a=1/2] in two columns facing each other s.t n=1,...,10. Jan 4 '14 at 9:43
• Excellent. Sorry, I wanna appear " $Y_n[x]=$ " in the Output. TableForm[Table[TraditionalForm /@ {Print["Y", n, "[x]="]Y[x, n, a, r],Print["Y", n, "[x]="] Y[x, n, 1/2, r]}, {n, 1, 10}]] Jan 4 '14 at 11:06
• This way of formatting may look nicer: TableForm[Table[{TraditionalForm[Subscript[y, n][x] == Y[x, n, 1/2, r]]}, {n, 1, 10}]] Jan 4 '14 at 11:12