2646.sagews

Solving an equation Picard's way

Introduction: Solving $x^2-x-1=0$

If we are trying to solve an equation like $x^2-x-1=0$, we could try $x^2=x+1$, so that $x=1+\frac{1}{x}$

And we can substitute the value of x again $x=1+\frac{1}{1+\frac{1}{x}}$

And again, and again...

$x+1+\frac{1}{1+\frac{1}{1+\frac{1}{x}}}$

We could, finally, now write

$x=1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{1 + \ldots}}}}$

Picard's way

Hence to solve $x^2-x-1=0$

We first write it as $x=1+\frac{1}{x}$

And then we transform it as $x_1=1$ and $x_{n+1}=1+\frac{1}{x_n}$

And then we get...

$x_1=1$; $x_2=1+\frac{1}{x_1}$; $x_3=1+\frac{1}{x_2}$; $x_4=1+\frac{1}{x_3}$

$x_1=1$; $x_2=1+1/1$; $x_3=1+\frac{1}{1+\frac{1}{1}}$; $x_4=1+\frac{1}{1+\frac{1}{1+\frac{1}{1}}}$

If the sequence $\{x_n\}_{n\in\mathbb{N}}$ converges, an exercise in calculus shows that $\displaystyle{\lim_{x\rightarrow \infty} x_n=x}$ is a solution to the problem.

In this case the sequence converges, other problems like $x=1-x$ lead to non-convergent sequences, and thus they do not lead to solutions...

Now, an interactive demonstration

#So we apply the next one.
#Let's see how it works:
def fromoldtonew(xold):
    return 1+1/xold    
    
@interact
def f(initial=1,n=(1,20,1)):
    html('This program starts with $x_1='+str(initial)+'$ and then iterates applying the function $x_{n+1}=1+\\frac{1}{x_n}$')
    xold=initial
    for i in range(n-1):
        xnew=fromoldtonew(xold)
        xold=xnew
    print html(xold)
    print xold.n()
    print html('Notice that regardless of the initial value, the process approximates the real root $1.618$')

Solving a differential equation

In a similar way, lets try to solve the initial value problem

$\frac{dy}{dx}=y$

$y(0)=1$

This corresponds to $y(x)-y(0)=\int_{0}^{x}y(s)ds$

Solving, and substituing the value of y(0), we get:

$y(x)=\int_{0}^{x}y(s)ds+1$

And we can substitute again...

$y(s)=\int_{0}^{x}\left(\int_{0}^{s}y(s_1)ds_1+1\right)+sds$

Solving a Differential Equation in Picard's way

We can do again the same thing we did for the other equation

$y_1=1$; the constant function 1, which does satisfy that $y(0)=1$

$y_{n+1}=\int_{0}^{x}y_n(s)ds+1$

Now, an interactive demonstration

#Restart the worksheet
s=var('s')
y=var('y')
#So we apply the next one.
#Let's see how it works:    
@interact
def g(fxy='y',x0=0,y0=1,initial=1,n=(1,9,1),assume_x_positive=True):
    fxy=sage_eval(fxy,locals={'x':x,'y':y})
    if assume_x_positive:
        assume(x>0)
    fxy=fxy+x-x+y-y
    r=str(latex(fxy))
    html('We are trying to solve the initial value problem')
    html(r'$\frac{dy}{dx}='+str(latex(fxy))+'$')
    html(r'$y('+str(x0)+')='+str(y0)+'$')
    html(r'Which is equivalent to the problem')
    html(r'$y=\int_{'+str(x0)+'}^x'+r.replace('y','y(x)').replace('x','s')+'ds+'+str(y0))
    r=r.replace('y','y_n(x)')
    r=r.replace('x','s')
    def fromoldtonew(yold):
        return integrate(fxy(y=yold(x=s),x=s),s,x0,x)+y0    
    html('This program starts with $x_1='+str(initial)+'$')
    html(r'and then iterates applying $y_{n+1}=\int_{'+str(x0)+'}^x'+r+'ds+'+str(y0))
    xold=initial+x-x
    print html(r'Initial=')
    show(initial)
    for i in range(n-1):
        xnew=fromoldtonew(xold)
        xold=xnew
        show(xold)
    print html(r'after $%s$ iterates we get'%str(n))
    show(xold)
    print html(r'That is, we get')
    show(xold)

Notice that this, indeed converges to the power series expansion of $e^x$, regardless of the initial functi

html(r'$\frac{1}{x}')

\frac{1}{x}