2 added 403 characters in body edited Jan 2 '16 at 19:12 Stelios 1,33111 gold badge99 silver badges1313 bronze badges One approach to this kind of problems is to discretize the continuous domain to a finite (but sufficiently large) number of (not necessarily equispaced) samples and then approximate the original problem on this grid. The resulting formulation is a linear programming problem for which efficient solvers are available. f[x_] := Exp[x]^(1/2) + 2 - Exp[x] + x^5; l = 5000; (*number of samples of the [0,1] interval*) xRange = Range[0., 1., 1./(l - 1)]; (*grid points*) (*auxiliary matrix; note that {a,b,c}.mA is ax^2+bx+c evaluated on the grid points*) mA = {xRange^2, xRange, Table[1., {l}]}; (*auxiliary matrix; note that {a,b,c}.mB is the derivative of ax^2+bx+c evaluated on the grid points*) mB = {2 xRange, Table[1., {l}]}; sol = (* note introduction of an additional optimization variable delta *) FindMinimum[{delta, And @@ Thread[-delta <= {a, b, c}.mA - f[xRange] <= delta] && And @@ Thread[{a, b}.mB <= 0]}, {delta, a, b, c}]  {0.181843, {delta -> 0.181843, a -> 0.0695606, b -> -0.139121, c -> 1.81816}}  Plot[{f[x], Evaluate[a x^2 + b x + c /. sol[[2]]]}, {x, 0, 1}]  The optimal polynomial approximation in this case is pretty far from close to the actual function due to the non-increasing constraint. Actually, generalizing this approach to larger order polynomials shows that the approximation gets only slightly better. Of course, removing this constraint results in much better approximations that rapidly converge to the actual function with increasing order. One approach to this kind of problems is to discretize the continuous domain to a finite (but sufficiently large) number of (not necessarily equispaced) samples and then approximate the original problem on this grid. The resulting formulation is a linear programming problem for which efficient solvers are available. f[x_] := Exp[x]^(1/2) + 2 - Exp[x] + x^5; l = 5000; (*number of samples of the [0,1] interval*) xRange = Range[0., 1., 1./(l - 1)]; (*grid points*) (*auxiliary matrix; note that {a,b,c}.mA is ax^2+bx+c evaluated on the grid points*) mA = {xRange^2, xRange, Table[1., {l}]}; (*auxiliary matrix; note that {a,b,c}.mB is the derivative of ax^2+bx+c evaluated on the grid points*) mB = {2 xRange, Table[1., {l}]}; sol = (* note introduction of an additional optimization variable delta *) FindMinimum[{delta, And @@ Thread[-delta <= {a, b, c}.mA - f[xRange] <= delta] && And @@ Thread[{a, b}.mB <= 0]}, {delta, a, b, c}]  {0.181843, {delta -> 0.181843, a -> 0.0695606, b -> -0.139121, c -> 1.81816}}  Plot[{f[x], Evaluate[a x^2 + b x + c /. sol[[2]]]}, {x, 0, 1}]  One approach to this kind of problems is to discretize the continuous domain to a finite (but sufficiently large) number of (not necessarily equispaced) samples and then approximate the original problem on this grid. The resulting formulation is a linear programming problem for which efficient solvers are available. f[x_] := Exp[x]^(1/2) + 2 - Exp[x] + x^5; l = 5000; (*number of samples of the [0,1] interval*) xRange = Range[0., 1., 1./(l - 1)]; (*grid points*) (*auxiliary matrix; note that {a,b,c}.mA is ax^2+bx+c evaluated on the grid points*) mA = {xRange^2, xRange, Table[1., {l}]}; (*auxiliary matrix; note that {a,b,c}.mB is the derivative of ax^2+bx+c evaluated on the grid points*) mB = {2 xRange, Table[1., {l}]}; sol = (* note introduction of an additional optimization variable delta *) FindMinimum[{delta, And @@ Thread[-delta <= {a, b, c}.mA - f[xRange] <= delta] && And @@ Thread[{a, b}.mB <= 0]}, {delta, a, b, c}]  {0.181843, {delta -> 0.181843, a -> 0.0695606, b -> -0.139121, c -> 1.81816}}  Plot[{f[x], Evaluate[a x^2 + b x + c /. sol[[2]]]}, {x, 0, 1}]  The optimal polynomial approximation in this case is pretty far from close to the actual function due to the non-increasing constraint. Actually, generalizing this approach to larger order polynomials shows that the approximation gets only slightly better. Of course, removing this constraint results in much better approximations that rapidly converge to the actual function with increasing order. 1 answered Jan 2 '16 at 18:47 Stelios 1,33111 gold badge99 silver badges1313 bronze badges One approach to this kind of problems is to discretize the continuous domain to a finite (but sufficiently large) number of (not necessarily equispaced) samples and then approximate the original problem on this grid. The resulting formulation is a linear programming problem for which efficient solvers are available. f[x_] := Exp[x]^(1/2) + 2 - Exp[x] + x^5; l = 5000; (*number of samples of the [0,1] interval*) xRange = Range[0., 1., 1./(l - 1)]; (*grid points*) (*auxiliary matrix; note that {a,b,c}.mA is ax^2+bx+c evaluated on the grid points*) mA = {xRange^2, xRange, Table[1., {l}]}; (*auxiliary matrix; note that {a,b,c}.mB is the derivative of ax^2+bx+c evaluated on the grid points*) mB = {2 xRange, Table[1., {l}]}; sol = (* note introduction of an additional optimization variable delta *) FindMinimum[{delta, And @@ Thread[-delta <= {a, b, c}.mA - f[xRange] <= delta] && And @@ Thread[{a, b}.mB <= 0]}, {delta, a, b, c}]  {0.181843, {delta -> 0.181843, a -> 0.0695606, b -> -0.139121, c -> 1.81816}}  Plot[{f[x], Evaluate[a x^2 + b x + c /. sol[[2]]]}, {x, 0, 1}]