Differential Equations Introduction

Goals

Solutions: Theoretical computation
Visualization
Properties of Equations: Stability, time behavior
Machine learning: Iuput equations, output numerical scheme

Consider \begin{eqnarray} \frac{d y}{d x} &=& x-y,\\ y(x_0)&=&y_0 \end{eqnarray}

Solution is sum of general solution, $y_g$ which satisfies $y_g'+y_g=0$, and special solution, $y_s$ which satisfies $y_s'+y_s=x$ since \begin{eqnarray} y'&=&(y_g+y_s)' \\ &= &y_g'+y_s'=-y_g+x-y_s=x-y \end{eqnarray}

$\frac{d}{d x} =D$ implies $ (D+1)y=x$.
$y_g$ satisfies $ (D+1)y_g=0$. \begin{eqnarray} (D+1)y_g&=&0 \\ (r+1)=0 &\to &r=-1 \\ y_g &= & C_1 \exp(-1\times x)=C_1 \exp(-x) \end{eqnarray}
$y_s$ can be calculated by the follows: \begin{eqnarray} (D+1)y_s&=&x \\ y_s &= &\frac{1}{1+D} x \\ &= & (1-D+D^2-D^3+\cdots) x \\ &= & x-1 \end{eqnarray}
$y(x_0)=y_0$ implies $y(x)=e^{x_0} (y_0+1-x_0) e^{-x}+x-1$: \begin{eqnarray} y(x_0)_s&=& y_0 \
```
  &=& y_g(x_0)+y_s(x_0) \\
  &=& C_1 \exp(-x_0)+x_0-1 \\
```
&\to & C_1=\exp(x_0)(y_0+1-x_0) | \end{eqnarray}

Conclusion

$$ y(x)=e^{x_0} (y_0+1-x_0) e^{-x}+x-1$$



In [1]:

    
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline



In [2]:

    
def f(x,y):                # DE is y'=f(x,y) for 
   return x-y              # this function f

def yexact(x,x0,y0):       # exact solution y(x) that satisfies y(x0)=y0
   y=y0*exp(-x + x0) + (x*exp(x) - x0*exp(x0) - exp(x) + exp(x0))*exp(-x)
   return y



In [13]:

    
fig=plt.figure(figsize=(5,5))          # fig is figure variable
ax=fig.add_subplot(1,1,1) # in 1x1 array of subplots, ax is subplot 1

#					  I. Slope Field

xmesh=linspace(-4.,4.,9); ymesh=linspace(-4.,4.,9) # x and y mesh values

#                                       plot mini-tangent for
#                                       y'=f(x,y) at each mesh point
for xp in xmesh:
   for yp in ymesh:
      m=f(xp,yp)
      h=0.25/sqrt(1.+m**2)
      #ax.plot([xp-h,xp+h],[yp-m*h,yp+m*h],'b')
      ax.arrow(xp,yp,2*h,2*m*h, head_width=0.2, head_length=0.2, fc='r', ec='r') 
                          
#					II. Dot at given (x0,y0) and solution curve through (x0,y0)
x0=-4.; y0=3.5
ax.plot([x0],[y0],'mo') # 'm'agenta d'o't
X=linspace(-4.,4.,101)  # X[0]=-4., .., X[100]=4.
Y=yexact(X,x0,y0)
ax.plot(X,Y,'m',linewidth=2) # lines joining points (X[i],Y[i]), 'm'agenta

ax.set_xlabel('x'); ax.set_ylabel('y')
ax.grid(True)           # add a grid to the plot
ax.set_aspect(1.)       # require that (scaling factor for y)/(scaling factor for x), is 1.
plt.title("Slope Field for y'=x-y and\nsolution curve through ({0:g},{1:g})".format(x0,y0))
#plt.savefig("306ch1_slope_field_soln_1.3.5.png")
#plt.show() # to display immediately









    Out[13]:





<matplotlib.text.Text at 0x107929610>



In [ ]:

Symbol DE solver

Solve $y'=x-y, y(0)=y_0$



In [33]:

    
from sympy import Symbol , dsolve , Function , Derivative , Eq, lambdify, symbols,solve



In [34]:

    
x=Symbol('x')
x0,y0=symbols('x0 y0')
y=Function('y')



In [35]:

    
f_=Derivative(y(x),x)-x+y(x)
sol=dsolve(f_,y(x))



In [36]:

    
print sol









    



y(x) == (C1 + (x - 1)*exp(x))*exp(-x)



In [52]:

    
yy=sol.rhs
C1=Symbol('C1')



In [53]:

    
z= lambdify(x,sol.rhs-y0)
z(0)









    Out[53]:





1.0*C1 - y0 - 1.0



In [55]:

    
c=solve(z(0),C1)



In [56]:

    
ysol=yy.subs({C1:c[0]})
print ysol









    



(y0 + (x - 1)*exp(x) + 1.0)*exp(-x)



In [ ]:

Steps for Scheme Auto-generated

SymPy -> theoretical solution
NumPy -> numerical computation
Matplotlib -> visualization



In [5]:

    
def ode1(f,x0,y0):
    sol=dsolve(f,y(x))
    
    yy=sol.rhs
    C1=Symbol('C1')
    z= lambdify(x,sol.rhs-y0)
    c=solve(z(x0),C1)
    ysol=yy.subs({C1:c[0]})
    return ysol



In [6]:

    
ysol=ode1(f_,0,1)



In [7]:

    
print ysol









    



((x - 1)*exp(x) + 2.0)*exp(-x)



In [8]:

    
z=lambdify(x,ysol)
z(1)









    Out[8]:





0.7357588823428847



In [9]:

    
t=linspace(-4.,4.,101)



In [67]:

    
def makepic(t,f_,xin,yin):
    yy=np.array([])
    ysol=ode1(f_,xin,yin)
    z=lambdify(x,ysol)
    for t1 in t:                
        yy=append(yy,z(t1))
    plt.plot(t,yy,'r--')
        
def yexact(t,f_,xmesh,ymesh): 
    plt.figure(figsize=(5,5)) 
    plt.xlim(-4,4)
    plt.ylim(-4,4)
    for xp in xmesh:
        for yp in ymesh:
            makepic(t,f_,xp,yp)
    yy=np.array([])
    ysol=ode1(f_,0,1)
    z=lambdify(x,ysol) 
    for t1 in t:                
        yy=append(yy,z(t1))   
    plt.plot(t,yy,'b')



In [43]:

    
xmesh=linspace(-4.,4.,9); 
ymesh=linspace(-4.,4.,9);



In [68]:

    
yy=yexact(t,f_,xmesh,ymesh)

Notes

If $x-y=0$, $y'=0$ implies $ (x_0,y_0)=(x_0,x_0) $ being equilibrium.
If $ x_0 >y_0$, $y'(x_0^+)>0 \to y$ increaseing and $y\sim x (\to \infty)$ as $x\to \infty$
If $ x_0 <y_0$, $y'(x_0^+)<0 \to y$ decreasing and then increasing ( since $C_1\exp(-x)<x$ eventually and referece above graph) and $y\to \infty$ as $x\to -\infty$

Question:

Explain the results in above note from the theoretical solution of $y(x)$.



In [ ]: