SNB I-3 Euclidean phase portrait.ipynb

Kernel: SageMath (stable)

9 Euclidean signature phase portrait

This SageMath notebook uses the optional TIDES package which provides an arbitrary precision ode solver.

To add the TIDES package to SageMath (on a unix machine) run in a shell     $ sage -i tides as instructed at http://doc.sagemath.org/html/en/reference/misc/sage/misc/package.html

The euclidean signature ode

$$\partial_{\tau}f_{\tau} = 2-2 f_{\tau}^{2} -2 f_{\tau}v_{\tau} \qquad \partial_{\tau}v_{\tau} = -3+3 f_{\tau}^{2}+ 4f_{\tau} v_{\tau}$$

can be regarded as a dynamical system. A solution of the ode is a trajectory in the $f_\tau, v_\tau$ plane parametrized by $\tau$. We investigate the phase portrait of this dynamical system.

The ode is invariant under the time-reflection symmetry

$$\tau \rightarrow -\tau \qquad f_\tau,v_\tau \rightarrow -f_\tau,-v_\tau$$

There are two fixed points, at $(1, 0)$ and at its reflection $(-1,0)$. The linearized ode at $(1,0)$ is

$$\partial_{\tau} \begin{pmatrix} \Delta f_{\tau}+\Delta v_{\tau}\\ 3\Delta f_{\tau}+\Delta v_{\tau} \end{pmatrix} = \begin{pmatrix} 2 & 0\\ 0 & -2 \end{pmatrix} \begin{pmatrix} \Delta f_{\tau}+\Delta v_{\tau}\\ 3\Delta f_{\tau}+\Delta v_{\tau} \end{pmatrix}$$

so the fixed point has an unstable manifold tangent to $\Delta v_\tau = -3 \Delta f_\tau$ and a stable manifold tangent to $\Delta v_\tau = - \Delta f_\tau$.

Expand around $f_{\tau}=\pm\infty$

$\newcommand{\ftau}{\mathbf{f}} \def\vtau{\mathbf{v}}$ The trajectories at large $f_{\tau}$ are asymptotic to one of four lines in the $f_{\tau},v_{\tau}$ plane

$$\left(f_{\tau}(\tau) ,v_{\tau}(\tau)\right) \;\rightarrow\; \frac1\tau \left(\frac{1+\epsilon\sqrt3}{2}, \frac{-\epsilon\sqrt3}2 \right) \qquad \epsilon = \pm 1 \quad \tau \rightarrow 0^{\pm}$$

This can be seen by changing variable from $v_\tau$ to $v_\tau/f_\tau$ and then taking the limit $f_\tau \rightarrow \infty$.

Consider the trajectories with $\tau\rightarrow 0^{+}$, i.e., $\tau\in(0,\cdot)$. These are the trajectories incoming from $f_\tau=\pm\infty$. Their reflections are the trajectories outgoing towards $f_\tau=\pm\infty$.

The asymptotic trajectories are parametrized by a real number $\beta$ which ranges over some interval around $0$. The trajectory is given by an expansion around $\tau=0$

$$f_{\tau,\beta\epsilon}(\tau) = \tau^{-1}\sum_{m,n=0} (\beta \tau^{\alpha})^{m}\tau^{2n} \ftau_{m,n} \qquad v_{\tau,\beta\epsilon}(\tau) = \tau^{-1}\sum_{m,n=0} (\beta \tau^{\alpha})^{m}\tau^{2n} \vtau_{m,n}$$

with $$ \alpha = 3+\epsilon\sqrt3 \qquad $$\begin{pmatrix} \ftau_{0,0}\\ \vtau_{0,0} \end{pmatrix}$$

$$\begin{pmatrix} \left(1+\epsilon\sqrt3\right)/2 \\ -\epsilon\sqrt3/2 \end{pmatrix}$$

