计算科学 - 使用具有自适应步长的 Runge-Kutta 将径向薛定谔方程与 Lennard-Jones 势积分最终步长为零 - 吾爱随笔录

使用具有自适应步长的 Runge-Kutta 将径向薛定谔方程与 Lennard-Jones 势积分最终步长为零

计算科学量子力学龙格库塔

2021-12-23 18:40:09

我目前正在学习计算物理课程。我是计算物理学和一般编程的新手。我正在使用数值方法来尝试将径向薛定谔方程与 Lennard-Jones 势相结合。

[\frac{ℏ^{2}}{2 m} \frac{d^{2}}{d r^{2}} + (E - V (r) - \frac{ℏ^{2} l (l + 1)}{2 m r^{2}})] u_{l} (r) = 0

$\left[ \frac{\hbar^2}{2m}\frac{d^2}{dr^2} + \left( E-V(r)-\frac{\hbar^2 l (l+1)}{2mr^2}\right) \right] u_l(r)=0$

V (r) = ϵ [{(\frac{ρ}{r})}^{12} - 2 {(\frac{ρ}{r})}^{6}]

$V(r)= \epsilon \left[ \left( \frac{\rho}{r} \right)^{12}-2\left(\frac{\rho}{r}\right)^6 \right]$

Numerical recipes 有一个名为 odeint 的函数，它将使用五阶 Runge-Kutta 算法为您积分一个常微分方程。该函数具有可调整的步长，这似乎导致我的代码出现问题。也就是说，我的步长将变为零，这会导致数值配方抛出错误并过早退出。

我在数值上从到进行积分，并在分析上达到我的最小值，以处理零处的奇点。我在下面附上了我的代码，有关 odeint 及其工作原理的更多信息，请访问：http ://www.itp.uni-hannover.de/Lehre/cip/odeint_c.pdf 。谁能帮我理解我哪里出错了？ $r_{min}=\frac{\rho}{2}$ $r_{max}=5\rho$

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

#define N 2

float dxsav,*xp,**yp;  /* defining declarations */
int kmax,kount;

int nrhs;   /* counts function evaluations */

/* Schrodinger equation and L-J parameters */
double alpha = 6.12; 
double rho = 3.57;
double epsilon = 5.9;
double l = 1.0;
double energy = 3;

double rmin=1.785;
double rmax=17.85;


void derivs(float x,float y[],float dydx[])
{
    nrhs++;
    printf("xodeint check: x=%f\n", x);
    dydx[1] = y[2];
    dydx[2]=alpha*((l*(l+1)/(alpha*x*x))+epsilon*(pow(rho/x,12)-2*pow(rho/x,6))-energy)*y[1];
}

int main(void)
{
    int i,nbad,nok;
    float eps=1.0e-4,h1=0.1,hmin=0,x1=rmin,x2=rmax,*ystart;

    ystart=vector(1,N);
    xp=vector(1,200);
    yp=matrix(1,10,1,200);
    ystart[1]=0.93583;
    ystart[2]=0.17385;
    nrhs=0;
    kmax=100;
    dxsav=(x2-x1)/20.0;
//  printf("%f\n", h1);
    odeint(ystart,N,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqs);
    printf("\n%s %13s %3d\n","successful steps:"," ",nok);
    printf("%s %20s %3d\n","bad steps:"," ",nbad);
    printf("%s %9s %3d\n","function evaluations:"," ",nrhs);
    printf("\n%s %3d\n","stored intermediate values:    ",kount);
    printf("\n%8s %18s %15s\n","r","integral","x^2");
    for (i=1;i<=kount;i++)
        printf("%10.4f %16.6f %14.6f\n",xp[i],yp[1][i],xp[i]*xp[i]);
    free_matrix(yp,1,10,1,200);
    free_vector(xp,1,200);
    free_vector(ystart,1,N);
    return 0;
}
#undef NRANSI

这输出

Numerical Recipes run-time error...
stepsize underflow in rkqs
...now exiting to system...

2个回答

我从我自己的研究中很清楚这个问题：这是因为方程非常僵硬。因此，您很可能没有做错任何事情（尽管我没有检查您的代码）。那么僵硬从何而来？

问题

