# Square root of a function "misbehaves" near the x-axis

I'm trying to plot height lines of a function of the form $H(x,y)=y^2+u(x)$, by taking multiple plots of $\pm\sqrt{c-u(x)}$. The problem is that for some values $c$ the domain is tricky and not easily computable. I tried to bypass this problem by plotting $\pm\max{\sqrt{c-f(x)},0}$ instead (which is "absorbed" by the x-axis) but it results in those visible, annoying "peaks". I need a clever idea to get rid of them.

My code:

\documentclass{article}
\usepackage{tikz}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}[
declare function = {u(\x) = 0.5*pow(\x-2,3)-\x+4;}
]
\draw[->] (-0.5,0)--(4.5,0);
\draw[->] (0,-2.2)--(0,2.2);

\clip (-0.5,-2.2) rectangle (4,2.2);
\def\samp{200}
\foreach \c in {0,0.8,...,3.6}
{
\draw plot[domain=-0.5:4, samples=\samp] (\x, {sqrt(max(\c-u(\x),0)))});
\draw plot[domain=-0.5:4, samples=\samp] (\x, {-sqrt(max(\c-u(\x),0)))});
}
\end{tikzpicture}
\end{document}


Result:

• Welcome to TeX.SX! Please make your code compilable (if possible), or at least complete it with \documentclass{...}, the required \usepackage's, \begin{document}, and \end{document}. That may seem tedious to you, but think of the extra work it represents for TeX.SX users willing to give you a hand. Help them help you: remove that one hurdle between you and a solution to your problem. Commented Sep 7, 2021 at 19:32
• Hope it's better now Commented Sep 7, 2021 at 19:37
• Just wondering: is it more a mathematical question (e.g. providing a better algorithm in certain regions), or more a Latex resp. tikz related one? Commented Sep 7, 2021 at 19:41
• @MS-SPO Thanks! It is an example of a phase portrait, which demonstrates conservation of energy in a one-dimensional mechanical system. For further information you can check out Mathematical Methods of Classical Mechanics by VI Arnold, Chapter 2 Commented Sep 7, 2021 at 20:01
• Thanks @35T41 for the pointer :) Exciting subject. // Glad my little trick helped :) Commented Sep 7, 2021 at 20:03

Looks like all you need to do is increasing your sample size. Because this takes a bit longer, perhaps you can do it only near y==0 , i.e. have a course (200 samples) and a limited fine region (2000 samples). Or you just spend the waiting time.

\documentclass{article}
\usepackage{tikz}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}[
declare function = {u(\x) = 0.5*pow(\x-2,3)-\x+4;}
]
\draw[->] (-0.5,0)--(4.5,0);
\draw[->] (0,-2.2)--(0,2.2);

\clip (-0.5,-2.2) rectangle (4,2.2);
\def\samp{2000}% <<--
\foreach \c in {0,0.8,...,3.6}
{
\draw plot[domain=-0.5:4, samples=\samp] (\x, {sqrt(max(\c-u(\x),0)))});
\draw plot[domain=-0.5:4, samples=\samp] (\x, {-sqrt(max(\c-u(\x),0)))});
}
\end{tikzpicture}
\end{document}


P.S.: If you have more drawings like this and are worried to combine them into a larger Latex document, here is a way to do it:

• create the drawings by at least one separate Latex-document
• save them as .eps

See Chapter "7.1 Export to pdf/eps" in the big tikz manual (pgfplots.pdf).

• Wouldn't such a large sample size be "too heavy" to compute, somehow? Should I ever even worry about something like that? Commented Sep 7, 2021 at 19:58
• Well, my notebook ran your code, and the increased sample size between the two comments ... :) So just try it, I'd say ... Commented Sep 7, 2021 at 19:59

For comparison purpose, compile with Asymptote.

size(7cm);

import contour;

// H(x,y)=y^2+u(x) with u(x)=0.5*(x-2)^3-x+4;
real f(real x, real y){ return y^2+0.5*(x-2)^3-x+4;}

// Array c={0,0.8,...,3.6}
real[] c;
for (int i=0; i<=36; i=i+8) {  c.push(i/10); }

