Kdtyjb

I am trying to incorporate the solution to the nonlinear differential equations

x' = x + y + x^2 + y^2

y' = x - y - x^2 + y^2

The solution to these nonlinear differential equations were found with a saddle point at (0,0) and an equilibrium point at (-1, -1). The equation was solved using Matlab and produced this result:

I have been trying to plot some of the lines and simplified the plot to produce these points with

x0=[-2.5,-1.5,-1, -0.5]; and

y0=-2.5.

This is the result:

I am trying to elegantly plot this and came across pst-ode. I attempted to get a plot to match but so far have failed miserably! I followed the code given here: Differential Equation direction plot with pgfplots but still no luck. Can you help me get the correct plot to match the original plot showing the lines.

I'm uncertain how to incorporate the second differential equation to come up with the correct plot. Thanks for your help!

Here is my attempt this far:

documentclass[border=10pt]{standalone}

usepackage{pst-plot,pst-ode}

begin{document}



psset{unit=3}

begin{pspicture}(-1.2,-1.2)(1.1,1.1)

psaxes[ticksize=0 4pt,axesstyle=frame,tickstyle=inner,subticks=20,

        Ox=-2.5,Oy=-2.5](-1,-1)(1,1)

psset{arrows=->,algebraic}

psVectorfield[linecolor=black!60](-0.9,-0.9)(0.9,0.9){ x + y + x^2 + y^2 }

%y0_a=-2.5

pstODEsolve[algebraicOutputFormat]{y0_a}{t | x[0]}{-1}{1}{100}{-2.5}{t + x[0] + t^2 + x[0]^2}

%y0_b=-1.5

pstODEsolve[algebraicOutputFormat]{y0_b}{t | x[0]}{-1}{1}{100}{-1.5}{t + x[0] + t^2 + x[0]^2}

%y0_c=-1

pstODEsolve[algebraicOutputFormat]{y0_c}{t | x[0]}{-1}{1}{100}{-1}{t + x[0] + t^2 + x[0]^2}

%y0_c=-0.5

pstODEsolve[algebraicOutputFormat]{y0_d}{t | x[0]}{-1}{1}{100}{-0.5}{t + x[0] + t^2 + x[0]^2}





psset{arrows=-,linewidth=1pt}%

listplot[linecolor=red  ]{y0_a}

listplot[linecolor=green]{y0_b}

listplot[linecolor=blue ]{y0_c}

listplot[linecolor=blue ]{y0_d}

end{pspicture}



end{document} 

The right-hand side of the given system of ODEs

dx/dt = x + y + x^2 + y^2

dy/dt = x − y − x^2 + y^2

translates into algebraic notation, used by pstODEsolve as its last argument, as

defodeRHS{

  x[0] + x[1] + x[0]^2 + x[1]^2

|

  x[0] − x[1] − x[0]^2 + x[1]^2

}

Here, in the given case, the RHS does not depend on t, the independent integration parameter.

We solve the set of ODEs many times with varying initial conditions, that is, starting points x0 and y0 in the x-y plane. Those starting points are aligned on a 0.1 by 0.1 raster which is shifted by 0.05 in x and y to avoid the saddle point as starting point. (This would not do any harm, but produce a coloured dot at (0,0).) The starting points are specified in algebraic notation as the last-but-one argument of pstODEsolve.

The example needs pst-ode v0.8 [2017/06/12] in order to compile.

The code is quite demanding in terms of TeX memory due to the nested loops. The TeX memory ↗can be enlarged for the latex engine, but it may be easier to just use lualatex instead:

lualatex phasePlane

dvips phasePlane

ps2pdf phasePlane.ps

The code:

ifdefinedoutputmodeoutputmode=0fi %enforce DVI output



documentclass{article}

usepackage[tightpage,active]{preview}

usepackage{pst-plot,pst-ode,multido}



%right hand side in algebraic notation

defodeRHS{

  x[0] + x[1] + x[0]^2 + x[1]^2

|

  x[0] - x[1] - x[0]^2 + x[1]^2

}



deftEnd{10} % end of integration interval; in- or de-crease to have longer/shorter lines

defoutputSteps{500} %increase for smoother lines (and larger PDF file)



definecolor[ps]{random}{rgb}{Rand Rand Rand} %define random colour



begin{document}

begin{preview}



% Solve ODE for initial conditions (starting points) on a 0.1 by 0.1 raster

multido{ii=0+1,rx=-2.55+0.1}{42}{

  multido{ij=0+1,ry=-2.55+0.1}{42}{

    pstODEsolve[algebraicAll]{line_ii_ij}{x[0] | x[1]}{0}{tEnd}{outputSteps}{rx | ry}{odeRHS}

  } 

}



psset{unit=3}

begin{pspicture}(-2.95,-2.8)(1.6,1.6)

  psaxes[ticksize=0 4pt,axesstyle=frame,tickstyle=inner,

    dx=0.5,dy=0.5,Dx=0.5,Dy=0.5,Ox=-2.5,Oy=-2.5,subticks=5

  ](-2.5,-2.5)(1.5,1.5)



  rput(-0.5,-2.75){$x$}

  rput(-2.9,-0.5){$y$}



  psclip{psframe[linestyle=none](-2.5,-2.5)(1.5,1.5)}%



  %plot the lines previously computed

  multido{ii=0+1}{42}{

    multido{ij=0+1}{42}{

      listplot[linecolor=random]{line_ii_ij}

    }  

  }



  endpsclip

end{pspicture}

end{preview}



end{document}