Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]]]]
share|improve this question
1  
1. $y_0(x)=1+r$ or $y_0(x)=1+x$? 2. Does $y_n(x)$ depend on $a$ or not? –  Artes Jan 3 at 22:00

1 Answer 1

up vote 4 down vote accepted
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}]]
share|improve this answer
    
@herppup, Thanks, $Y[x_, 0] = 1+r$. How to plot all together in a coordinate axis when $a=1/2, r=2$? –  Angel Jan 4 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. –  Angel Jan 4 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}]] –  Angel Jan 4 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}]] –  heropup Jan 4 at 11:12
    
Thanks a million dear heropup. –  Angel Jan 4 at 11:24

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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