Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have discrete collection of data points (10 to 10^4). I want to describe them by a continuous function and find a x value z, where areas S1 and S2 bounded by {x, f[x], yBot, yTop} are equal. Furthermore, I think it is called Matano plane.

There is 1 important special case: if function values dips beyond {yBot, yTop}, then that area should be added to the other area.

Shaded regions on this plot should be equal (S1=blue,S2=yellow):

enter image description here

Example data for this plot: n = 3; and z is the divider what I am looking for

rc = RandomChoice[{Sech[2 x - 0.5]/1.5 + 1.5, 2}]
data = Transpose[{
   # + RandomReal[] & /@ Range[-10, 30, 2],
   Tanh[#] + rc /. x -> # & /@ Range[-4, 4, 0.4]}]
f = GetRLine3[{data}, 3][x];
{xLimits, yLimits} = {
  {Min@#1, Max@#1}, 
  {Mean@#2[[1 ;; n]], Mean@#2[[-n ;; -1]]}} & @@ Transpose[data]

Where GetRLine3 is:

GetRLine3[MMStdata_, IO_: 1][x_: x] :=ListInterpolation[#, InterpolationOrder -> IO,
  Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);

I can write a adhoc solution that finds z, but then it takes about 5.953s to find.

z = step = xLimits[[1]] + (xLimits[[2]] - xLimits[[1]])/2;
While[{part1, part2} = Flatten[{Integrate[f - yLimits[[1]], {x, xLimits[[1]], z}],
     -Integrate[(f - yLimits[[2]]), {x, z, xLimits[[2]]}]}] // N; part1 != part2,
 If[part1 > part2, step /= 2; z -= step, step /= 2; z += step]]

In this case z=8.075996010318093.

I can't seem to be able to use solve to find z. I have tried other regression curves

  • Fit with Table[x^k,{k,0,100}]
  • NonlinearModelFit with a/2 Erf[(x-\[Mu])/(Sqrt[2]\[Sigma])]+y0

and they are about 10x faster, but they did't work for both small and large sets of data. Finding accurate model or set of functions seems to be beyond me.

Is there a better or faster way to solve for z in range xLimits.

share|improve this question
up vote 7 down vote accepted

Maybe I'm missing something but isn't the value for z just the solution of the equation

(z - xLimits[[1]])*yLimits[[1]] + (xLimits[[2]] - z)*yLimits[[2]] == area

where area is the total area under the graph, i.e.

area = NIntegrate[GetRLine3[{data}, 3][x], {x, xLimits[[1]], xLimits[[2]]}][[1]]

Therefore, z is equal to

z = (xLimits[[2]]*yLimits[[2]] - xLimits[[1]]*yLimits[[1]] - area)/
      (yLimits[[2]] - yLimits[[1]])
share|improve this answer
This seems to be correct and it's x24 times faster. – Margus Apr 28 '12 at 8:21

For a start it seems to me that you can use NIntegrate:

z = step = xLimits[[1]] + (xLimits[[2]] - xLimits[[1]])/2;

  {part1, part2} = 
    NIntegrate[f - yLimits[[1]], {x, xLimits[[1]], z}],
   -NIntegrate[(f - yLimits[[2]]), {x, z, xLimits[[2]]}]
  part1 != part2, 
  If[part1 > part2, step /= 2; z -= step, step /= 2; z += step]
] // Timing

{0.14, Null}

share|improve this answer
You are right, change form Integrate + //N to NIntegrate does give 2.5x speed up in time. – Margus Apr 27 '12 at 21:58

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.