10

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:

enter image description here

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:

enter image description here

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} 
3
  • You mean you tried quiver from pgfplots? Can you show us a/the failing code?
    – Symbol 1
    Jun 8, 2017 at 2:48
  • I did not try quiver. I made my first attempt here with pst-ode. Am I going down the wrong path?
    – Joe
    Jun 8, 2017 at 2:53
  • Probably the right path. PS can do more computations than PGF. But I am not familiar with the magical PStricks family.
    – Symbol 1
    Jun 8, 2017 at 3:00

1 Answer 1

9

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

\def\odeRHS{
  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.

enter image description here 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 dvilualatex instead:

dvilualatex phasePlane
dvips phasePlane
ps2pdf phasePlane.ps

Direct PDF output with lualatex is possible nowadays, thanks to Marcel Krüger's PostScript interpreter written in Lua. (It is just a bit slower than ps2pdf):

lualatex phasePlane

The code:

\documentclass{standalone}

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

%right hand side in algebraic notation
\def\odeRHS{
  x[0] + x[1] + x[0]^2 + x[1]^2
|
  x[0] - x[1] - x[0]^2 + x[1]^2
}

\def\tEnd{10} % end of integration interval; in- or de-crease to have longer/shorter lines
\def\outputSteps{500} %increase for smoother lines (and larger PDF file)

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

\begin{document}

\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)}%

  % 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}
    } 
  }
  %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{document}
17
  • 1
    I may need to revise the pst-ode code in order catch the numerical underflow error, preventing the code from chrashing but just stopping the integration instead.
    – AlexG
    Jun 8, 2017 at 18:04
  • 1
    Ok, pst-ode has been updated (v0.8, [2017/06/12]). This should answer most of your questions.
    – AlexG
    Jun 13, 2017 at 14:26
  • 1
    Sincerest Thanks for providing this elegant solution!
    – Joe
    Jun 14, 2017 at 5:05
  • 1
    I was looking for a way to randomly set the line colours. Now we get a result which is close to the Matlab output. I updated pst-ode once more (v 0.9), but the code should happily compile with v 0.8. Occasionally, the alogrithm entered an infinite loop. Thanks for the interesting topic which helped me to improve the pkg.
    – AlexG
    Jun 15, 2017 at 8:28
  • 1
    also interesting: \definecolor[ps]{random}{cmyk}{0.2 0.2 0.2 Rand}
    – user2478
    Jun 16, 2017 at 20:04

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .