# Finding a monotonic polynomial over an interval

I want to find a polynomial of specified degree $$d$$ defining the function $$f:[0,1]\to[0,1]$$ satisfying $$f(0)=0,f(1)=1$$ and $$f$$ is monotonically increasing over the interval. This seems like an easy enough task especially since there's the trivial solution $$f(x)=x^d$$, and although my code seems to work decently well for $$d=2,3$$, it takes around a minute on my machine to produce an answer for $$d=4$$ and it couldn't get anything in the time I set it running for $$d=5$$. I'm looking to test this for values of $$d$$ up to $$10^3$$ and find multiple instances for $$f$$, so this is certainly not going to work. What am I doing wrong and how can I make this efficient?

Here is my code:

deg = 4;
f[x_] := Sum[Subscript[c, k] x^k, {k, 1, deg}];
FindInstance[{
f[1] == 1,
ForAll[x, 0 <= x <= 1, f'[x] >= 0]
}, Table[Subscript[c, i], {i, 1, deg}], Reals]

• The space of these monotonic functions is quite large. Are there other considerations one might use to narrow the search space? – Daniel Lichtblau Jan 2 at 2:42
• @DanielLichtblau Yes, for my applications I'm actually interested in such functions satisfying $\int_0^1xf(x)dx=K$ where $K$ is a fixed constant. (Note that not every value of $K$ works, a safe choice would be somewhere around $1/3$.) – YiFan Jan 2 at 2:45
• If you are willing to restrict the polynomial to having positive coefficients, that will guarantee monotonicity. Then FindInstance, for example, can be used: In[50]:= d = 5; coeffs = Array[a, d + 1, 0]; integral = Integrate[x*coeffs.x^Range[0, d], {x, 0, 1}]; k = 1/3; FindInstance[Flatten@{integral == 1/3, Thread[0 <coeffs]}, coeffs] Out[54]= {{a[0] -> 1/48, a[1] -> 1/2, a[2] -> 1/3, a[3] -> 5/24, a[4] -> 1/8, a[5] -> 7/96}} – Daniel Lichtblau Jan 2 at 15:16

Try this

deg=4;
c=Join[{0.},Sort[RandomReal[{0,1},deg-1]],{1.}];
d=Rest[c]-Most[c];
v=Table[x^n,{n,1,deg}];
poly=Dot[d,v]
Plot[poly,{x,0,1},PlotRange->All]


That finds a list of positive Real coefficients which sum to 1 and builds your polynomial from them. Thus it monotonically increases from 0 to 1.

Numerical accuracy may become an issue with degree 1000 so you may need to modify this to generate a list of positive exact rational coefficients which sum to 1 to avoid that problem. But the current form completes the calculations in a few seconds with deg=1000.

• Thank you! Do you think it's possible to modify this idea to take into account a constraint of the form $\int_0^1xf(x)dx=K$ for a constant $K$? (Equivalently we would need $\sum a_n/(n+2)=K$ where $a_n$ are the coefficients, but I don't know how to take into account both this and $\sum a_n=1$.) – YiFan Jan 2 at 3:39
• I do not see a way to find exactly K. If an approximation would do then deg=4;k=1/3;besterr=Infinity; Do[c=Join[{0.},Sort[RandomReal[{0,1},deg-1]],{1.}];d=Rest[c]-Most[c];v=Table[x^n,{n,1,deg}];poly=Dot[d,v];int=NIntegrate[x*poly,{x,0,1}];If[besterr>Abs[int-k], besterr=Abs[int-k];Print[besterr];bestpoly=poly],{100}]; bestpoly can try to find a polynomial close to your condition. Change the number of iterations if needed. I do not see how using any Mathematica minimization function would do better than this. – Bill Jan 2 at 4:17