Janky pgfplot behaviour with arcsinh function

Question

I'm trying to draw a titration curve using pgfplot. Based off this paper, there is an equation for the titration curve:

I've attempted to implement this within pgfplots, and I get a very strange result:

I have no idea why this would occur; even attempting to recreate the exact formula I used within desmos did not work.

For reference, here is the code, and also a link to an overleaf project.

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\pgfkeys{
    /pgf/declare function={
        arcsinh(\x) = ln(\x + sqrt(\x^2+1));
    },
    /pgf/declare function={
        Va = 0.025;
        Ma = 0.1;
        Mb = 0.1;
        V(\x) = \x / 1000;
        Kw = 1*10^(-14);
        p(\o) = -ln(\o)/ln(10);
    }
}
\begin{center}
\begin{tikzpicture}
    \begin{axis}[
        xlabel = {Solution Added (mL)},
        ylabel = {pH},
        ymin=0,
        ymax=14,
        ytick distance=7,
        xtick distance=10,
    ]
    \addplot[%
        samples=100,
        color=red,
        domain=0:50,
    ]{%
    7 + 1/ln(10) * arcsinh( 1/(2sqrt(Kw))    (MbV(x) - MaVa) / (Va + V(x)) )
    };
    \end{axis}
\end{tikzpicture}
\end{center}
\end{document}

~~endnote: if there's a way I can have variables e.g. V_a within the pgfplot function, that'd be much nicer than having constants throughout.~~ thanks Torjorn

@egreg If you're referring to the division by ln(10) int the definition of arcsinh, it's just been shifted from outside to inside (if you look at the reference function it's outside). I don't think it makes any difference, and I've tried having it both ways. It just makes it easier to read imo, even if it's not strictly an arcsinh function anymore. — IllustriousMagenta, Sep 02 '20 at 14:33
@Modelmat That's what I wanted to underline: you know, a mathematician's hyperbolic sine uses e. ;-) — egreg, Sep 02 '20 at 14:58
Nope, although I do see a large number of NOTE: coordinate (1Y9.999739e1],3Y0.0e0]) has been dropped because it is unboun ded (in y). (see also unbounded coords=jump). — IllustriousMagenta, Sep 03 '20 at 03:44
I've updated the question with some cleaner code and also a link to an overleaf project demonstrating this. Any help would be appreciated! — IllustriousMagenta, Sep 03 '20 at 12:56

Stefan Pinnow · Accepted Answer · 2020-09-04T09:16:52.023

One gets this strange result because TeXs limits are reached -- as can also be seen by the dropped coordinates -- and thus the "zigzag" results from precision limits (red line). If you use either gnuplot (green dots) or Lua (blue line) as calculation engine, it works as expected. Of course for the Lua solution you have to use LuaLaTeX as TeX engine.

Sidenote:
If you want to avoid to use such a high number of samples, consider reformulating the eqution to make use of non-linear spacing. For that see e.g.

% used PGFPlots v1.17
\documentclass[border={5pt}]{standalone}
\usepackage{pgfplots}
    % use this `compat` level or higher to use LUA backend for calculation (if possible)
    \pgfplotsset{compat=1.12}
\begin{document}
    % for gnuplot solution
    \newcommand*{\Kw}{1e-14}
    \newcommand*{\Ma}{0.1}
    \newcommand*{\Mb}{0.1}
    \newcommand*{\Va}{0.025}
\begin{tikzpicture}[
    % from https://tex.stackexchange.com/q/144778
    /pgf/declare function={
        % for LUA solution
        arcsinh(\x) = ln(\x + sqrt(\x^2+1));
        Kw = 1e-14;
        Ma = 0.1;
        Mb = 0.1;
        Va = 0.025;
        V(\x) = \x / 1000;
        p(\o) = -0.5*ln(\o)/ln(10);
    },
]
    \begin{axis}[
        xlabel={Solution Added (mL)},
        ylabel={pH},
        domain=0:50,
        samples=201,
    ]
        % gnuplot
        \addplot+ [ultra thick,green,mark size=1pt,only marks,opacity=0.5] gnuplot
            {-0.5*log10(\Kw) + asinh(1/(2*sqrt(\Kw)) * (\Ma*\Va - \Mb*x/1000)/(\Va + x/1000) )/log(10)};
        % Lua(TeX)
        \addplot [thick,blue]
            {p(Kw) + arcsinh(1/(2*sqrt(Kw)) * (Ma*Va - Mb*V(x))/(Va + V(x)) )/ln(10)};
        % TeX
        \addplot [red,opacity=0.75]
            {p(\Kw) + arcsinh(1/(2*sqrt(Kw)) * (Ma*Va - Mb*V(x))/(Va + V(x)) )/ln(10)};
    \end{axis}
\end{tikzpicture}
\end{document}

Thanks heaps! I'm unsure how replacing -0.5log10(K_w) with 7 would make the curve "incorrect"? K_w will always be what it is, and thus 0.5*pK_w will always be 7. — IllustriousMagenta, Sep 03 '20 at 20:11
You are right. I got distracted from just not using your definition of p ... I removed that comment. — Stefan Pinnow, Sep 04 '20 at 09:15

Janky pgfplot behaviour with arcsinh function

1 Answers1

Linked