Differential Equation direction plot with pgfplots

In earlier thread Jake provided some code whom successfully draws the following differential equations in the range [0; 1]

dy/dx=2*x
dy/dx=x*sqrt(x)


I have not been able to determine how to change the range to [-1; 1] for both dimensions. My attempt

xmin=-1.1, xmax=1.1, % Axis limits
ymin=-1.1, ymax=1.1,
domain=-1:1, y domain=-1:1,


Also, I had some trouble with the notation for a differential equation that consists of y, example (dy/dx=x^2+y^2-1).

My attempt:

declare function={f(\x) = \x^2 + f(\x)^2 - 1;}


Using PGFPlotstable, you can use a naive numerical integration scheme to find the function directly within LaTeX:

\documentclass{article}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=1.8}

\usepackage{amsmath}

\pgfplotstableset{
create on use/x/.style={
create col/expr={
\pgfplotstablerow/201*2-1
}
},
create on use/y/.style={
create col/expr accum={
\pgfmathaccuma+(2/201)*(abs(\pgfmathaccuma^2)+abs(\thisrow{x}^2)-1)
}{0.6}
}
}

\begin{document}
\begin{tikzpicture}
\begin{axis}[
view={0}{90},
domain=-1:1,
y domain=-1:1,
xmax=1, ymax=1,
samples=21
]
\addplot3 [gray, quiver={u={1}, v={x^2+y^2-1}, scale arrows=0.075, every arrow/.append style={-latex}}] (x,y,0);
\end{axis}
\end{tikzpicture}

\end{document}

• How to you change the x and y values for the viewing window? – MaoYiyi Oct 17 '13 at 13:00

(Sets of) ordinary differential equations (ODE) can be numerically solved using \pstODEsolve from the pst-ode Plain-TeX / LaTeX package.

Integration of the ODEs is done using the Runge-Kutta-Fehlberg method of 4th order with adaptive step size control (RKF45) during the ps2pdf conversion step.

For plotting packages other than PSTricks, such as pgfplots to be used here, LaTeX must be run twice. The first run, must go the dvips+ps2pdf route, in order to write the solution table into a text file, which is then read and transformed into a chart by your favourite plotting package and TeX engine.

The first compilation is done like this (option -dNOSAFER is necessary to allow Ghostscript to write files):

latex myfile
dvips myfile
ps2pdf -dNOSAFER myfile.ps


The second run can involve your preferred toolchain pdflatex, xelatex, latex/dvips/ps2pdf, whatever. It plots the results using pgfplots.

Alternatively, a single call to your favourite LaTeX engine with option --shell-escape is sufficient if the second code box below is used. (Under the hood, the chain of commands tex, dvips, ps2pdf -dNOSAFER is executed before the document is actually typeset.) Solution provided by Herbert.

\documentclass{article}
\usepackage{amsmath}
\usepackage{pgfplots}
\pgfplotsset{compat=1.8}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{ifpdf}
\ifpdf
\IfFileExists{y0=-0.5.dat}{}{
\GenericError{}{\MessageBreak%
===================================\MessageBreak
First, you have to run\MessageBreak\MessageBreak
latex \jobname\MessageBreak
dvips \jobname\MessageBreak
ps2pdf -dNOSAFER \jobname.ps\MessageBreak\MessageBreak
to solve the differential equation.\MessageBreak
===================================}{}{}
}
\else
\usepackage{pst-ode}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%solve dy/dx=x^2 + y^2 - 1 numerically for different initial values of y in the
%interval x=[-1.1,1.1]; write resulting curves as tables with 100 output points
%into text files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%y0=-0.5 --> y0=-0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{%
y0=-0.5%  %name of the output file y0=-0.5.dat'
}{
t | x[0]  %table format in y0=-0.5.dat': x y
}{
-1.1      %integration domain x_min=-1.1
}{
1.1       %integration domain x_max=1.1
}{
100       %number of output points
}{
-0.5      %initial value y0(x_min)
}{
t^2+x[0]^2-1  % right hand side of ODE, note the special notation:  x --> t, y --> x[0]
}

