2
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ The space of these monotonic functions is quite large. Are there other considerations one might use to narrow the search space? $\endgroup$ – Daniel Lichtblau Jan 2 at 2:42
  • $\begingroup$ @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$.) $\endgroup$ – YiFan Jan 2 at 2:45
  • $\begingroup$ 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}} $\endgroup$ – Daniel Lichtblau Jan 2 at 15:16
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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$.) $\endgroup$ – YiFan Jan 2 at 3:39
  • $\begingroup$ 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. $\endgroup$ – Bill Jan 2 at 4:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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