# What's the equation behind this strange curve?

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... but that gave me no clue about the very nature of this curve. @user21820: the data can be found on pastebin: here

This sequence of points seems to grow with a slope of exactly 1 until the slope decreases to zero quite sharply... 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... This is numerically verified by the look of the second derivative: the third derivative, if of any help? (note that the grey line here should not to be considered) Well... Hope you can help me !

NB: I can provide the raw data if necessary, but I felt this would render the question unreadable (2000 values...)

edit: @Ron: The plot of $f(x) = 500 * \tanh(x/500)$ in red; and $f(x) = 500 * \sqrt(\tanh((x/500)^2))$ in green: edit2: @user21280: This is also very close but there seems to be a little disagreement, again around the bending point. edit3: 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 compartment. 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 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.

edit4: 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 <- B0

for(i in 1:niequ){

Fab <- .001*A
Fba <- 0.1*B*(1-A/Amax)
B <- B - Fba + Fab
A <- A  + Fba - Fab

}
rec <- 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...

• Where does the data come from? I don't think one always should expect that collected data can be represented by some very simple function. Of course one could come up with some "monster" that will fit it pretty good, but as I read it this is not what you look for? – mickep Dec 29 '14 at 12:18
• @mickep: You are right, this is important and I will try to explain how I came across this curve in a few minutes. – Rodolphe Dec 29 '14 at 12:55
• Could you provide the numerical data? Use pastebin or something else. – XYZT Dec 29 '14 at 13:09
• State your source of the data.. The data fits too perfectly to a curve to be from real-world measurements. – user21820 Dec 29 '14 at 14:04
• Equilibrium occurs at $F_{ab}=F_{ba}$. Solve for $a$ and $b$ subject to the conservation of stock, $a+b=a_0+b_0=x$. You get $$a= \frac{1}{2} \left(x+505-\sqrt{x^2-990 x+255025}\right),\\b= \frac{1}{2} \left(x-505+\sqrt{x^2-990 x+255025}\right).$$ – user856 Dec 29 '14 at 15:08

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_{eq} = \frac12\left(x+1.01A_{max}\pm\sqrt{x^2-1.98A_{max}x+1.0201A_{max}^2}\right)\\ A_{eq} = \frac12\left(x+505\pm\sqrt{x^2-990x+255025}\right)\\ \text{(probably the negative square-root)}\\ y=x-A_{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)}$$

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.

• Wait a minute, neither does your function fit. Your function fits the first derivative very well according to my experience of sigmoid functions. – user21820 Dec 29 '14 at 12:29
• @user21820: a plot of the second derivative reveals a bump below the axis as illustrated. The only difference is that the second derivative has a nonzero slope at zero, but I feel this is a minor point. The third derivative has the same general behavior as pictured. The OP may wish to introduce additional parameters to improve the fit, but this is a good start. – Ron Gordon Dec 29 '14 at 12:31
• @RonGordon: Your were closer than i never was with this equation, however it does not seem to "bend fast enough" when compared with the actual data. I added a plot of your equation for comparison. – Rodolphe Dec 29 '14 at 12:37
• @Rodolphe: see my edit, which accounts for the bending. – Ron Gordon Dec 29 '14 at 12:40
• @RonGordon: edited the graph with adjusted equation in green. – Rodolphe Dec 29 '14 at 12:47

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!

• This function would not have the second derivative profile shown. – Ron Gordon Dec 29 '14 at 12:23
• @RonGordon: Oh yes.. Let me fix the differential equation then.. – user21820 Dec 29 '14 at 12:24
• @RonGordon: See this now. =) – user21820 Dec 29 '14 at 12:44
• @Rodolphe: How does my guess compare? You might as well upload the data somewhere where we can get it... – user21820 Dec 29 '14 at 12:48
• @user21820: Yes I would like to upload it but I don't know how nor where, nor how to format it for you to use it (I use R programming language). As for your guess, I am currently working on it. back in a few minutes :-) – Rodolphe Dec 29 '14 at 12:51