# Stream plot in TikZ with singularities

I have the vector field F(x,y)=[F1(x,y),F2(x,y)] defined by

F1(x,y)= -(2*y/x^2)*(2*ln(cosh(x)) + (y-2*x)*tanh(x) - x*y*(1-(tanh(x))^2))
F2(x,y)= -(4/x)*(x*y-ln(cosh(x)) - y*tanh(x))


In Mathematica, I'm able to easily plot the corresponding stream plot via the instruction

StreamPlot[{F1[a, c], F2[a, c]}, {a, -10, 10}, {c, -10, 10}, AspectRatio -> 1/GoldenRatio]


which produces

I've been trying for a while to produce something similar with Tikz. Indeed, I tried with \addplot3 and quiver with normalization, but the singularities broke the code (e.g., in the [-10,10]x[-10,10] domain). Besides, I tried to split everything into non-singular subdomains, but although compiling, the appearance is horrible due to the "straight arrows". Find below a corresponding piece of code that doesn't break and its output:

\documentclass{standalone}

\usepackage{amsmath, tikz, pgfplots}
\usepgfplotslibrary{colormaps}

\begin{document}
\begin{tikzpicture}
\begin{axis}[
zmin = 0, zmax = 1,
axis equal image,
xtick distance = 1,
ytick distance = 1,
view = {0}{90},
scale = 1.25,
title = {$F$},
height=7cm,
xlabel = {$x$},
ylabel = {$y$},
colormap/viridis,
colorbar,
colorbar style = {
ylabel = {Vector Length}
}
]
point meta = {sqrt(( (2*y / x^2)*( 2*ln(cosh(x)) + (y-2*x)*tanh(x) - x*y*(1-(tanh(x))^2) ) )^2 + ( (4 / x)*( x*y - ln(cosh(x)) - y*tanh(x) )^2)},
quiver = {
u = {-(2*y / x^2)*( 2*ln(cosh(x)) + (y-2*x)*tanh(x) - x*y*(1-(tanh(x))^2) ) / sqrt(( (2*y / x^2)*( 2*ln(cosh(x)) + (y-2*x)*tanh(x) - x*y*(1-(tanh(x))^2) ) )^2 + ( (4 / x)*( x*y - ln(cosh(x)) - y*tanh(x) )^2)},
v = {-(4 / x)*( x*y - ln(cosh(x)) - y*tanh(x) ) / sqrt(( (2*y / x^2)*( 2*ln(cosh(x)) + (y-2*x)*tanh(x) - x*y*(1-(tanh(x))^2) ) )^2 + ( (4 / x)*( x*y - ln(cosh(x)) - y*tanh(x) )^2)},
scale arrows = 0.15,
},
quiver/colored = {mapped color},
-stealth,
domain = 0:5,
domain y = 0:5,
] {0};

\end{axis}
\end{tikzpicture}
\end{document}


Note that the above code normalizes the que vector field in F1=u, F2=v by normalizing the length of the arrow (which thereafter is indicated in the color-scale legend via point meta).

I want to obtain in Tikz something like in the first case (arrowed curves, instead of many straight arrows) that is able to deal with the "singularities".

• Hi, sometimes, TikZ isn't the best tool to do it.
– anis
May 12, 2023 at 18:32
• Maybe this can help: tex.stackexchange.com/q/339243/47927. I fear PGFplots can't do this out of the box. May 12, 2023 at 18:32

I tried to obtain a "stream plot" using basic TikZ. Each trajectory (or orbit) is constructed starting from a point and going forward and respectively backward in time. Considering your specific problem, my code does not cover the case of closed orbits. A different branch must be added inside the command \streamPlot.

The above drawing uses to sets of orbits constructed with stream plot by name and stream plot TB (top bottom). Note that these are pic objects so they should be handled with care in case scaling is used globally. They are based on \streamPlot.

I needed both since I consider that one needs to go closer to the singular points of the vector field to have a good image. For example, the image constructed with Mathematica doesn't really show what is happening near the origin.

Some remarks about the vector field that is considered: some rewriting is in order to arrive at a form that might be used with TikZ (a low computation accuracy, but enough for drawing). Below is one way of looking at the two components. I added a Taylor polynomial in x needed to understand the behavior for small values of x.

Remarks

1. Looking at the stream plot in the first drawing, it can be noticed that something is happening where trajectories "meet". To understand the issue, we can draw the vector field together with the curves U(x, y)=0 (in blue) and V(x, y)=0 (in red). See also the last image below build with Python (this is why it is at the end of my answer).

2. The singularities of the vector field do not follow the blue curve (let's say) completely. But we cannot see this in the stream plot (Fig. 1). This is why I added the solution of the corresponding differential equation in red for a certain initial condition. It uses the function ode-trajectory.

3. The algorithm behind all the computations is the classical Runge-Kutta developed in the style runge-kutta.

The code The code is rather long mainly because the vector field components are defined three times each (needed for different contexts), and because the graphical elements are of three types: vector field, stream plot, trajectory of a Cauchy problem.

\documentclass[11pt, margin=17pt]{standalone}
\usepackage{amsmath, amssymb}
\usepackage{ifthen}

\usepackage{tikz}
\usetikzlibrary{math, calc, arrows.meta}
\usetikzlibrary{decorations.markings}

\pgfkeys{/tikz/.cd,
x bound start/.store in=\xBoundStart,
x bound end/.store in=\xBoundEnd,
y bound start/.store in=\yBoundStart,
y bound end/.store in=\yBoundEnd
}

\makeatletter
\tikzset{%
bounds/.style n args={4}{%
evaluate={%
\xBoundStart = #1;
\xBoundEnd = #2;
\yBoundStart = #3;
\yBoundEnd = #4;
}
},
pics/vector field background/.style={% axes + grid for a vector field
code={%
\begin{scope}[every node/.style={scale=.8},
evaluate={%
\iIni = {int(\xBoundStart)};
\iEnd = {int(\xBoundEnd)};
\jIni = {int(\yBoundStart)};
\jEnd = {int(\yBoundEnd)};
}]
\draw[gray!50] (\xBoundStart, \yBoundStart)
grid (\xBoundEnd, \yBoundEnd);
\fill[black] (0, 0) circle (1pt);
\foreach \i in {\iIni, ..., \iEnd}{%
\path (\i, \yBoundStart) node[below=1ex] {$\i$};
}
\foreach \j in {\jIni, ..., \jEnd}{%
\path (\xBoundStart, \j) node[left=1ex] {$\j$};
}
\end{scope}
}
},
flow/.style={%  create a directional arrow along a trajectory
decoration={
markings,% switch on markings
mark=at position .5 with {\arrow[#1]{Latex}}
}, postaction=decorate
},
-flow/.style={%
decoration={
markings,% switch on markings
mark=at position .5 with {\arrowreversed[#1]{Latex}}
}, postaction=decorate
},
runge-kutta/.style={%  to be invoked before using ode-trajectory, ...
evaluate={%
int \N@ode, \i, \j;
\N@ode = 4;
for \i in {1, ..., \N@ode}{%
for \j in {1, ..., \N@ode}{%
\a@ode{\i,\j} = 0;
};
};
\a@ode{2,1} = 1/2; \a@ode{3,2} = 1/2; \a@ode{4,3} = 1;
\c@ode{1} = 0; \c@ode{2} = 1/2; \c@ode{3} = 1/2; \c@ode{4} = 1;
\b@ode{1} = 1/6; \b@ode{2} = 1/3; \b@ode{3} = 1/3; \b@ode{4} = 1/6;
}
},
% for ode in two dimensions
ode-trajectory/.style args={RHTx=#1, RHTy=#2, from=#3, to=#4, steps=#5}{%
insert path={coordinate (tmp) let \p1 = (tmp) in},
evaluate={%
real \showpointr, \testTmp;
\showpointr = 1;
int \steps@ode, \i, \j, \k, \s;
\t@ode{0} = #3;
\u@ode{0} = \x1*0.0352778;  % converts pt to cm
\v@ode{0} = \y1*0.0352778;
\steps@ode = int(#5);
\t@ode{\steps@ode} = #4;
\h = (\t@ode{\steps@ode} -\t@ode{0})/\steps@ode;
for \s in {1, 2, ..., \steps@ode}{
\k = int(\s -1);
\t@ode{0,1} = \t@ode{\k};
\x@ode{0,1} = \u@ode{\k};
\y@ode{0,1} = \v@ode{\k};
\qx@ode{0,1} = #1(\t@ode{0,1}, \x@ode{0,1}, \y@ode{0,1});
\qy@ode{0,1} = #2(\t@ode{0,1}, \x@ode{0,1}, \y@ode{0,1});
for \i in {2, ..., \N@ode}{%
\t@ode{0,\i} = \t@ode{\k} +\h*\c@ode{\i};
\x@ode{0,\i} = \u@ode{\k};
\y@ode{0,\i} = \v@ode{\k};
for \j in {1, ..., {int(\i -1)}}{%
\x@ode{0,\i} = \x@ode{0,\i} +\h*\a@ode{\i,\j}*\qx@ode{0,\j};
\y@ode{0,\i} = \y@ode{0,\i} +\h*\a@ode{\i,\j}*\qy@ode{0,\j};
};
\qx@ode{0,\i} = #1(\t@ode{0,\i}, \x@ode{0,\i}, \y@ode{0,\i});
\qy@ode{0,\i} = #2(\t@ode{0,\i}, \x@ode{0,\i}, \y@ode{0,\i});
};
\t@ode{\s} = \t@ode{\k} +\h;
\u@ode{\s} = \u@ode{\k};
\v@ode{\s} = \v@ode{\k};
for \i in {1, 2, ..., \N@ode}{%
\u@ode{\s} = \u@ode{\s} +\h*\b@ode{\i}*\qx@ode{0,\i};
\v@ode{\s} = \v@ode{\s} +\h*\b@ode{\i}*\qy@ode{0,\i};
};
};
\k = \steps@ode-1;
\testTmp = pow(pow(\u@ode{\steps@ode} -\u@ode{\k}, 2)
+pow(\v@ode{\steps@ode} -\v@ode{\k}, 2), .5);
},
insert path={%
\foreach \s in {1, ..., \steps@ode}{
-- (\u@ode{\s}, \v@ode{\s})
}
}
},
m vector field/.style n args={5}{% u / v / Nx / Ny / arrow scale
insert path = {%
\pgfextra{%
\let\tikz@mode@save=\tikz@mode
\let\tikz@options@save=\tikz@options%
}
},
evaluate={%
coordinate \tmp@P, \tmp@V;
integer \Nx, \Ny;
\Nx = #3;
\Ny = #4;
real \dx, \dy, \tmp, \uVF, \vVF, \nVF, \s@cst;
\s@cst = .49;
\dx = (\xBoundEnd -\xBoundStart)/\Nx;
\dy = (\yBoundEnd -\yBoundStart)/\Ny;
for \i in {1, ..., \Nx}{%
for \j in {1, ..., \Ny}{%
\uVF = #1(\xBoundStart -\dx/2 +\i*\dx, \yBoundStart -\dy/2 +\j*\dy);
\vVF = #2(\xBoundStart -\dx/2 +\i*\dx, \yBoundStart -\dy/2 +\j*\dy);
\nVF = pow(\uVF*\uVF +\vVF*\vVF, .5);
if  \nVF > \s@cst then {%
\uVF = \uVF/\nVF*.5;
\vVF = \vVF/\nVF*.5;
};
\tmp@P = (\xBoundStart +\i*\dx, \yBoundStart +\j*\dy);
\tmp@V = (#5*\uVF, #5*\vVF);
{
\draw[arrows={-Latex[width=2.5pt, length=2pt]}]
\pgfextra{\let\tikz@mode=\tikz@mode@save
\let\tikz@options=\tikz@options@save}
($(\tmp@P) -.5*(\tmp@V)$) -- ++(\tmp@V);
};
};
};
}
}
}
\newcommand{\streamPlot}[7]{% RHTx, RHTy, h, px, py, min, color
\tikzmath{%
integer \steps@ode, \i, \j, \k, \s\, flag;
real \h, \nor;
\flag = 0;
\s = 0;
\k = 0;
\h = #3;
\u@ode{0} = #4;
\v@ode{0} = #5;
\t@ode{0} = 0;
},
\whiledo{\equal{\flag}{0}}{%
\tikzmath{
\s = \s +1;
\t@ode{0,1} = 0;
\x@ode{0,1} = \u@ode{\k};
\y@ode{0,1} = \v@ode{\k};
\qx@ode{0,1} = #1(\t@ode{0,1}, \x@ode{0,1}, \y@ode{0,1});
\qy@ode{0,1} = #2(\t@ode{0,1}, \x@ode{0,1}, \y@ode{0,1});
for \i in {2, ..., \N@ode}{%
\t@ode{0,\i} = \t@ode{\k} +\h*\c@ode{\i};
\x@ode{0,\i} = \u@ode{\k};
\y@ode{0,\i} = \v@ode{\k};
for \j in {1, ..., {int(\i -1)}}{%
\x@ode{0,\i} = \x@ode{0,\i} +\h*\a@ode{\i,\j}*\qx@ode{0,\j};
\y@ode{0,\i} = \y@ode{0,\i} +\h*\a@ode{\i,\j}*\qy@ode{0,\j};
};
\qx@ode{0,\i} = #1(\t@ode{0,\i}, \x@ode{0,\i}, \y@ode{0,\i});
\qy@ode{0,\i} = #2(\t@ode{0,\i}, \x@ode{0,\i}, \y@ode{0,\i});
};
\t@ode{\s} = \t@ode{\k} +\h;
\u@ode{\s} = \u@ode{\k};
\v@ode{\s} = \v@ode{\k};
for \i in {1, 2, ..., \N@ode}{%
\u@ode{\s} = \u@ode{\s} +\h*\b@ode{\i}*\qx@ode{0,\i};
\v@ode{\s} = \v@ode{\s} +\h*\b@ode{\i}*\qy@ode{0,\i};
};
\nor = pow(%
pow(\u@ode{\s} -\u@ode{\k}, 2) +pow(\v@ode{\s} -\v@ode{\k}, 2),
.5);
if \nor<#6 then {\flag = 1;};
if \u@ode{\s}<\xBoundStart then {\flag = 1;};
if \u@ode{\s}>\xBoundEnd then {\flag = 1;};
if \v@ode{\s}<\yBoundStart then {\flag = 1;};
if \v@ode{\s}>\yBoundEnd then {\flag = 1;};
\k = \k +1;
}
},
\tikzmath{%
if \s>1 then {%
if #3>0 then {%
{%
\draw[#7, flow={#7}]  (\u@ode{0}, \v@ode{0})
\foreach \i [parse=true] in {1, ..., \s-1}{
-- (\u@ode{\i}, \v@ode{\i})
};
};
} else {%
{%
\draw[#7, -flow={#7}]  (\u@ode{0}, \v@ode{0})
\foreach \i [parse=true] in {1, ..., \s-1}{
-- (\u@ode{\i}, \v@ode{\i})
};
};
};
};
}
}
\makeatother

\tikzset{
pics/stream plot by name/.style args={RHTx=#1, RHTy=#2, dT=#3, pName=#4,
pNumber=#5, min=#6, color=#7}{%
code={
\foreach \i in {1, ..., #5}{%
\tikzmath{%
coordinate \P;
\P = (#4-\i);
\px = \Px*1pt/1cm;  % converts pt to cm
\py = \Py*1pt/1cm;
}
% \draw (\px, \py) circle (2pt);
\streamPlot{#1}{#2}{#3}{\px}{\py}{#6}{#7};
\streamPlot{#1}{#2}{-#3}{\px}{\py}{#6}{#7};
}
}
},
pics/stream plot TB/.style args={%
RHTx=#1, RHTy=#2, dT=#3, uSteps=#4, min=#5, color=#6}{%
code={
\tikzmath{%
real \px, \py, \dx;
\dx = (\xBoundEnd -\xBoundStart)/#4;
}
\foreach \i in {0, ..., #4}{%
\tikzmath{%
\px = \xBoundStart +\i*\dx;
\py = \yBoundStart;
}
\streamPlot{#1}{#2}{#3}{\px}{\py}{#5}{#6};
\streamPlot{#1}{#2}{-#3}{\px}{\py}{#5}{#6};
}
\foreach \i in {0, ..., #4}{%
\tikzmath{%
\px = \xBoundStart +\i*\dx;
\py = \yBoundEnd;
}
\streamPlot{#1}{#2}{#3}{\px}{\py}{#5}{#6};
\streamPlot{#1}{#2}{-#3}{\px}{\py}{#5}{#6};
}
}
}
}

\begin{document}

%\iffalse
\begin{tikzpicture}[every node/.style={%
right, align=left, inner sep=3ex
}]
\path (-1.3, 0) node {$U(x, y)$}
(0, 0) node[] {$\displaystyle = y\left[ -2\,\frac{\ln(\cosh x)}{x^2} +\frac{x -\tanh x}{x^2}\,y +\frac{\tanh x}{x}\,(2 -y\,\tanh x) \right]$}
++(0, -1) node[] {$\displaystyle = y\left[ 1 +\frac{2}{9}\,x^4 -\bigg(\frac{2}{3}\,x -\frac{8}{15}\,x^3\bigg)y \right] +[5]$};

\path (-1.3, -3) node {$V(x, y)$}
(0, -3) node[] {$\displaystyle = 2\,\frac{\ln(\cosh x)}{x} +2\,\frac{\tanh x -x}{x}\,y$}
++(0, -1) node[] {$\displaystyle = x\left[ 1 -\frac{1}{3}\,x^2 -\bigg( \frac{2}{3}\,x -\frac{4}{15}\,x^3 \bigg)y \right] +[5]$};
\end{tikzpicture}
%\fi

\tikzmath{%
function Ux(\t, \u, \v) {%
real \a, \b, \c;
if abs(\u)>.05 then {%
\a = -2*ln(cosh(\u))/(\u*\u);
\b = (\u -tanh(\u))/(\u*\u)*\v;
\c = (tanh(\u)/\u)*(2 - tanh(\u)*\v);
} else {%
\a = 1 -2/3*\u*\v;
\b = 8/15*pow(\u, 3)*\v;
\c = 2/9*pow(\u, 4);
};
return {(\a +\b +\c)*\v};
};
function Uy(\t, \u, \v) {%
real \a, \b;
if abs(\u)>.05 then {%
\a = ln(cosh(\u))/\u;
\b = (tanh(\u) -\u)/\u;
return {2*\a +2*\b*\v};
} else {
\a = 1 -pow(\u, 2)/3;
\b = -(2/3)*\u +(4/15)*pow(\u, 3);
return {\u*(\a +\b*\v)};
};
};
function fUy(\u) {%
real \a, \b;
\a = ln(cosh(\u))/\u;
\b = (\u -tanh(\u))/\u;
return \a/\b;
};
function UxVF(\u, \v) {%
real \a, \b, \c;
if abs(\u)>.05 then {%
\a = -2*ln(cosh(\u))/(\u*\u);
\b = (\u -tanh(\u))/(\u*\u)*\v;
\c = (tanh(\u)/\u)*(2 - tanh(\u)*\v);
} else {%
\a = 1 -2/3*\u*\v;
\b = 8/15*pow(\u, 3)*\v;
\c = 2/9*pow(\u, 4);
};
return {(\a +\b +\c)*\v};
};
function fUx(\u) {%
real \a, \b;
\a = 2*(tanh(\u)/\u) -2*ln(cosh(\u))/(\u*\u);
\b = -(\u -tanh(\u))/(\u*\u) +(tanh(\u)/\u)*tanh(\u);
return \a/\b;
};
function UyVF(\u, \v) {%
real \a, \b;
if abs(\u)>.05 then {%
\a = ln(cosh(\u))/\u;
\b = (tanh(\u) -\u)/\u;
return {2*\a +2*\b*\v};
} else {
\a = 1 -pow(\u, 2)/3;
\b = -(2/3)*\u +(4/15)*pow(\u, 3);
return {\u*(\a +\b*\v)};
};
};
}

\iffalse
\begin{tikzpicture}
\draw[gray!50, very thin] (-3.5, -3.25) grid (3.5, 3.25);
\draw[->] (-3.5, 0) -- (3.5, 0) node[above right] {$x$};
\draw[->] (0,-3.5) -- (0, 3.5) node[left] {$y$};

\draw[blue!70!black, thick, variable=\t, domain=.5:3.25, samples=200]
plot (\t, {fUx(\t)});
\draw[blue!70!black, thick] (-3.25, 0) -- (3.25, 0);
\draw[blue!70!black, thick, variable=\t, domain=-3.25:-.5, samples=200]
plot (\t, {fUx(\t)});
\draw[red, thin, variable=\t, domain=.5:3.25, samples=200]
plot (\t, {fUy(\t)});
\draw[red] (-0, -3.25) -- (0, 3.25);
\draw[red, thin, variable=\t, domain=-3.25:-.5, samples=200]
plot (\t, {fUy(\t)});
\end{tikzpicture}
\fi

\iffalse
\begin{tikzpicture}[runge-kutta, bounds={-4}{4}{-4}{4}, scale=1.2]
\path pic[scale=1.2] {vector field background};

\draw[blue!70!black, thick,
variable=\t, domain=.35:\xBoundEnd, samples=200]
plot (\t, {fUx(\t)});
\draw[blue!70!black, thick] (\xBoundStart, 0) -- (\xBoundEnd, 0);
\draw[blue!70!black, thick,
variable=\t, domain=\xBoundStart:-.35, samples=200]
plot (\t, {fUx(\t)});
\draw[red, thin, variable=\t, domain=.35:\xBoundEnd, samples=200]
plot (\t, {fUy(\t)});
\draw[red] (-0, \xBoundStart) -- (0, \xBoundEnd);
\draw[red, thin, variable=\t, domain=\xBoundStart:-.35, samples=200]
plot (\t, {fUy(\t)});
\draw[green!50!black, thick, m vector field={UxVF}{UyVF}{35}{35}{.4}];
\end{tikzpicture}
\fi

\iffalse
\begin{tikzpicture}[runge-kutta, bounds={-4}{4}{-3.5}{3.5}]
\path pic {vector field background};

% A points
\path (0, 1.5) coordinate (A-1);
\path (0, -1.5) coordinate (A-2);
\path (2.2, \yBoundEnd) coordinate (A-3);
\path (-2.2, \yBoundStart) coordinate (A-4);

% B points
\tikzmath{
integer \N;
\N = 4;
}
\foreach \i in {1, ..., \N}{%
\path (\i*360/\N: .5) coordinate (B-\i);
}

\path pic {stream plot by name={RHTx=Ux, RHTy=Uy, dT=.05,
pName=A, pNumber=4, min=.005, color=blue!60!black}};
\path pic {stream plot by name={RHTx=Ux, RHTy=Uy, dT=.05,
pName=B, pNumber=\N, min=.005, color=blue!60!black}};

\path pic {stream plot TB={RHTx=Ux, RHTy=Uy, dT=.05,
uSteps=11, min=.005, color=blue!60!black}};

\draw[red, very thick] (2.7, -3)
[ode-trajectory={RHTx=Ux, RHTy=Uy, from=0, to=15, steps=1500}];
\end{tikzpicture}
\fi

\end{document}


The level curves of the U component of the vector field. The red curve from Fig. 1 is the black trajectory here!