\qquad $$\begin{pmatrix} \ftau_{1,0}\\ \vtau_{1,0} \end{pmatrix}$$

$$\begin{pmatrix} 1\\ -1-\epsilon \sqrt3 \end{pmatrix}$$

$$

The coefficients for $(m,n)\ne (0,0),(1,0)$ are determined by recursion relations

$$A_{m,n} \begin{pmatrix} \ftau_{m,n}\\ \vtau_{m,n} \end{pmatrix} = \delta_{m,0}\delta_{n,1} \begin{pmatrix} 2 \\ -3 \end{pmatrix} + \begin{pmatrix} -2 & -2\\ 3 &4 \end{pmatrix} \sum_{\substack{m'+m''=m\\ n'+n''=n\\ \ne (0,0),(m,n)}}\ftau_{m',n'} \begin{pmatrix} \ftau_{m'',n''}\\ \vtau_{m'',n''} \end{pmatrix}$$$$\begin{aligned} A_{m,n} &=\begin{pmatrix} \alpha (m+1)+2n-2 &\alpha-2\\ -\alpha &\alpha (m-2)+2n+3 \end{pmatrix} $1ex] \det A_{m,n} &= (\alpha m+2n+1) \left(\alpha (m-1) + 2n\right) \end{aligned}$$

The phase portrait is obtained by combining the expansion at large $f_\tau$ with numerical solutions of the ode.

Calculate the large $f$ expansion numerically

First, enter three precision parameters.

decimal_precision = 120        # number of decimal digits in floating point numbers
m_order =  60                   # m = 0,1,..., m_order-1
n_order =  60                   # n = 0,1,..., n_order-1
#
binary_precision = ceil(decimal_precision*log(10,2))
R = RealField(binary_precision); RealNumber = R
binary_precision

399

Construct the expansion as function

$$\mathbf{fvoftau}(\tau,\epsilon,\beta) = \left[\mathbf{foftau}, \: \mathbf{voftau}\right]\,(\tau,\epsilon,\beta)$$

zero2vector = vector(R,2)
Rxy.<x,y> = R[]        # a polynomial ring in x,y
fvtau,eps_sqrt3,alpha,fpols=var('fvtau eps_sqrt3 alpha fpols')
MatrixRHS = matrix(R,[[-2,-2],[3,4]])
A = lambda m,n: matrix(R,[[(m+1)*alpha+2*n-2,alpha-2],
                     [-alpha,(m-2)*alpha+2*n+3]])
def fvSumvect(m,n):
    fvSumvec = zero2vector
    for mp,np in mrange([m+1,n+1]):
        fvSumvec += fvtau[mp][np][0]* fvtau[m-mp][n-np]
    return(fvSumvec)
def make_fvtau(epsilon):
    global alpha,fvtau
    eps_sqrt3 = epsilon*sqrt(3.0)
    alpha = 3.0+eps_sqrt3
    fvtau=[[zero2vector for m in range(m_order)] for n in range(n_order)]
    fvtau[0][0] = vector(R,[(1+eps_sqrt3)/2,-eps_sqrt3/2])
    fvtau[1][0] = vector(R,[1,-eps_sqrt3 -1])
    for m in range(2,m_order):
        fvtau[m][0] = (A(m,0)^(-1))*MatrixRHS*fvSumvect(m,0)
    fvtau[0][1]=(A(0,1)^(-1))*(MatrixRHS*fvSumvect(0,1)+ vector(R,[2,-3]))
    for n in range(2,n_order):
        fvtau[0][n] = (A(0,n)^(-1))*MatrixRHS*fvSumvect(0,n)
    for m in range(2,m_order):
        for n in range(2,n_order):
            fvtau[m][n] = (A(m,n)^(-1))*MatrixRHS*fvSumvect(m,n)
