Picard.sagews

def picard_iteration(f, a, c, N):
    '''
    Computes the N-th Picard iterate for the IVP  

        x' = f(t,x), x(a) = c.

    EXAMPLES:
        sage: var('x t s')
        (x, t, s)
        sage: a = 0; c = 2
        sage: f = lambda t,x: 1-x
        sage: picard_iteration(f, a, c, 0)
        2
         sage: picard_iteration(f, a, c, 1)
        2 - t
        sage: picard_iteration(f, a, c, 2)
        t^2/2 - t + 2
        sage: picard_iteration(f, a, c, 3)
        -t^3/6 + t^2/2 - t + 2
        sage: var('x t s')
        (x, t, s)
        sage: a = 0; c = 2
        sage: f = lambda t,x: (x+t)^2
        sage: picard_iteration(f, a, c, 0)
        2
        sage: picard_iteration(f, a, c, 1)
        t^3/3 + 2*t^2 + 4*t + 2
        sage: picard_iteration(f, a, c, 2)
        t^7/63 + 2*t^6/9 + 22*t^5/15 + 16*t^4/3 + 11*t^3 + 10*t^2 + 4*t + 2

    '''
    if N == 0:
        return c*t**0
    if N == 1:
        #print integral(f(s,c*s**0), s, a, t)
        x0 = lambda t: c + integral(f(s,c*s**0), s, a, t)
        return expand(x0(t))
    for i in range(N):
        x_old = lambda s: picard_iteration(f, a, c, N-1).subs(t=s)
        #print x_old(s)
        x0 = lambda t: c + integral(f(s,x_old(s)), s, a, t)
    return expand(x0(t))
v=var('x t s')
a = 0; c = 1; N=5; b=3;
f = lambda t,x: x*t
z=[picard_iteration(f, a, c, i) for i in range(N)]
for i in range(N):
    show(z[i])
from sage.plot.colors import rainbow
c=rainbow(N)
where = [x,-b,b]
p=plot(exp(t^2/2),where,color='gray')                     #Solución exacta.
p+=plot(z[0],where,gridlines=True)
for i in range(1,N):
    p+=plot(z[i],where,color=c[i])

show(p)
    

$\displaystyle 1$

$\displaystyle \frac{1}{2} \, t^{2} + 1$

$\displaystyle \frac{1}{8} \, t^{4} + \frac{1}{2} \, t^{2} + 1$

$\displaystyle \frac{1}{48} \, t^{6} + \frac{1}{8} \, t^{4} + \frac{1}{2} \, t^{2} + 1$

$\displaystyle \frac{1}{384} \, t^{8} + \frac{1}{48} \, t^{6} + \frac{1}{8} \, t^{4} + \frac{1}{2} \, t^{2} + 1$