1
$\begingroup$

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]]]]
$\endgroup$
  • 1
    $\begingroup$ 1. $y_0(x)=1+r$ or $y_0(x)=1+x$? 2. Does $y_n(x)$ depend on $a$ or not? $\endgroup$ – Artes Jan 3 '14 at 22:00
4
$\begingroup$
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}]]
$\endgroup$
  • $\begingroup$ @herppup, Thanks, $Y[x_, 0] = 1+r$. How to plot all together in a coordinate axis when $a=1/2, r=2$? $\endgroup$ – Angel Jan 4 '14 at 8:56
  • $\begingroup$ 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. $\endgroup$ – Angel Jan 4 '14 at 9:43
  • $\begingroup$ 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}]] $\endgroup$ – Angel Jan 4 '14 at 11:06
  • $\begingroup$ This way of formatting may look nicer: TableForm[Table[{TraditionalForm[Subscript[y, n][x] == Y[x, n, 1/2, r]]}, {n, 1, 10}]] $\endgroup$ – heropup Jan 4 '14 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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