Why does it matter at what point I replace?

Question

I need to find four orthogonal linear combinations of complicated functions, that vanish at four different points.

I use (LK4 is defined below, but its shape should have nothing to do with my problem)

NMinimize[LK4[{a, b, c, d}, I, 0, 0, 0, 0], {a, b, c, d}]

and get a result, as replacement rule (using MachinePrecision)

{3.3493*10^-17, {a -> 2.34594, b -> -1.80385, c -> 2.51873, d -> 0.2406}}

Now if I create a vector for these coefficients

P1 = {a, b, c, d}/.%[[2]]

and call the function with this vector I get a different result than by replacing directly in the function call:

LK4[{a, b, c, d}, I, 0, 0, 0, 0]/.%%[[2]] returns 3.3493*10^-17, while LK4[P1, I, 0, 0, 0, 0] returns 0.373542.

What's worse is that LK4[({a, b, c, d}/.%%%%[[2]]), 0, 0, 0, 0] also returns the inaccurate result 0.373542.

What I wanted to do is find a linear combination that minimizes at the first spot, then define a coefficient vector from that and minimize at the second spot with the constraint that the coefficient vector has to be orthogonal to the first. LK4[{a, b, c, d}, tau, xi1, xi2, x, y] normalizes the coefficient vector btw.

How can I get Mathematica to assign appropriately precise values upon using the replacement rule?

I tried setting the precision to 100 instead of MachinePrecision, but with the very same results.

Definition of LK4:

    Phi[j_, qM_, x_, y_, tau_, xi1_, xi2_] := 
  Surd[2 qM Im[tau]/Abs[tau]^2, 4] Exp[
    I 2 Pi (xi1 x + xi2 y)]  Exp[- Pi I Conjugate[tau]/Abs[tau]^2 ( 
       j^2/qM + qM x^2) + 
     2 Pi I j Conjugate[tau]/Abs[tau]^2 (x + tau y)] * 
   N[EllipticTheta[3, 
     Pi Conjugate[tau]/Abs[tau]^2 (qM (x + tau y) - j), 
     Exp[ -Pi I Conjugate[tau]/Abs[tau]^2 qM]], 20];
PhiEven[j_, qM_, x_, y_, tau_, xi1_, xi2_] := 
  Phi[j, qM, x, y, tau, xi1, xi2] /; j == 0;
PhiEven[j_, qM_, x_, y_, tau_, xi1_, xi2_] := 
  Phi[j, qM, x, y, tau, xi1, xi2] /;  
   qM/2 \[Element] Integers && j == qM/2;
PhiEven[j_, qM_, x_, y_, tau_, xi1_, xi2_] := 
  1/Sqrt[2] (Phi[j, qM, x, y, tau, xi1, xi2] + 
     Phi[qM - j, qM, x, y, tau, xi1, xi2]);
PhiOdd[j_, qM_, x_, y_, tau_, xi1_, xi2_] := 
  1/Sqrt[2] (Phi[j, qM, x, y, tau, xi1, xi2] - 
      Phi[qM - j, qM, x, y, tau, xi1, xi2]) /; j != 0 && j != qM/2;
LK4[coeff_, tau_, xi1_, xi2_, x_, y_] := 
  Sum[coeff[[a]]/(Sqrt@Total[coeff^2]) PhiEven[a - 1, 6, x, y, tau, 
     xi1, xi2], {a, 1, 4}];

Answer 1

Look at

LK4[{a, b, c, d}, I, 0, 0, 0, 0]

What has happened is that the a in the argument {a, b, c, d} has been replaced by {1, 2, 3, 4} in the Sum[..., {a, 1, 4}] code in the definition of LK4.

If you change the definition of LK4 to use a different iterator, you get consistent results:

LK4[coeff_, tau_, xi1_, xi2_, x_, y_] := 
 Sum[coeff[[a1]]/(Sqrt@Total[coeff^2]) PhiEven[a1 - 1, 6, x, y, tau, xi1, xi2],
 {a1, 1, 4}]

This happens because, according to the docs, Sum effectively uses Block to "localize" the iterator variable a. This means that any instances of a in the expression of the (evaluated) summand, such as an a in coeff argument {a, b, c, d}, will be replaced by 1, 2, etc.