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
Gradient 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)
