Ask Your Question
0

Problem with minimization and evaluation

asked 2022-04-19 23:17:23 +0200

Steve S gravatar image

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)

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2022-04-20 16:35:42 +0200

Max Alekseyev gravatar image

updated 2022-04-20 16:36:32 +0200

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, ...})
edit flag offensive delete link more

Comments

Thanks, it works fine !

Steve S gravatar imageSteve S ( 2022-04-20 20:00:11 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2022-04-18 16:19:19 +0200

Seen: 41 times

Last updated: Apr 20