def fvoftau(tau,epsilon,beta):
    alpha = 3.0+epsilon*sqrt(3.0)
    fpol,vpol = fvpols[epsilon]
    betataualpha = beta*(tau^alpha)
    return ([fpol(tau,betataualpha)/tau,vpol(tau,betataualpha)/tau])
# construct polynomials fpol(y,x) = sum_{m,n} ftau_{m,n} y^m x^n 
#                       vpol(y,x) = sum_{m,n} vtau_{m,n} y^m x^n
fvpols = {}
for epsilon in {+1,-1}:
    make_fvtau(epsilon)
    fpol= sum( [(y^m) * (x^(2*n))*fvtau[m][n][0]
                for m,n in mrange([m_order,n_order])])
    vpol= sum( [(y^m) * (x^(2*n))*fvtau[m][n][1]
                for m,n in mrange([m_order,n_order])])
    fvpols[epsilon] = (fpol,vpol)

Calculate trajectories numerically

Use the TIDES numerical ode solver to calculate trajectories with initial conditions given by the large $f$ expansion.

sage.calculus.desolvers.desolve_tides_mpfr(f, ics, initial, final, delta, tolrel=1e-16, tolabs=1e-16, digits=50) Solve numerically a system of first order differential equations using the taylor series integrator in arbitrary precision implemented in tides.

INPUT:

    f – symbolic function. Its first argument will be the independent variable. Its output should be derivatives of the dependent variables.
    ics – a list or tuple with the initial conditions.
    initial – the starting value for the independent variable.
    final – the final value for the independent value.
    delta – the size of the steps in the output.
    tolrel – the relative tolerance for the method.
    tolabs – the absolute tolerance for the method.
    digits – the digits of precision used in the computation.

Set: 1. the step size in $\tau$ for the trajectory returned by the ode solver as a list of points $(f,v)$. 2. the tolerances for the ode solver.

solver_step_size = 0.01   # step_size in \tau for the trajectories returned by the ode solver
solver_tolrel = 1e-60      # the relative tolerance for the solver
solver_tolabs = 1e-60      # the absolute tolerance for the solver

Define a function $\mathbf{calc\_traj}\,(epsilon,\,beta,\,tau\_i,\,tau\_f)$ which returns the trajectory as a list of points in the form $\left [\tau,f_{\tau,\beta\epsilon}(\tau),v_{\tau,\beta\epsilon}(\tau)\right]$ with $\tau$ ranging from $tau\_i$ to $tau\_f$ in steps of size $solver\_step\_size$.

ode_function(tau,f,v) = [2-2*f^2-2*v*f,-3+3*f^2+4*v*f]
def calc_traj(epsilon,beta,tau_i,tau_f):
    [f_i,v_i]=fvoftau(tau_i,epsilon,beta)
    traj_list=desolve_tides_mpfr(ode_function, [f_i,v_i], tau_i, tau_f, solver_step_size,solver_tolrel, solver_tolabs, decimal_precision)
    return(traj_list)

fmax = 3.0
vmax = 6.0
def choptraj(traj):
    traj_list = []
    for tau,f,v in traj:
        if abs(f) <= fmax and abs(v) <= vmax:
            traj_list.append([f,v])
    return(traj_list)
my_arrow = lambda p1,p2: arrow(p1,p2,arrowsize=2.5,width=0.1)
def plot_traj(traj,arrow1ratio=1/40,arrow2ratio=1/3,labelratio=1/100,trajlabel=""):
    traj_list=choptraj(traj)
    LP=list_plot(traj_list,size=1)
    npoints=len(traj_list)
    arrow1index=int(npoints*arrow1ratio)
    arrow2index=int(npoints*arrow2ratio)
    labelindex=int(npoints*labelratio)
    LP+= my_arrow(traj_list[arrow1index],traj_list[arrow1index+1])
    LP+= my_arrow(traj_list[arrow2index],traj_list[arrow2index+1])
    offset=fmax/15
    temp=traj_list[labelindex]
    labelpos=(temp[0]+offset,temp[1]+offset)
