0

I'm trying to minimize a quadratic energy with linear equality constraints in Maple. Currently I'm issuing things like:

with(Optimization):
p := (t) -> c3*t^3 + c2*t^2;
m := Minimize(int(diff(p(t),t)^2,t=0..1),{eval(p(t),t=1)=1,eval(diff(p(t),t),t=1)=0});

but this seems to give me a numerically optimized solution complete with floating point error:

m := [1.19999999999997, [c2 = 3.00000000000000, c3 = -2.00000000000000]]

(The correct answer is m:= [6/5,[c2=3,c3=-2]])

Is there a way to compute the solution symbolically using maple?

I'd rather not have to work out the Lagrangian myself. I'm hoping for a flag like symbolic=true.

Mr.Wizard
  • 24,179
  • 5
  • 44
  • 125
Alec Jacobson
  • 6,032
  • 5
  • 51
  • 88

1 Answers1

1

Apply Lagrange multipliers step-by-step. It's fairly easy to do.

p:= t-> c3*t^3+c2*t^2:
Obj:= int(D(p)^2,  0..1):
Con:= {p(1)=1, D(p)(1)=0}:
L:= Obj - l1*lhs(Con[1]) - l2*lhs(Con[2]):
solve(Con union {diff(L,c3), diff(L,c2)});
              /                      12       -1\ 
             { c2 = 3, c3 = -2, l1 = --, l2 = -- }
              \                      5        5 / 
eval(Obj, %);
                               6
                               -
                               5
Carl Love
  • 1,461
  • 8
  • 7
  • So I also wrote out the Lagrangian http://www.alecjacobson.com/weblog/?p=3425 , but I noticed that with polynomials of degree ~15 or so the solutions coming out of solve are no good: they don't seem to satisfy the constraints anymore if I evaluate a polynomial with those coefficients. Is `solve` not robust to handle larger degree polys? (BTW is it easy to load the actually values of c2,c3,etc. found by solve into p so that I can evaluate: `p(1)`?) – Alec Jacobson Aug 21 '13 at 19:56
  • 1
    I followed your web link, but I just see the original problem, not one with degree 15. Can you post an example where solve doesn't seem to work? You can use the *assign* command to assign the values of c2 and c3. – Carl Love Aug 22 '13 at 15:17
  • 1
    I don't see how you can get a symbolic solution with any practical value from **solve** for a high-degree polynomial, but you shouldn't get something that doesn't satisfy the constraints; that would be a bug. – Carl Love Aug 22 '13 at 15:21
  • Here's the 15-degree for reference `p := (t) -> c02*t^2+c03*t^3+c04*t^4+c05*t^5+c06*t^6+c07*t^7+c08*t^8+c09*t^9+c10*t^10+c11*t^11+c12*t^12+c13*t^13+c14*t^14+c15*t^15:` After using the `assign` command (thank you), I see that `solve` is indeed satisfying the constraints. I was evaluating it numerically and plotting the function. This gave me garbage. Interestingly, the resulting coeffs from `Minimize` is pretty different. Though `Minimize` only satisfies the constraints up to about single precision. – Alec Jacobson Aug 22 '13 at 16:01
  • 1
    High-degree polynomial solutions are notoriously unstable numerically. But just because it doesn't appear to satisfy the polynomial when you plug it back in does not mean that the root was not accurate to the full number of digits specified. As an experiment to see what I mean, try `expand((x-10.01)^15); eval(%, x= 10.01);` – Carl Love Aug 22 '13 at 17:27
  • Agreed. Nice example. I should have held my breath until checking the symbolic solution. – Alec Jacobson Aug 22 '13 at 19:49