# Problem with minimization and evaluation

Hi, sorry for my English, i'm not a native speaker and a very beginner in Sage.
I have a problem with the minimize function, my goal is to minimize the Hamiltonian of a certain function psi, to find the energy of a helium atom.
So I have a function psi depending of three variables (s t u) and 19 coefficients alpha, c1 to c18 and I build the Hamiltonian Hmoy with this function (Hmoy depends only of the alpha, c1 ... c18 and Z, which is fixed to 2) and I try to minimize this.

The minimization function return me 19 numerical values for the coefficients and a evaluated minimized Hmoy value of -2.903293 which is the expected value, but the problem happen when I evaluate manually Hmoy with the given coefficients, the numerical evaluation is systematically different from what the minimize function return. This is a problem because I need to reuse the coefficients from another calculus.
Do you have an idea of what's going on here ?
Here's my code :

 R = RealField(100)
RealNumber = R """Variables declaration"""
s,t,u = var('s t u')
alpha,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18 = var('alpha c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18')
Z = var('Z') """Wave-function definition (minimizing c1 ... c18)"""
psi(s,t,u) = exp(-alphas) * (c1 + c2su + c3ut + c4st + c5stu + c6s^2tu + c7st^2u + c8stu^2 + c9s^2t^2u^2 + c10s + c11t + c12u + c13s^2t + c14s^2u + c15t^2s + c16t^2u + c17u^2s + c18u^2*t) """Assumptions to evaluate integrals"""
assume(alpha>0)
assume(beta>0)
assume(c1>0)
assume(c2>0)
assume(Z>1)
assume(s>0)
assume(u>0)
assume(u-s>0) """Integrals"""
norme = integrate(integrate(integrate(pi^2(s^2-t^2)u*psi(s,t,u)^2 ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo)) I1 = integrate(integrate(integrate(pi^2(s^2-t^2)upsi(s,t,u)^2(-2Z/(s+t)) ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo))
I2 = integrate(integrate(integrate(pi^2(s^2-t^2)upsi(s,t,u)^2(-2Z/(s-t)) ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo)) I3 = integrate(integrate(integrate(pi^2(s^2-t^2)u * psi(s,t,u)^2 * 1/u ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo)) I4 = integrate(integrate(integrate(-1/2 * pi^2(s^2-t^2)u psi(s,t,u) * ( diff(psi(s,t,u),s,2) + diff(psi(s,t,u),t,2) + diff(psi(s,t,u),u,2)+ 4/(s+t) * ( diff(psi(s,t,u),s) + diff(psi(s,t,u),t)) + 2/u * diff(psi(s,t,u),u) + 2 * diff(diff(psi(s,t,u),s),t) + 2 * (st+u^2)/(u(s+t)) * ( diff(diff(psi(s,t,u),s),u) + diff(diff(psi(s,t,u),t),u) ) ) ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo))
I5 = integrate(integrate(integrate(-1/2 * pi^2(s^2-t^2)u *psi(s,t,u) * ( diff(psi(s,t,u),s,2) + diff(psi(s,t,u),t,2) + diff(psi(s,t,u),u,2)+ 4/(s-t) * ( diff(psi(s,t,u),s) - diff(psi(s,t,u),t)) + 2/u * diff(psi(s,t,u),u) - 2 * diff(diff(psi(s,t,u),s),t) + 2 * (st-u^2)/(u*(s-t)) * ( diff(diff(psi(s,t,u),t),u) - diff(diff(psi(s,t,u),s),u) ) ) ,(t,-u,u)) ,(u,0,s)) ,(s,0,+oo)) """Function to minimize (c1...c18 the minimizing coefficients, Z fixed to 2)"""
Hmoy(alpha,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,Z) = (I1 + I2 + I3 + I4 + I5) /norme forget(alpha>0) """Minimization"""
S = minimize(Hmoy(alpha,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,2), [2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],algorithm="bfgs", verbose=True) """printing results"""
print("Manually evaluated Hmoy value :")
print(Hmoy(2.091464938388313, 9.267331957693454, 2.7093098187035247, -0.004815455491397233, 3.789378019957407, -0.0017404538769011718, 0.9076364671216567, 1.5328777169654666, 0.2982270361252481, -0.4727098102678102, -0.0019807632284596236, -0.17237817441579442, 0.0038481273111783993, 0.004918946275046137, -0.0009974484597963381, 0.00041078260111615007, -0.44701927735295105, 0.0003337602854370373, 0.014689271189557896,2))
print("List of minimized coefficients alpha, c1 ... c18 :")
print(S)


And the results :


Optimization terminated successfully.
Current function value: -2.903293
Iterations: 147
Function evaluations: 162
Manually evaluated Hmoy value :
-2.0749434605617761790619679333
List of minimized coefficients alpha, c1 ... c18 :
(2.091464938388313, 9.267331957693454, 2.7093098187035247, -0.004815455491397233, 3.789378019957407, -0.0017404538769011718, 0.9076364671216567, 1.5328777169654666, 0.2982270361252481, -0.4727098102678102, -0.0019807632284596236, -0.17237817441579442, 0.0038481273111783993, 0.004918946275046137, -0.0009974484597963381, 0.00041078260111615007, -0.44701927735295105, 0.0003337602854370373, 0.014689271189557896)


edit retag close merge delete

Sort by » oldest newest most voted

The issue here is that the meaning of c1, c2, ... is intermixed between parameters and function arguments. Take a look - you define:

Hmoy(alpha,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,Z) = (I1 + I2 + I3 + I4 + I5) /norme


but the right-hand side here does not take into account the given values of alpha,c1,....

So, either psi, I1, ..., I5, norme need to be defined as functions taking arguments alpha,c1,... explicitly, OR keep them as symbolic expressions and define Hmoy as such:

Hmoy = (I1 + I2 + I3 + I4 + I5) /norme


and evaluate it via substitution of parameters:

Hmoy.subs({alpha:2.091464938388313, c1:9.267331957693454, ...})

more

Thanks, it works fine !

( 2022-04-20 20:00:11 +0200 )edit