Following the approach as given in the answer of Belisarius here I just wonder how one could incorporate some more sophisticated requirements on the approximating function.

Suppose the approximating model is as follows:

 funcr[al_?NumericQ, be_?NumericQ, ga_?NumericQ, de_?NumericQ, 
       k1_?NumericQ, k2_?NumericQ, x_?NumericQ] := 
   (x^2 al + x be + k1) (x^2 ga + x de + k2)

And there are two extra constraints on the approximating function:

  • part (x^2 ga + x de + k2) should integrate to 1 (on $x\in[0,1]$) and
  • part (x^2 al + x be + k1) should be non-decreasing (on [0,1]).

I can easily incorporate the first requirement: just integrate this function and get the direct constraint: de/2 + ga/3 + k2 == 1 that I can introduce into NMinimize[..], so the final code is something like this:

piece[x_] := Piecewise[{{x^3, 1 >= x >= 0}, {1, x > 1}}, 0]; *"data"*

funcr[al_?NumericQ, be_?NumericQ, ga_?NumericQ, de_?NumericQ, 
      k1_?NumericQ, k2_?NumericQ, x_?NumericQ] := 
  (x^2 al + x be + k1) (x^2 ga + x de + k2);

norm[al_, be_, ga_, de_, k1_, k2_,x_] := 
  (funcr[al, be, ga, de, k1, k2, x] - piece[x])^2

int[al_?NumericQ, be_?NumericQ, ga_?NumericQ, de_?NumericQ, 
    k1_?NumericQ, k2_?NumericQ] := 
  NIntegrate[norm[al, be, ga, de, k1, k2, x], {x, 0, 1}]

  nm = 
      {int[al, be, ga, de, k1, k2], de/2 + ga/3 + k2 == 1}, 
      {al, be, ga, de, k1, k2}, 
      Method -> "NelderMead", AccuracyGoal -> 3]

But how can I introduce the condition on the first part, the requirement that (x^2 al + x be + k1) is non-decreasing on [0,1]? This involves $x$ now (namely derivative of (x^2 al + x be + k1) is positive) while NMinimize has no idea that $x$ exists.

I think this should be a part of the definition of the model, i.e. funcr but my attempts write it as a list of conditions all fail (neither makes sense to me.. nor works once in the program).

  • $\begingroup$ In this case that derivative is just 2*a1*x+be and insisting it be nonnegative on x in [0,1] reduces to checks at the endpoints since it is linear in x. So you have as constraints be>=0 and 2*a1+be>=0. $\endgroup$ – Daniel Lichtblau Nov 10 '15 at 17:01
  • $\begingroup$ Thank you! Yes, this works surely. But I would like to know how to incorporate such constraints generally, for more complicated models, when such tricks are unavailable... The above model just a toy example. $\endgroup$ – Kass Nov 10 '15 at 17:17
  • $\begingroup$ It can become quite difficult of course. One possibility, not foolproof, is to break the interval into smallish parts and just check endpoints of those subintervals for nonnegativity of the first derivative (or check that the function itself is nondecreasing as you go from 0 to 1 at those particular points). If the function can evaluate Interval values then a sufficient condition is that the derivative thus eval'd be nonnegative (but this is not a necessary condition, and so could cause one to lose possible optimizing values). Requires some finesse, I think. $\endgroup$ – Daniel Lichtblau Nov 10 '15 at 17:25
  • $\begingroup$ Thanks yes it would work as well.. But for a multidimensional model, with higher degrees of polynomials etc that would be a nightmare. Moreover I am sure there should be some simple and efficient way to add those extra constraints, given which type of other complex problems could be solved almost automatized in Mathematica.. $\endgroup$ – Kass Nov 10 '15 at 18:25
  • $\begingroup$ Hi Daniel again, in fact after some thinking I realize that the method you suggest is pretty good one. As I am trying to approximate by polynomials (and so rather smooth functions), imposing constraints on a grid of values should be rather sufficient. I am currently trying to implement it with 10 values on $[0,1]$, which must pin down the polynomial pretty well. Of course the question of general approach remain valid, not least as the matter of curiosity! $\endgroup$ – Kass Nov 11 '15 at 0:43

Your Answer

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

Browse other questions tagged or ask your own question.