%y0=0.0 --> y0=0.0.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.0}{t | x[0]}{-1.1}{1.1}{100}{0.0}{t^2+x[0]^2-1}

%y0=0.5 --> y0=0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.5}{t | x[0]}{-1.1}{1.1}{100}{0.5}{t^2+x[0]^2-1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\fi

\begin{document}
\IfFileExists{y0=-0.5.dat}{}{dummy text\end{document}}

\begin{tikzpicture}
\begin{axis}[
axis equal image, % Unit vectors for both axes have the same length
xmin=-1.1, xmax=1.1, % Axis limits
ymin=-1.1, ymax=1.1,
xtick={-1,-0.5,0,0.5,1}, ytick={-1,-0.5,0,0.5,1}, % Tick marks
no markers,
title={$\dfrac{\mathrm{d}y}{\mathrm{d}x}=x^2+y^2-1$},
view={0}{90},samples=21,domain=-1.1:1.1, y domain=-1.1:1.1, %for direction field
]
\addplot3 [gray, quiver={u={1}, v={x^2+y^2-1}, scale arrows=0.075, every arrow/.append style={-latex}}] (x,y,0);
\end{axis}
\end{tikzpicture}

\end{document}


Code for typesetting at one sweep using option --shell-escape with any LaTeX engine:

\documentclass{article}
\usepackage{amsmath}
\usepackage{pgfplots}
\pgfplotsset{compat=1.8}

\usepackage{filecontents}
\begin{filecontents}{xyz.tex}
\input pst-ode
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%solve dy/dx=x^2 + y^2 - 1 numerically for different initial values of y in the
%interval x=[-1.1,1.1]; write resulting curves as tables with 100 output points
%into text files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%y0=-0.5 --> y0=-0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=-0.5}{t | x[0]}{-1.1}{1.1}{100}{-0.5}{t^2+x[0]^2-1}

%y0=0.0 --> y0=0.0.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.0}{t | x[0]}{-1.1}{1.1}{100}{0.0}{t^2+x[0]^2-1}

