from resources306 import *
%config InlineBackend.figure_format = 'retina'

def eigen(F,X):  # compute eigenvalues and eigenvectors of linearization of F at X
    x0,y0 = X
    x,y = sp.symbols('x,y')
    f,g = F((x,y))
    J = sp.Matrix([[sp.diff(f,x),sp.diff(f,y)],
                   [sp.diff(g,x),sp.diff(g,y)]]).subs({x:x0,y:y0})
    return J.eigenvects()

# Runge-Kutta
# more accurate replacement for Euler

b21 = 1/2
b31 = 0
b32 = 1/2
b41 = 0
b42 = 0
b43 = 1
b51 = 1/6
b52 = 1/3
b53 = 1/3
b54 = 1/6

def rk4_step(f,t,y,h):
    #print('rk4_step',y,type(y))
    K1 = f(t,y)  
    y1 = y + h*b21*K1
    K2 = f(t+h*b21,y1)
    y2 = y + h*(b31*K1 + b32*K2)
    K3 = f(t+h*(b31    + b32   ), y2)
    y3 = y + h*(b41*K1 + b42*K2 + b43*K3)
    K4 = f(t+h*(b41    + b42    + b43   ), y3)
    y4 = y + h*(b51*K1 + b52*K2 + b53*K3 + b54*K4)
    return y4
        
def rk4(f,t0,t1,y0,h,box=2):
    #print('rk4')
    t = t0
    y = np.array(y0)
    datay = [y]
    datat = [t0]
    #print(t,y,max(np.abs(y)))
    while t<t1 and max(np.abs(y))<=box:
        y = rk4_step(f,t,y,h)
        datay.append(y)
        t += h
        datat.append(t)
    return np.array(datat),np.array(datay)

def curve(f,t0,t1,y0,h,box=2,*args):
    datat,datay = rk4(f,t0,t1,np.array(y0),h,box=box)
    plt.plot(datay[:,0],datay[:,1],lw=3,*args)

Example 4 from Day 18

Linearization at (0,0) has a center. Is (0,0) a center of the nonlinear system?

Numerically it appears so - see figure below.

# another center in the linearization
def F(X):
    x,y = X
    return y,-x + x**2 

x0,y0 = 0,0
print(f'eigendata at {x0},{y0}:')
display(eigen(F,(x0,y0)))

x0,y0 = 1,0
print(f'eigendata at {x0},{y0}:')
display(eigen(F,(x0,y0)))

def FF(t,X): return np.array(F(X))
sin = np.sin
plt.figure(figsize=(10,10))
b = 1.5
fieldplot(F,-b,b,-b,b,alpha=0.5,color='gray')
def circle(n):
    thetas = np.linspace(0,2*np.pi,n,endpoint=False)
    return np.array([[np.cos(theta),np.sin(theta)] for theta in thetas])

for y0 in np.linspace(-1.5,0,10):
    curve(FF,0,10,[1.5,y0],.05)
for x0 in np.linspace(0,1.,10):
    curve(FF,0,10,[x0,0],.05)

for x,y in [(0,0),(1,0)]:
    p, = plt.plot(x,y,'o',markersize=10,alpha=.9)
plt.title('a center in the linearization ... and apparently also in the nonlinear system');
plt.xlabel('x'); plt.ylabel('y',rotation=0);
plt.savefig('day18_s25/phase_portrait_example_3.png')

eigendata at 0,0:

$\displaystyle \left[ \left( - i, \ 1, \ \left[ \left[\begin{matrix}i\\1\end{matrix}\right]\right]\right), \ \left( i, \ 1, \ \left[ \left[\begin{matrix}- i\\1\end{matrix}\right]\right]\right)\right]$

eigendata at 1,0:

$\displaystyle \left[ \left( -1, \ 1, \ \left[ \left[\begin{matrix}-1\\1\end{matrix}\right]\right]\right), \ \left( 1, \ 1, \ \left[ \left[\begin{matrix}1\\1\end{matrix}\right]\right]\right)\right]$

