Runge-Kutta 四阶方法。向后整合答案

【问题标题】：Runge-Kutta fourth order method. Integrating backwardsRunge-Kutta 四阶方法。向后整合
【发布时间】：2019-10-14 13:12:07
【问题描述】：

我正在使用 Runge-Kutta 四阶方法来数值求解具有四次势的弯曲时空中背景标量场的通常运动方程：

$\phi^{''}=-3\left(1+\frac{H^{'}}{3H}\right)\phi^{'}-\lambda\phi^3/H^2$,

$'$ denoting the derivative w.r.t. the e-folds number $\textrm{d}N=H\textrm{d}t$ and, from the Friedmann equation:

$H^2=\frac{\lambda \phi^4}{4}\frac{1}{3M_{Pl}^2-(1/2)\phi^{'2}}$;

$H^{'}=-\frac{1}{2M_{Pl}^2}H\phi^{'2}$.

当使用前向积分后得到的最终值作为初始条件进行反向积分时，问题就出现了。向前积分时，结果与之前获得的值不匹配。我根本不明白问题出在哪里，因为方程和代码都不是未知的。首先，我整合了 0 到 64 个电子折叠。然后我简单地反转整合方向。

我也附上代码：

def rk4trial(f,v0,t0,tf,n,V):  
    t=np.linspace(t0,tf,n)
    h=t[1]-t[0]
    v=np.array((n+1)*[v0])
    for j in range(n):  
        V.append(v[j])
        V[j]=copy.deepcopy(V[j])
        k1=f(v[j],t[j])*h
        k2=f(v[j]+(1/2)*k1,t[j]+(1/2)*h)*h
        k3=f(v[j]+(1/2)*k2,t[j]+(1/2)*h)*h
        k4=f(v[j]+k3,t[j]+h)*h
        v[j+1]=v[j]+(k1+2*k2+2*k3+k4)/6
    return V, t, h


def Fdet(v,t):
    phi, sigma = v
    H=(((lamb/4)*phi**4)/(3*mpl**2-(1/2)*sigma**2))**(1/2)
    HH=-((1/2)*sigma**2)*(1/mpl**2)
    return np.array([sigma,-3*(1+HH/3)*sigma-lamb*phi**3/(H**2)])

PS：这里也发过这个问题：https://scicomp.stackexchange.com/questions/33583/runge-kutta-fourth-order-method-integrating-backwards，里面有详细的公式。

编辑：删除了代码中不必要的部分。

编辑：作为对@LutzL 的回应，我附上了 \phi/M_{Pl} 和 \phi^{'} 在向前（实线）和向后（虚线）积分后的图，通过做什么（s ）他说。如您所见，我无法解释前向积分结果的突然偏差。

【问题讨论】：

我想，这也不是问的正确地方。尝试将其移至scicomp.stackexchange.com
您能否更深入地解释数组V 和v 的作用？为什么V是函数调用中的参数？ V[j]=copy.deepcopy(V[j])这行背后的意图是什么？
@panch ：这并不是真正的算法设计问题。数值算法的所有部分都存在且正确。值得怀疑的是函数接口的设计，输入输出的假设，这是一个编程问题。
@LutzL：感谢您的评论。那是我编写的更复杂代码的一部分，我必须在其中执行多个平均值并保存结果。您实际上可以忽略它，它不会在此处修改结果。
scicomp.stackexchange.com/q/33583/11850 上的 cmets 听起来很明智，您不只是在尝试做一些对该方程组无效的事情

标签： python physics runge-kutta

【解决方案1】：

我会将 RK4 方法更改为必要的最低限度。没有必要让v 数组部分复制V 数组的内容，所以

def rk4trial(f,v0,t0,tf,n,V):  
    t=np.linspace(t0,tf,n)
    h=t[1]-t[0]
    v=v0
    for j in range(n):  
        V.append(v)
        k1=f(v,t[j])*h
        k2=f(v+0.5*k1,t[j]+0.5*h)*h
        k3=f(v+0.5*k2,t[j]+0.5*h)*h
        k4=f(v+k3,t[j]+h)*h
        v=v+(k1+2*k2+2*k3+k4)/6
    return V, t, h

这里不存在复制问题，因为v在每一步都是重新构造的，因此附加到返回数组的对象都是独立的。

后向集成应该和前向集成一样简单，

V1, t1, h1 = rk4trial(Fdet,v0,t0,tf,n,[])
V2, t2, h2 = rk4trial(Fdet,V1[-1],tf,t0,n,[])

和V2[k] 应该在方法错误的范围内与V1[-k-1] 相同。只有在刚性 ODE 中才会出现较大差异。

【讨论】：