%y0=0.5 --> y0=0.5.dat
\pstODEsolve[algebraicOutputFormat,algebraic,saveData]{y0=0.5}{t | x[0]}{-1.1}{1.1}{100}{0.5}{t^2+x[0]^2-1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bye
\end{filecontents}

\immediate\write18{tex xyz}
\immediate\write18{dvips xyz}
\immediate\write18{ps2pdf -dNOSAFER xyz.ps}

\begin{document}

\begin{tikzpicture}
\begin{axis}[
axis equal image, % Unit vectors for both axes have the same length
xmin=-1.1, xmax=1.1, % Axis limits
ymin=-1.1, ymax=1.1,
xtick={-1,-0.5,0,0.5,1}, ytick={-1,-0.5,0,0.5,1}, % Tick marks
no markers,
title={$\dfrac{\mathrm{d}y}{\mathrm{d}x}=x^2+y^2-1$},
view={0}{90},samples=21,domain=-1.1:1.1, y domain=-1.1:1.1, %for direction field
]
\addplot3 [gray, quiver={u={1}, v={x^2+y^2-1}, scale arrows=0.075, every arrow/.append style={-latex}}] (x,y,0);
\end{axis}
\end{tikzpicture}

\end{document}

• See my edit for running everything from within one file. Then you do not need the \if... structures – user2478 Oct 17 '13 at 9:17
• @Herbert: Thanks a lot. I assimilated your suggestion in the second code block. – AlexG Oct 17 '13 at 10:41
• I don't get the three lines, I get the arrows to show, but not the lines. I get an error about not being able to read the tables. – MaoYiyi Oct 17 '13 at 13:03
• @MaoYiyi: (1) Copy contents of 2nd code block into a file, say myfile.tex. (2) Run pdflatex --shell-escape myfile. (3) Open myfile.pdf in a PDF reader of your choice. Option --shell-escape is crucial here, as it allows for automatically executing the tex,dvips,ps2pdf chain for producing the data tables in the files y0=*.dat. – AlexG Oct 17 '13 at 13:14

MWE with Asymptote

% odeslope.tex:
%
\documentclass{article}
\usepackage[inline]{asymptote}
\usepackage{lmodern}
\begin{document}
\begin{figure}
\centering
\begin{asy}
import graph;
import slopefield;
import fontsize;
defaultpen(fontsize(9pt));
size(200);
real dy(real x,real y) {return x^2+y^2-1;}
real xmin=-1, xmax=1;
real ymin=-0.2, ymax=1;

pair C=(0.5,0.4);
draw(curve(C,dy,(xmin,ymin),(xmax,ymax)),deepblue+1bp);

label("$C$",C,NE,UnFill);
dot(C,UnFill);

xaxis(YEquals(ymin),xmin,xmax,LeftTicks());
xaxis(YEquals(ymax),xmin,xmax);
yaxis(XEquals(xmin),ymin,ymax,RightTicks());
yaxis(XEquals(xmax),ymin,ymax);
\end{asy}
\caption{$\frac{\mathrm{d}y}{\mathrm{d}x}=x^2+y^2-1$}
\end{figure}
\end{document}
%
% Process:
%
% pdflatex odeslope.tex
% asy odeslope-*.asy
% pdflatex odeslope.tex

• @AlexG: The documentation is humble about it, but the code in slopefield.asy looks more like the Runge–Kutta method. – g.kov Oct 16 '13 at 14:59
• Yes and it seems to apply some automatic step size adaption. Thus, it is quite accurate. – AlexG Oct 16 '13 at 15:41
• @Marc van Dongen: Please, not again. It was once concluded after long and painful discussion that TikZ/PSTriks/Metapost/Asymptote are somewhat equivalent and the community can only benefit from answers that provide alternative approaches to get the desired high quality graphics with TeX and friends. – g.kov Oct 17 '13 at 10:40
• @MarcvanDongen: There was a bit of discussion about the relevance of Asymptote in the comments to tex.stackexchange.com/a/117436/2552, and the votes and lack of opposition on meta.tex.stackexchange.com/q/3518/2552 also seem to be quite a clear indicator that Asymptote can be considered equivalent to TikZ/PSTricks for the purposes of this site. I think the community has kind of come to the conclusion that questions about generating graphics can be answered with different tools (even if the question asks about a particular one) without having to list the pros and cons of the ... – Jake Oct 17 '13 at 10:55
• @AlexG: I know, but as I said, many questions get answers using different tools, even if the question asks for a particular one. The community has more or less agreed on this being a welcome thing (see meta.tex.stackexchange.com/questions/3408/…): Even if the answers might not be relevant to the asker because they're bound to a particular tool, other people will definitely benefit from seeing a range of different approaches. – Jake Oct 17 '13 at 11:25

A PSTricks solution, it can be run with xelatex but that takes a lot of time in difference to latex->dvips->ps2pdf

\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=-1,Oy=-1](-1,-1)(1,1)
\psset{arrows=->,algebraic}
\psVectorfield[linecolor=black!60](-0.9,-0.9)(0.9,0.9){ x^2+y^2-1 }
%y0_a=-0.5
\pstODEsolve[algebraicOutputFormat]{y0_a}{t | x[0]}{-1}{1}{100}{-0.5}{t^2+x[0]^2-1}
%y0_b=0.0
\pstODEsolve[algebraicOutputFormat]{y0_b}{t | x[0]}{-1}{1}{100}{0.0}{t^2+x[0]^2-1}
%y0_c=0.5
\pstODEsolve[algebraicOutputFormat]{y0_c}{t | x[0]}{-1}{1}{100}{0.5}{t^2+x[0]^2-1}

\psset{arrows=-,linewidth=1pt}%
\listplot[linecolor=red  ]{y0_a}
\listplot[linecolor=green]{y0_b}
\listplot[linecolor=blue ]{y0_c}
\end{pspicture}

\end{document}