We find that $H(x,y) = \frac{1}{2}x^2 + \frac{1}{2}y^2 + \frac{1}{3}x^3$ is a conserved quantity for this system, and its level curves near (0,0) are closed curves (approximately circles):

delta = 0.025
x = np.arange(-1.5, 1.5, delta)
y = np.arange(-1.5, 1.5, delta)
X, Y = np.meshgrid(x, y)

H = X**2/2 + Y**2/2 - X**3/3

fig, ax = plt.subplots(figsize=(10,10))
ax.set_aspect(1)
CS = ax.contourf(X, Y, H, levels=np.linspace(0,2,21)**2,cmap='cool')
CS = ax.contour(X, Y, H, levels=np.linspace(0,1,21)**2,colors='k')
plt.title('level curves of $H(x,y) = x^2/2 + y^2/2 - x^3/3$');
plt.xlabel('x'); plt.ylabel('y',rotation=0);

Text(0, 0.5, 'y')

Back to example 5, the undamped pendulum

Noting that $$𝐻(𝑥,𝑦)=12𝑦^2−\cos 𝑥$$ is conserved on solution curves, we can plot level curves of H to get solution curves:

def H(X,Y):
    return Y**2/2 - np.cos(X)
x = np.linspace(-3*np.pi,4*np.pi,211)
y = np.linspace(-3,3,30)
X,Y = np.meshgrid(x,y)
plt.figure(figsize=(21,6))
plt.contour (X,Y,H(X,Y),levels=10,colors='k')
plt.contourf(X,Y,H(X,Y),levels=10,cmap='cool')
plt.axis('equal');

Rabbits and Foxes

The rabbits-and-foxes predator-prey model is another nonlinear system with an equilibrium with a center in the linearization.

Indeed numerical solutions suggest it is a center in the nonlinear system, with closed-loop solution curves:

# Rabbits and Foxes
a,b,c,d = .4,.3,.01,.003
def F(X):
    R,F = X
    return a*R - c*R*F, -b*F + d*R*F
def FF(t,X): return np.array(F(X))
sin = np.sin
plt.figure(figsize=(10,10))

fieldplot(F,0,300,0,100,alpha=0.5,color='gray')

for R0 in np.linspace(0,100,20):
    curve(FF,0,20,[R0,40],.05,box=300)

for x,y in [(0,0),(100,40)]:
    p, = plt.plot(x,y,'o',markersize=10,alpha=.9,clip_on=False)
plt.title('rabbits and foxes'); plt.xlabel('rabbit population, R'); plt.ylabel('fox population, F');

Here we also find a conserved quantity, which is (not obvious how to find this):

$H(x,y)= y^a x^b e^{-dx-cy}$

and we can plot its contours:

def H(x,y): return  y**a * x**b * np.exp(-d*x-c*y)
x = np.linspace(50,150,31)
y = np.linspace(20,60,31)
X,Y = np.meshgrid(x,y)
plt.figure(figsize=(10,10))
plt.contour (X,Y,H(X,Y),levels=np.linspace(8,9.2,20),colors='k')
plt.contourf(X,Y,H(X,Y),levels=np.linspace(8,9.2,20),cmap='cool');
plt.plot(b/d,a/c,'o',color='orange',markersize=10,alpha=1,clip_on=False)

$\displaystyle 8.64605231880377$

Other qualitative aspects of 2D nonlinear systems

Limit cycle

A closed solution curve that nearby solutions approach asymptotically.

Example: the Van der Pol oscillator

$$ x^\prime = y $$$$ y^\prime = (1-x^2)y - x$$

from IPython.display import Image
Image('day20_s25/day20_vanderpol_portrait.png')

What is the scope of behavior in 2D autonomous systems?

We have a theorem that in a 2D autonomous system, solutions can have only one of the following asymptotic behaviors:

• fly off to infinity
• approach an equilibrium
• be a cycle
• approach a limit cycle
• approach a cycle linking a set of equilibria

There are no other possibilities in 2D.

In 3D (and higher)

there are other possibilities ... including perpetually irregular (chaotic) trajectories, as in the Lorenz system: