Revision history [back]

Hello, @sha! I am quite fascinated by @juanjo's answer. I suggest you to take a careful look at his reasoning and derivations. However, I would also like to give you an answer more similar to your original formulation.

But first, I would like to point out that I separated every term in your definition of u[m], and I experimented with each one of them. I found out that the computation stops after the first iteration do to the integral

$$\int_0^t(t-\tau)^{\eta-\alpha-1}u_{m-1}(\tau)P(\tau)\;d\tau$$

The problem is that $\eta-\alpha-1=-0.9$, and Sage seems to be unwilling to compute the integral in the second iteration in that case. (I really don't know why!) I tried to mitigate this problem by computing the indefinite integral first, and that helped a little bit. That's the reason why I used the aux (the indefinite integral) and aux1 (the definite integral) temporary variables in the code that I sent you in my previous comment. However, it didn't help much.

Once again, I don't know why, but using rational numbers instead of floating point numbers did work. (This is also pointed out by @juanjo's answer.) As a consequence, I have replaced the definitions of $\alpha$ and $\eta$ using rational numbers:

alpha = 6/5
eta = 13/10

That tells Sage to use rational numbers, not floating point numbers, and that should avoid the problems with the integral.

Anyway, in the following code, I have deleted the kernel k (since you don't use it in your calculations). I have also deleted your definition of the gamma function, since it is predefined in Sage. On the other hand, notice that

$$\sum_{k=0}^{n-1}\frac{b_k}{k!}t^k$$

and

$$\sum_{k=0}^{n-1}\frac{u_{m-1}^{(k)}}{k!}t^k$$

are the Taylor Polynomials of degree $n-1$ in $t=0$ of the functions $b_ke^t$ and $u_{m-1}(t)$.

Now, I think, @juanjo made a mistake when he said that $\chi$ is unnecessary. I think the problem is that he wrote the original formula (which was written in terms of $m$ and $n$) in terms of just $n$. If I understood correctly, $n$ is a constant and $m$ is your iteration variable (am I right?). If that is the case, the $\chi$ is needed in front of u[m-1], but not in front of your second summation, which is the Taylor Polynomial for $u_{m-1}$, which turns out to be zero for $m=1$, since $u_0(t)=1$. However, notice that $\chi(m)$ is equivalently defined as $\operatorname{sign}(m-1)$, so I used sign(m-1).

Finally, Sage asked me to specify if $t>0$ or $t<0$ (so it can compute the integral). I assumed $t>0$. That's the command assume(t>0).

If I did everything right, the following code should work:

from __future__ import print_function

NUMBER_OF_ITERATIONS = 4

alpha = 6/5
n = 2
eta = 13/10
lamda = 1
h = 1
T = 1

var('h t x s j tau')
U(t) = 1/3*t^3+1
P(t) = 2*t-lamda*(-t^4/24+t^5/15-t^7/126+t^8/72)
b_k = 1

u = vector(SR, NUMBER_OF_ITERATIONS)
u[0] = 1
U = u[0]
assume(t>0)
for m in range(1, NUMBER_OF_ITERATIONS):
    u[m] = b_k*taylor(e^t, t, 0, n-1) + sign(m-1)*u[m-1] - taylor(u[m-1], t, 0, n-1) + h/gamma(eta-alpha)*integrate((t-tau)^(eta-alpha-1)*(u[m-1].subs(t=tau))*P(tau),tau, t, 0); print(u[m])
    U += u[m]

forget()
show(U(t))
print('U(0.1) =', N(U(t=0.1,h=1)))
plot(U(t).subs(h=1), t, 0, 1)

I hope this helps! As always, if the answer needs to be modified or re-explained you can ask as a comment or send me a mail.

Hello, @sha! I am quite fascinated by @juanjo's answer. I suggest you to take a careful look at his reasoning and derivations. However, I would also like to give you an try to answer more similar to your original formulation.post, since I got different results (and formulation, see below), and I have some tiny improvements.

But first, I would like to point out that I separated every term in your definition of u[m], and I experimented with each one of them. I found out that the computation stops after the first iteration do to the integral

$$\int_0^t(t-\tau)^{\eta-\alpha-1}u_{m-1}(\tau)P(\tau)\;d\tau$$

The problem is that $\eta-\alpha-1=-0.9$, and Sage seems to be unwilling to compute the integral in the second iteration in that case. (I really don't know why!) I tried to mitigate this problem by computing the indefinite integral first, and that helped a little bit. That's the reason why I used the aux (the indefinite integral) and aux1 (the definite integral) temporary variables in the code that I sent you in my previous comment. However, it didn't help much.

Once again, I don't know why, but using rational numbers instead of floating point numbers did work. (This is also pointed out by @juanjo's answer.) As a consequence, I have replaced the definitions of $\alpha$ and $\eta$ using rational numbers:

alpha = 6/5
eta = 13/10

That tells Sage to use rational numbers, not floating point numbers, and that should avoid the problems with the integral.

Anyway, in the following code, I have deleted the kernel k (since you don't use it in your calculations). I have also deleted your definition of the gamma function, since it is predefined in Sage. On the other hand, notice that

$$\sum_{k=0}^{n-1}\frac{b_k}{k!}t^k$$

and

$$\sum_{k=0}^{n-1}\frac{u_{m-1}^{(k)}}{k!}t^k$$

are the Taylor Polynomials of degree $n-1$ in $t=0$ of the functions $b_ke^t$ and $u_{m-1}(t)$.

Now, I think, @juanjo made a mistake when he said that $\chi$ is unnecessary. I think the problem is that he wrote the original formula (which was written in terms of $m$ and $n$) in terms of just $n$. If I understood correctly, $n$ is a constant and $m$ is your iteration variable (am I right?). If that is the case, the $\chi$ is needed in front of u[m-1], but not in front of your second summation, which is the Taylor Polynomial for $u_{m-1}$, which turns out to be zero for $m=1$, since $u_0(t)=1$. However, notice that $\chi(m)$ is equivalently defined as $\operatorname{sign}(m-1)$, so I used sign(m-1).

Finally, Sage asked me to specify if $t>0$ or $t<0$ (so it can compute the integral). I assumed $t>0$. That's the command assume(t>0).

If I did everything right, the following code should work:

from __future__ import print_function

NUMBER_OF_ITERATIONS = 4

alpha = 6/5
n = 2
eta = 13/10
lamda = 1
h = 1
T = 1

var('h t x s j tau')
U(t) = 1/3*t^3+1
P(t) = 2*t-lamda*(-t^4/24+t^5/15-t^7/126+t^8/72)
b_k = 1

u = vector(SR, NUMBER_OF_ITERATIONS)
u[0] = 1
U = u[0]
assume(t>0)
for m in range(1, NUMBER_OF_ITERATIONS):
    u[m] = b_k*taylor(e^t, t, 0, n-1) + sign(m-1)*u[m-1] - taylor(u[m-1], t, 0, n-1) + h/gamma(eta-alpha)*integrate((t-tau)^(eta-alpha-1)*(u[m-1].subs(t=tau))*P(tau),tau, t, 0); print(u[m])
    U += u[m]

forget()
show(U(t))
print('U(0.1) =', N(U(t=0.1,h=1)))
plot(U(t).subs(h=1), t, 0, 1)

I hope this helps! As always, if the answer needs to be modified or re-explained you can ask as a comment or send me a mail.

Hello, @sha! I am quite fascinated by @juanjo's answer. I suggest you to take a careful look at his reasoning and derivations. However, I would also like to try to answer your post, since I got different results (and formulation, see below), and I have some tiny improvements.

But first, I would like to point out that I separated every term in your definition of u[m], and I experimented with each one of them. I found out that the computation stops after the first iteration do to the integral

$$\int_0^t(t-\tau)^{\eta-\alpha-1}u_{m-1}(\tau)P(\tau)\;d\tau$$

The problem is that $\eta-\alpha-1=-0.9$, and Sage seems to be unwilling to compute the integral in the second iteration in that case. (I really don't know why!) I tried to mitigate this problem by computing the indefinite integral first, and that helped a little bit. That's the reason why I used the aux (the indefinite integral) and aux1 (the definite integral) temporary variables in the code that I sent you in my previous comment. However, it didn't help much.

Once again, I don't know why, but using rational numbers instead of floating point numbers did work. (This is also pointed out by @juanjo's answer.) As a consequence, I have replaced the definitions of $\alpha$ and $\eta$ using rational numbers:

alpha = 6/5
eta = 13/10

That tells Sage to use rational numbers, not floating point numbers, and that should avoid the problems with the integral.

Anyway, in the following code, I have deleted the kernel k (since you don't use it in your calculations). I have also deleted your definition of the gamma function, since it is predefined in Sage. On the other hand, notice that

$$\sum_{k=0}^{n-1}\frac{b_k}{k!}t^k$$

and

$$\sum_{k=0}^{n-1}\frac{u_{m-1}^{(k)}}{k!}t^k$$

are the Taylor Polynomials of degree $n-1$ in $t=0$ of the functions $b_ke^t$ and $u_{m-1}(t)$.

Now, I think, @juanjo made a mistake when he said that $\chi$ is unnecessary. I think the problem is that he wrote the original formula (which was written in terms of $m$ and $n$) in terms of just $n$. If I understood correctly, $n$ is a constant and $m$ is your iteration variable (am I right?). If that is the case, the $\chi$ is needed in front of u[m-1], but not in front of your second summation, which is the Taylor Polynomial for $u_{m-1}$, $u_{m-1}(t)$ (as I just explained), which turns out to be zero for $m=1$, since $u_0(t)=1$. However, notice that $\chi(m)$ is equivalently defined as $\operatorname{sign}(m-1)$, so I used sign(m-1).

Finally, Sage asked me to specify if $t>0$ or $t<0$ (so it can compute the integral). I assumed $t>0$. That's the command assume(t>0).

If I did everything right, the following code should work:

from __future__ import print_function

NUMBER_OF_ITERATIONS = 4

alpha = 6/5
n = 2
eta = 13/10
lamda = 1
h = 1
T = 1

var('h t x s j tau')
U(t) = 1/3*t^3+1
P(t) = 2*t-lamda*(-t^4/24+t^5/15-t^7/126+t^8/72)
b_k = 1

u = vector(SR, NUMBER_OF_ITERATIONS)
u[0] = 1
U = u[0]
assume(t>0)
for m in range(1, NUMBER_OF_ITERATIONS):
    u[m] = b_k*taylor(e^t, t, 0, n-1) + sign(m-1)*u[m-1] - taylor(u[m-1], t, 0, n-1) + h/gamma(eta-alpha)*integrate((t-tau)^(eta-alpha-1)*(u[m-1].subs(t=tau))*P(tau),tau, t, 0); print(u[m])
    U += u[m]

forget()
show(U(t))
print('U(0.1) =', N(U(t=0.1,h=1)))
plot(U(t).subs(h=1), t, 0, 1)

I hope this helps! As always, if the answer needs to be modified or re-explained you can ask as a comment or send me a mail.