Sharedcloud-examples / sage / Picard.ipynbOpen in CoCalc
Jupyter notebook cloud-examples/sage/Picard.ipynb
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 = 4; N=10; b=2;
f = lambda t,x: sin(t)-2*x; assume(t>0) 
z=[picard_iteration(f, a, c, i) for i in range(N+1)]
for i in range(N+1):
    show(z[i])
from sage.plot.colors import rainbow
c=rainbow(N+1)
where = [x,-b+1,b]
p=plot(-1/5*(cos(t)*e^(2*t) - 2*e^(2*t)*sin(t) - 21)*e^(-2*t),where,color='gray',gridlines=True)                     #Solución exacta.
#p+=plot(z[0],where,gridlines=True)
for i in range(7,N+1):
    p+=plot(z[i],where,color=c[i])

show(p)
44
8tcos(t)+5-8 \, t - \cos\left(t\right) + 5
8t210tcos(t)+2sin(t)+58 \, t^{2} - 10 \, t - \cos\left(t\right) + 2 \, \sin\left(t\right) + 5
163t3+10t210t+3cos(t)+2sin(t)+1-\frac{16}{3} \, t^{3} + 10 \, t^{2} - 10 \, t + 3 \, \cos\left(t\right) + 2 \, \sin\left(t\right) + 1
83t4203t3+10t22t+3cos(t)6sin(t)+1\frac{8}{3} \, t^{4} - \frac{20}{3} \, t^{3} + 10 \, t^{2} - 2 \, t + 3 \, \cos\left(t\right) - 6 \, \sin\left(t\right) + 1
1615t5+103t4203t3+2t22t13cos(t)6sin(t)+17-\frac{16}{15} \, t^{5} + \frac{10}{3} \, t^{4} - \frac{20}{3} \, t^{3} + 2 \, t^{2} - 2 \, t - 13 \, \cos\left(t\right) - 6 \, \sin\left(t\right) + 17
1645t643t5+103t443t3+2t234t13cos(t)+26sin(t)+17\frac{16}{45} \, t^{6} - \frac{4}{3} \, t^{5} + \frac{10}{3} \, t^{4} - \frac{4}{3} \, t^{3} + 2 \, t^{2} - 34 \, t - 13 \, \cos\left(t\right) + 26 \, \sin\left(t\right) + 17
32315t7+49t643t5+23t443t3+34t234t+51cos(t)+26sin(t)47-\frac{32}{315} \, t^{7} + \frac{4}{9} \, t^{6} - \frac{4}{3} \, t^{5} + \frac{2}{3} \, t^{4} - \frac{4}{3} \, t^{3} + 34 \, t^{2} - 34 \, t + 51 \, \cos\left(t\right) + 26 \, \sin\left(t\right) - 47
8315t8863t7+49t6415t5+23t4683t3+34t2+94t+51cos(t)102sin(t)47\frac{8}{315} \, t^{8} - \frac{8}{63} \, t^{7} + \frac{4}{9} \, t^{6} - \frac{4}{15} \, t^{5} + \frac{2}{3} \, t^{4} - \frac{68}{3} \, t^{3} + 34 \, t^{2} + 94 \, t + 51 \, \cos\left(t\right) - 102 \, \sin\left(t\right) - 47
162835t9+263t8863t7+445t6415t5+343t4683t394t2+94t205cos(t)102sin(t)+209-\frac{16}{2835} \, t^{9} + \frac{2}{63} \, t^{8} - \frac{8}{63} \, t^{7} + \frac{4}{45} \, t^{6} - \frac{4}{15} \, t^{5} + \frac{34}{3} \, t^{4} - \frac{68}{3} \, t^{3} - 94 \, t^{2} + 94 \, t - 205 \, \cos\left(t\right) - 102 \, \sin\left(t\right) + 209
1614175t104567t9+263t88315t7+445t66815t5+343t4+1883t394t2418t205cos(t)+410sin(t)+209\frac{16}{14175} \, t^{10} - \frac{4}{567} \, t^{9} + \frac{2}{63} \, t^{8} - \frac{8}{315} \, t^{7} + \frac{4}{45} \, t^{6} - \frac{68}{15} \, t^{5} + \frac{34}{3} \, t^{4} + \frac{188}{3} \, t^{3} - 94 \, t^{2} - 418 \, t - 205 \, \cos\left(t\right) + 410 \, \sin\left(t\right) + 209
N=50;b=1
from sage.plot.colors import rainbow
c=rainbow(N+1)
where = [x,0,b]
p=plot(x^0,where,color='gray',gridlines=True)
for i in range(1,N+1):
    p+=plot(x^i,where,color=c[i])
show(p)
x = var('x')
y = function('y')(x)
show(desolve(diff(y,x) - exp(x+y), y))
(e(x+y(x))+1)e(y(x))=C-{\left(e^{\left(x + y\left(x\right)\right)} + 1\right)} e^{\left(-y\left(x\right)\right)} = C
x = var('x')
y = function('y')(x)
f = desolve(diff(y,x) -exp(x+y), y, ics=[0,1]); show(f)
(e(x+y(x))+1)e(y(x))=(e+1)e(1)-{\left(e^{\left(x + y\left(x\right)\right)} + 1\right)} e^{\left(-y\left(x\right)\right)} = -{\left(e + 1\right)} e^{\left(-1\right)}
t = var('t')
x = function('x')(t)
f = desolve(diff(x,t) -sin(t) + 2*x, x, ics=[0,4]); f
-1/5*(cos(t)*e^(2*t) - 2*e^(2*t)*sin(t) - 21)*e^(-2*t)
t = var('t')
x = function('x')(t)
f = desolve(diff(x,t) -sin(t) + 2*x, x, ics=[0,4]); show(f)
15(cos(t)e(2t)2e(2t)sin(t)21)e(2t)-\frac{1}{5} \, {\left(\cos\left(t\right) e^{\left(2 \, t\right)} - 2 \, e^{\left(2 \, t\right)} \sin\left(t\right) - 21\right)} e^{\left(-2 \, t\right)}
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)
        assume(s>0) 
        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=2; b=.5;
f = lambda t,x: exp(x+t); assume(t>0) 
z=[picard_iteration(f, a, c, i) for i in range(N+1)]
for i in range(N+1):
    show(z[i])
from sage.plot.colors import rainbow
c=rainbow(N+1)
where = [x,-b,b]
p=plot(-log(abs(-1-exp(-1)+exp(t))),where,color='gray',gridlines=True)                     #Solución exacta.
#p+=plot(z[0],where,gridlines=True)
for i in range(N+1):
    p+=plot(z[i],where,color=c[i])

show(p)
11
e+e(t+1)+1-e + e^{\left(t + 1\right)} + 1
e(e+e(t+1))e^{\left(-e + e^{\left(t + 1\right)}\right)}