#    LP+=text(trajlabel,labelpos,fontsize='smaller')
    traj_list_refl = [[-f,-v] for f,v in traj_list]
    LP+=list_plot(traj_list_refl,size=1)
    LP+= my_arrow(traj_list_refl[arrow1index+1],traj_list_refl[arrow1index])
    LP+= my_arrow(traj_list_refl[arrow2index+1],traj_list_refl[arrow2index])
    temp=traj_list_refl[labelindex]
    labelpos=(temp[0]+offset,temp[1]+offset)
    #LP+=text(trajlabel,labelpos)
    return(LP)

def traj_plotted(epsilon,beta,tau_i,tau_f,arrow1ratio=1/40,arrow2ratio=1/3,labelratio=1/100,trajlabel=""):
    traj=calc_traj(epsilon,beta,tau_i,tau_f)
    return(plot_traj(traj,arrow1ratio,arrow2ratio,labelratio,trajlabel))

portrait = traj_plotted(epsilon=-1,beta=0.0,tau_i=0.13,tau_f=2.8,trajlabel=r"$-,0.0$")
portrait+= traj_plotted(epsilon=1,beta=0.0,tau_i=0.45,tau_f=2.8,trajlabel=r"$+,0.0$")
portrait+= traj_plotted(epsilon=-1,beta=0.5,tau_i=0.11,tau_f=5.5,trajlabel=r"$-,0.5$")
portrait+= traj_plotted(epsilon=-1,beta=0.05,tau_i=0.11,tau_f=5.5,trajlabel=r"$-,0.41$",
                    arrow1ratio=1/10,arrow2ratio=.4  )
portrait+= traj_plotted(epsilon=-1,beta=2.0,tau_i=0.17,tau_f=2.5,trajlabel=r"$-,2.0$")
portrait+= traj_plotted(epsilon=-1,beta=-0.4,tau_i=0.1,tau_f=1.5,trajlabel=r"$-,-.4$",
                       arrow1ratio=1/20,arrow2ratio=19/20)
portrait+= traj_plotted(epsilon=-1,beta=-1.109257,tau_i=0.01,tau_f=2.7,trajlabel=r"$-,-1.1$",
                    arrow1ratio=1/60,arrow2ratio=1/8  )
#good_portrait=portrait
#show(good_portrait)
portrait+=traj_plotted(epsilon=-1,beta=-1.3,tau_i=0.08,tau_f=1.05,trajlabel=r"$-,1.3$",
                    arrow1ratio=1/10,arrow2ratio=.8  )
#show(good_portrait)
portrait+= traj_plotted(epsilon=+1,beta=-0.1,tau_i=0.35,tau_f=4.0,trajlabel=r"$-,.03$",
                    arrow1ratio=1/10,arrow2ratio=.8  )
portrait+= traj_plotted(epsilon=+1,beta=-1.5,tau_i=0.35,tau_f=4.0,trajlabel=r"$-,.03$",
                    arrow1ratio=1/10,arrow2ratio=.8  )
portrait+= traj_plotted(epsilon=+1,beta=-10,tau_i=0.35,tau_f=4.0,trajlabel=r"$-,.03$",
                    arrow1ratio=1/10,arrow2ratio=.8  )
portrait+= traj_plotted(epsilon=+1,beta=-100,tau_i=0.3,tau_f=0.6,trajlabel=r"$-,.03$",
                    arrow1ratio=1/10,arrow2ratio=.8  )
portrait.fontsize(4)
#show(portrait)
plot(portrait).save('euclidean_phase_portrait.pdf',dpi=200,aspect_ratio=0.5,axes_labels=[r"$f_T$",r"$v_T$"],
    axes_labels_size=2)
show(portrait,dpi=200,aspect_ratio=0.5,axes_labels=[r"$f_T$",r"$v_T$"],
    axes_labels_size=2)