What’s the equation behind this strange curve?

Question

I have come across a strange kind of curve, which I can’t imagine what is the equation behind. I’ve tried a variety of approaches, from numerically calculating the derivative and second (and so on) derivative$\ldots$ but that gave me no clue about the very nature of this curve.

@user21820: The data can be found on pastebin [here].(http://pastebin.com/UbhmETsd)

This sequence of points seems to grow with a slope of exactly 1 until the slope decreases to zero quite sharply$\ldots$ I numerically calculated the first three derivatives of this curve:

The first derivative looks like:

Interestingly (I think), the first derivative’s “turning point” is at $x=500$, just like the $y$ maxima of the function I a looking for$\ldots$ This is numerically verified by the looks of the second derivative:

Is the third derivative of any help? (Note that the grey line here should not to be considered.)

Well$\ldots$ Hope you can help me!

N.B.: I can provide the raw data if necessary, but I felt this would render the question unreadable (2000 values$\ldots$)

Edit 1: @Ron: The plot of $f(x) = 500 \tanh(x/500)$ in red; and $f(x) = 500 \sqrt \tanh{((x/500)^2)}$ in green:

Edit 2: @user21280: This is also very close but there seems to be a little discrepancie, again around the bending point.

Edit 3: The story behind this mysterious curve: In a closed system with two compartments A and B, exchanges of some element E occurs between these two compartments. A is finite while B is not. Flux from A to B (Fab) in one unit of time is a fraction of stock of E in A (Fab = 0.001 * A). Flux from B to A is a fraction of stock of E in B times the fraction of free “space” in A (Fba = 0.1 * B * (Amax-A)/Amax).

Considering, at the start of simulation, that A = 0 and B = x, for x = 1:2000, what is the value y of B at the equilibrium? (given Amax = 500)

When B0 gets high enough, A is saturated and B_eq (y) increases as a linear function of B0. Hence the idea of substracting B0 from B_eq. Which gives away my mysterious curve.

Edit 4: The reproducible R code behind:

Amax <- 500 # max stock of element in solid phase
A0 <- 0 # stock of element in solid phase (A) at time zero
niequ <- 1000   # time allowed for reaching equilibrium between A and B
rec <- c(0)
for(B0 in 1:2000){  # B0 = stock of element in B at time zero
    A <- A0
B &lt;- B0

for(i in 1:niequ){

    Fab &lt;- .001*A
    Fba &lt;- 0.1*B*(1-A/Amax)
    B &lt;- B - Fba + Fab
    A &lt;- A  + Fba - Fab

}
rec &lt;- c(rec,mean(B))

}
rec <- rec[-1]
plot(rec~c(1:length(rec)))  # B at equilibrium as a function of B at time zero
rec2 <- (rec-c(1:length(rec)))  # substraction of linear increase
rec2 <- -rec2 # for working with positive values (easier for me)
plot(x=c(1:length(rec)), y=rec2)    # the mysterious curve...

Answer 1

This to me seems like a hyperbolic tangent function:

$$f(x) = a \tanh{b x}$$

A quick inspection reveals that $a$ should be about $500$. Perhaps the slope at $x=0$ being about $1.0$ would imply that $b=1/500$. Obviously, these numbers should be taken with a grain of salt but at least represent a starting point for the analysis. The functional form seems to at least match somewhat the derivative profiles.

EDIT

A comparison to the data showed that $f$ did not bend fast enough. Another possibility is

$$g(x) = 500 \sqrt{\tanh{\left [\left (\frac{x}{500} \right )^2\right ]}}$$

A quick inspection shows that there is a better fit to the bending, although the behavior of the second derivative is a bit different.

Answer 2

I don't know a good reason for the data, but it looks like $x \mapsto \int_0^x \left( \frac{1}{2} - \frac{1}{2}\tanh(\frac{t-500}{100}) \right)\ dt$. It matches all the graphs perfectly as far as I can tell.

[Edit: After Rodolphe finally told the origin of the data...]

It can be seen that it is just a differential equation $\frac{da}{dt} = -pa+q(1-\frac{a}{m})(c-a)$ where $a$ is the amount in $A$, $c$ is the total amount in both, $m$ is the capacity of $A$ and $p,q$ are fixed parameters. The differential equation is separable and just requires integrating the reciprocal of a quadratic function, which is easy. I'm too lazy to do it though.

[Edit: After Rodolphe posted the code, it is clear that I misunderstood the process...]

Michael did it so I don't have to. Great!

Answer 3

At equilibrium, $$F_{ab}=F_{ba}\\ 0.001A = 0.1B(1-A/A_{\max.})\\ 0.001A = 0.1(x-A)(1-A/A_{\max.})\\ 0.01AA_{\max.}=(A-x)(A-A_{\max.})\\ A^2-(x+1.01A_{\max.})A+xA_{\max.}=0\\ A_\text{eq.} = \frac12\left(x+1.01A_{\max.}\pm\sqrt{x^2-1.98A_{\max.}x+1.0201A_{\max.}^2}\right)\\ A_\text{eq.} = \frac12\left(x+505\pm\sqrt{x^2-990x+255025}\right)\\ \text{(probably the negative square-root)}\\ y=x-A_\text{eq.} = \frac12\left(x-1.01A_{\max.}\mp\sqrt{x^2-1.98A_{\max.}x+1.0201A_{\max.}^2}\right)\\ y=\frac12\left(x-505\mp\sqrt{x^2-990x+255025}\right)\\ \text{(probably the positive square-root)} $$

Answer 4

Here is IMHO a slightly more direct approach.

Indeed, this curve with 2 neat asymptotes invites to consider it as a (branch of) hyperbola having a horizontal asymptote with approx. equation $y-500=0$ and a slant asymptote with approx. equation $y-x=0$. The common equation of all hyperbolas with such asymptotes is:

$$(y-500)(y-x)=k$$

where $k$ is constant (where the LHS is the product of the LHS of the asymptotes).

Forcing the curve to pass by a certain point $A(x_0,y_0)$ clearly fixes the value of the constant. If we fix (approx.) $A(1,0.9900796)$ (according to the given data), we obtain the following equation:

$$(y-500)(y-x)=k_0\tag{1}$$

If we expand (1) as a polynomial in $y$ with parameter $x$, we get a quadratic equation

$$y^2-(x+500)y+(500x-k_0)=0$$

whose solutions are:

$$y=\tfrac12 \left(x+500 \pm \sqrt{(x+500)^2-(500x-k_0)} \right)$$

finding back the solution given by others.

Remark: in fact the $\pm$ sign accounts for the two branches of the hyperbola. Clearly, we need to replace it by a minus sign in order to get the desired branch.