draw((-0.5,0)--(4.5,0),Arrow);
draw((0,-2.2)--(0,2.2),Arrow);

// one-dimensional pen array
pen[] pp={red, green, blue, gray};
pp.cyclic=true; // See 6.12 Arrays

// two-dimensional guide array
guide[][] g=contour(f,(-0.5,-2.2),(4,2.2),c,400); // 400 is ngraph!
guide[] G=concat(... g); // See 6.12 Arrays

for (int i=0; i<G.length; ++i)
draw(G[i],pp[i]);


The default resolution, ngraph x ngraph (here ngraph defaults to 100) can be increased for greater accuracy.

See also 8.35 contour (Main documentation) for more information.

You can handle with each guide manually like this:

size(7cm);

import contour;

// H(x,y)=y^2+u(x) with u(x)=0.5*(x-2)^3-x+4;
real f(real x, real y){ return y^2+0.5*(x-2)^3-x+4;}

// Array c={0,0.8,...,3.6}
real[] c;
for (int i=0; i<=36; i=i+8) {  c.push(i/10); }

draw((-0.5,0)--(4.5,0),Arrow);
draw((0,-2.2)--(0,2.2),Arrow);

// two-dimensional guide array
guide[][] g=contour(f,(-0.5,-2.2),(4,2.2),c,400); // 400 is ngraph!
guide[] G=concat(... g); // See 6.12 Arrays

write(G.length); // Outputs: 7

draw(Label("$0$",EndPoint),G[0],red,Arrow);
draw(Label("$1$",BeginPoint),G[1],green,Arrow);
draw(Label("$2$",BeginPoint),G[2],blue+dashdotted,Arrow);
draw(Label("$3$",BeginPoint),G[3],red+dashed,Arrow);
draw(Label("$4$",EndPoint),G[4],gray,Arrow);
draw(Label("$5$",BeginPoint),G[5],magenta,Arrow);
draw(Label("$6$",EndPoint),G[6],cyan,Arrow);


We can do this with TikZ but using less samples too. But in this case we need to help TikZ to do the task.

Just take a numerical calculus software of your choice and find when the polynomial has multiple roots (Cubic equation) and then find its roots. Next, plot only the positive parts.

I think that in this case TikZ probably could handle the maths (Cardano and friends or Newton-Raphson?) but the time we saved in the samples will go to the maths so it doesn't worth the effort.

I made my example with at most 100+100 samples per plot. This number could be reduced if we split the curves in more pieces and put more samples near the roots.

This is my code (I used Python to find the roots):

\documentclass[tikz,border=2mm]{standalone}

\tikzset
{% abs is to prevent rounding errors causing a 'false negative' number
declare function={f(\x,\y)=sqrt(abs(-0.5*\x*\x*\x+3*\x*\x-5*\x+\y));},
my color/.style={blue, opacity={-0.4*abs(2-#1)+1.2}}
}

\begin{document}
\begin{tikzpicture}
\draw[-latex] (-0.5,0)-- (4.5,0);
\draw[-latex] (0,-3)  -- (0,3);
\def\samp{100}
% one piece
\foreach\d/\r in {0/0,0.4/0.0842,0.8/0.1786,1.2/0.2871,2.8/3.7129,3.2/3.8214,3.6/3.9158,4/4} \foreach\j in {-1,1}
\draw[my color=\d] plot[domain=-0.5:\r, samples=\samp] (\x, {\j*f(\x,\d)}) -- (\r,0);
% two pieces
\foreach\d/\r/\s/\t in {1.6/0.4171/2.4437/3.1392,2/0.5858/2/3.4142,2.4/0.8608/1.5563/3.5829,2.5443/1.1835/1.1835/3.633} \foreach\j in {-1,1}
\draw[my color=\d] plot[domain=-0.5:\r, samples=\samp] (\x, {\j*f(\x,\d)}) -- (\r,0)
(\s,0) -- plot[domain=\s:\t,   samples=\samp] (\x, {\j*f(\x,\d)}) -- (\t,0);
\fill[blue] (2.8165,0) circle (1pt);
\end{tikzpicture}
\end{document}


And the output: