Code to make a phase portrait of a 2x2 linear constant-coefficient system¶

In [23]:
%config InLineBackend.figure_format = 'retina'
from resources306 import *

from random import random
def randomhexcolor( pallor=0.0 ):
    s ='#'
    for i in [1,2,3]:
        if pallor >= 0:
            r = int( 256*pallor + 256*(1-pallor)*random() )  # pallor = whiteness
        else:
            r = int(              256*(1+pallor)*random() )  # pallor = whiteness
        s += hex(r)[2:].zfill(2)
    return s
rc = randomhexcolor

def vectors(A,color=None):
    #display(sp.Matrix(A))
    plt.figure(figsize=(10,10))
    if color is None: color = rc(pallor=0.7)
    plt.fill([-1,1,1,-1],[-1,-1,1,1],color=color,)
    color = rc(pallor=0.1)
    plt.plot(0,0,'o',color='k',alpha=0.8)
    fieldplotlinear(A,-1,1,-1,1,alpha=0.8,color='k',nx=16,boostarrows=1.25 )
    plt.xticks([])
    plt.yticks([])
    plt.title(f'{A[0][0]:2} {round(A[0][1],2):4}\n{A[1][0]:2} {round(A[1][1],2):4}',fontsize=30)
In [26]:
# version with numerical solutions instead of manually constructed formulas  3/24/2025

import os

def portrait(A,color=None,name=None):
    
    def F(X): 
        x,y = X
        return np.array([A[0][0]*x + A[0][1]*y, A[1][0]*x + A[1][1]*y])

    vectors(A,color=color)    

    eigvals,eigvecs = np.linalg.eig(A)
    #print(eigvals)
    #print(eigvecs)
    if np.allclose(np.imag(eigvals),0):  # real eigenvalues: draw eigenvectors
        for v in eigvecs.T:
            v *= np.sqrt(2)  # so goes all the way to edge of frame
            plt.plot([0, v[0]],[0, v[1]],'k',lw=5,alpha=0.5)
            plt.plot([0,-v[0]],[0,-v[1]],'k',lw=5,alpha=0.5)
        #assert 0
    for theta in np.linspace(0,2*np.pi,20,endpoint=False):
        x0,y0 = np.cos(theta),np.sin(theta)
        phaseportrait(F,[(x0,y0,-3,3)],color='m',alpha=0.5,lw=5)
    if name is None:
        filename = f'2x2_linear_{A[0][0]}_{A[0][1]}_{A[1][0]}_{A[1][1]}.png'
    else:
        filename = name+'.png' #'.svg'
    plt.savefig(filename,bbox_inches='tight')
In [27]:
portrait(
    [[1,2],
     [-3,1]] 
)
In [ ]: