如何将多变量函数集成到单个变量中，并最小化返回的单变量函数？ （蟒蛇）[关闭]

Question

我是编码和Python新手，最近被分配了一个项目来求解氢原子的变分波函数。我不需要对其进行归一化，只需求解波函数本身内的常数因子。 经过一些代数之后，我能够用函数简化哈密顿量定义的能量方程

E(k) = Ak^2 + BI(k)

对于最初猜测的波函数

ψ(r,k) = e^(-kr)

其中 E(k) 是哈密顿量定义的能量，A 和 B 是常数，k 是初始猜测的波函数内的常数因子，I(k) 是积分分数

I(k) = 积分((e^(-2kr)/r)dr) / 积分(e^-2kr)dr

（抱歉，我还没有自学 LaTeX！）

我的这段代码的目标是最小化 E(k) 来找到满足 ψ(r,k) 的 k 值。问题是我找不到任何方法返回 I(k) ，使其仍然是 k 的函数而不是常数，而且我对 python 缺乏经验，不知道如何真正发挥它的创意。

这是我当前的（完全不起作用）代码，供参考：

#import packages
import numpy as np
import scipy
import scipy.constants as constants
from scipy.integrate import quad 
import scipy.optimize
from scipy.optimize import minimize

#defining constants
A = (-43.9048056327210225) / (constants.pi**2 * 800 * 9.1093837)
B = (2.5669699537181569) / (8.854187817 * 4 * constants.pi)
k = 1 #this is my problem; quad won't run without this

def I_N(r): #numerator integrand
    return (np.exp(-2*k*r) * r**(-1))
def I_D(r): #denominator integrand
    return (np.exp(-2*k*r))
numerator, _ = quad(I_N, 0, 100)
denominator, _ = quad(I_D, 0, 100)
I = numerator/denominator

def E(k):
    return (A*k**2) + B*I

H_Atom = minimize(E, k, tol=1e-20)
optimal_k = H_Atom.x
print("Optimal k: " , optimal_k)

我考虑将积分 I(k) 重新定义为总和，但我不知道如何将 sigma 表示法转换为 python，因为积分((e^(-2kr)/r)dr) 是一个发散级数.

请记住，我不允许使用线性代数来求解，只能使用变分法。另外，请记住，我非常缺乏经验（这是我的第一个编码项目），所以请保持语言平易近人，以便我能够理解！谢谢你:)

注意，k 没有特别的限制，只要它返回初始波函数 e^(-kr) 的最小值即可。数学计算应该是正确的（为了计算机计算的准确性，已经删除了 10^-17 倍）。虽然我不会写出常数 A 和 B 的全部计算，但一般形式应该是正确的，因为它遵循 E = KE + PE 的一般格式，其中 KE 的平方是 a 的结果微分算子（哈密顿量被设计为返回遵循这种形式的微分方程，因为它是一个能量算子）。

Answer 1

如果我们重新表述您的问题，仅用

k

的函数来表达它，然后采取一些数字预防措施，就有可能解决它：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize 

A = (-43.9048056327210225) / (np.pi**2 * 800 * 9.1093837)
B = (2.5669699537181569) / (8.854187817 * 4 * np.pi)

def numerator(k):
    def integrand(r):
        return np.exp(-2 * k * r) / r
    return integrate.quad(integrand, 1e-8, np.abs(10. * k))[0]

def denominator(k):
    def integrand(r):
        return np.exp(-2 * k * r)
    return integrate.quad(integrand, 1e-8, np.abs(10. * k))[0]

def I(k):
    return numerator(k) / denominator(k)

def E(k):
    return A * k**2 + B * I(k)

solution = optimize.minimize(E, [1.], bounds=[(0, np.inf)])

#      fun: 0.29720690448818193
# hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
#      jac: array([-4.05231404e-07])
#  message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
#     nfev: 10
#      nit: 3
#     njev: 5
#   status: 0
#  success: True
#        x: array([0.25478549])

然后我们可以确认我们已经找到了正常数的最小值：

k = np.logspace(-2, 1, 200, base=10)
Ev = np.vectorize(E)

fig, axe = plt.subplots()
axe.loglog(k, Ev(k))
axe.scatter(solution.x, Ev(solution.x))
axe.grid()