为了对方程进行数值求解，通常选择一个网格并使用有限差分来离散化函数和算子。为简单起见，假设您的网格是间距为的空间。通常，为了获得良好的近似，人们想选择 $\Delta x$ $\Delta x$

现在对方程中的项进行分析：动能项（即导数）导致的上能量，作为可表示的最大频率网格上是。（我很懒所以我在这里写你肯定知道公式）。该项包含在任何薛定谔方程中，并不真正构成问题。 $\mathcal O(\frac{1}{\Delta x^2})$ $\mathcal O(\frac{1}{\Delta x})$ $\mathcal O$

现在让我们考虑潜在的功能。这些是网格上的对角线，因此函数将具有最大能量。在您的等式中，您的最大为。正是这个词主要扼杀了你的传播。 $1/r^k$ $1/(\Delta x^k)$ $k$ $12$

解释

为什么这个词会带来很大的问题？考虑您的完整离散哈密顿矩阵（对称），以及您的波函数在网格方面的扩展。薛定谔方程的显式解由下式给出 $H$

ψ (t) = \exp (i H t) ψ (t = 0) .

$\psi(t) = \exp(i H t) \psi(t=0)\,.$

请注意，没有时间排序算子，因为哈密顿量是与时间无关的。现在对您的哈密顿量进行特征分解。然后，指数可以表示为。如果您现在在对角矩阵，那么您的传播中有一个类似的项。 $H = U E U^\dagger$ $\exp(i H t) = U \exp(i E t) U^\dagger$ $e_{max}$ $E$ $exp(ie_{max} t)$

这里进入香农定理（或者傅里叶变换），它指出为了很好地采样这个函数，它应该是。否则，你会得到奇怪的东西，比如混叠效果，这会破坏你的解决方案。 $\Delta t \sim 1/e_{max}$

为方便起见，可以将完整哈密顿量的最大本征能量近似为其各项的最大本征能量之和。有了这个，您会看到时间步长需要像一样小。也就是说，如果你的是（这里没有单位，你知道它们），你需要一个的时间步长。自适应步长积分器将因此越来越小步长，直到它们达到其下溢常数。 $\Delta t \sim 1/\Delta x$ $\Delta x$ $10^{-2}$ $\sim 10^{-24}$

这是一个相当非数学的处理，所以不要引用我的数字，但它直观地解释了正在发生的影响。

建议的解决方案

我认为有一个简单的替代方法适合您的问题，通常称为光谱法。如果您的网格具有中小型尺寸 - 比如说最多个网格点 - 我猜这里就是这种情况，它是适用的。那么，事实上，它是您的问题的选择方法。 $N=10000$

解决方案就是简单地计算矩阵指数，如上图所示，并将其应用于初始波函数向量。对角化的成本是，但您只需要执行一次。另外，通过这种方式，您可以获得精确的数值解——所有这些 ODE 方法都试图实现的解。 $\exp(iHt)$ $\mathcal O(N^3)$

缺点是这种方法仅适用于哈密顿量与时间无关的情况。否则，当它随时间变化时，您需要经常重新对角化它。

我要尝试的第一件事是隐式积分器，以检查您的 ODE 是否僵硬。如果您使用 RK5（一种显式方法）来求解 ODE 并且步长变得非常小，则可能是稳定性考虑严重限制了您的步长。如果您想使用用 C 编写的 ODE 求解器，SUNDIALS 或 PETSc 是生产代码的不错选择。PETSc 得到很好的支持。SUNDIALS 有一个邮件列表，它也获得了相当多的流量。GNU 科学图书馆也有一些可能对学术使用或原型设计有益的例程。

您也可以尝试不同的 RK 实现。与从数字食谱转录代码的人相比，我更有可能信任人们使用和支持的库中的代码。关于数字食谱的意见各不相同。作者自己说，对库中错误的批评是没有根据的，并列出了当前的错误报告（现在已经过时，因为他们有一个论坛）。然而，我听过的关于数字食谱的大多数故事都是这样的——不是很好。你的旅费可能会改变。

最后，作为一个总括，您可以仔细检查（最好是测试）以确保在您的 ODE 右侧的实现中没有错误，但这一步可能是我排除第一步后的最后一步问题的两个来源。

其它你可能感兴趣的问题

上一篇Nesterov 的第一种分段线性凸优化问题方法中的 beta 下一篇C++中一般复矩阵的特征值