使用 scipy.solve_bvp 求解 BVP，其中函数返回一个数组答案

【问题标题】：Solving a BVP using scipy.solve_bvp where the function returns an array使用 scipy.solve_bvp 求解 BVP，其中函数返回一个数组
【发布时间】：2020-03-06 15:30:11
【问题描述】：

这是一个非常笼统的问题，因为我觉得我的错误是由于对 scipy.solve_bvp 工作原理的一些误解造成的。我有一个函数def，它接受一个包含 12 个数字的数组，并返回给定时间的微分方程组列表，形状为 (2,6)。我将有一个长度为 n 的一维数组作为我的时间步长，然后是一个数组 yof 形状为 (12,n) 的输入值。我的代码旨在模拟地球和火星在 1000 天内受边界条件影响的运动；在 t=0 位置 = rpast（对应的速度由函数 find_vel_past() 返回），在 t=1000 的位置和速度分别由 rs 和 vs 给出。我的代码位于底部，上面有我试图解决的两个函数：

from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from scipy import integrate
from scipy import signal

G       = 6.67408e-11 # m^3 s^-1 kg^-2
AU      = 149.597e9 # m
Mearth  = 5.9721986e24 # kg
Mmars   = 6.41693e23 # kg
Msun    = 1.988435e30 # kg
day2sec = 3600 * 24 # seconds in one day

rs = [[-4.8957151e10, -1.4359284e11, 501896.65],  # Earth
      [-1.1742901e11, 2.1375285e11, 7.3558899e9]] # Mars (units of m)
vs = [[27712., -9730., -0.64148], # Earth
      [-20333., -9601., 300.34]]  # Mars (units of m/s)
# positions of the planets at (2019/6/2)-1000 days
rspast = [[1.44109e11, -4.45267e10, -509142.],   # Earth
          [1.11393e11, -1.77611e11, -6.45385e9]] # Mars
def motions(t, y):

    rx1,ry1,rz1, rx2,ry2,rz2, vx1,vy1,vz1, vx2,vy2,vz2 = y
    drx1 = vx1
    dry1 = vy1
    drz1 = vz1
    drx2 = vx2
    dry2 = vy2
    drz2 = vz2

    GMmars  = G*Mmars
    GMearth = G*Mearth
    GMsun   = G*Msun

    rx12  = rx1 - rx2
    ry12  = ry1 - ry2
    rz12  = rz1 - rz2
    xyz12 = np.power(np.power(rx12,2) + np.power(ry12,2) + np.power(rz12,2), 1.5)
    xyz1  = np.power(np.power(rx1, 2) + np.power(ry1, 2) + np.power(rz1, 2), 1.5)
    xyz2  = np.power(np.power(rx2, 2) + np.power(ry2, 2) + np.power(rz2, 2), 1.5)

    dvx1 = -GMmars  * rx12 / xyz12 - GMsun * rx1 / xyz1
    dvy1 = -GMmars  * ry12 / xyz12 - GMsun * ry1 / xyz1
    dvz1 = -GMmars  * rz12 / xyz12 - GMsun * rz1 / xyz1
    dvx2 =  GMearth * rx12 / xyz12 - GMsun * rx2 / xyz2
    dvy2 =  GMearth * ry12 / xyz12 - GMsun * ry2 / xyz2
    dvz2 =  GMearth * rz12 / xyz12 - GMsun * rz2 / xyz2

    return np.array([drx1,dry1,drz1, drx2,dry2,drz2,
                     dvx1,dvy1,dvz1, dvx2,dvy2,dvz2])

def find_vel_past():
    daynum=1000
    ts=np.linspace(0,-daynum*day2sec,daynum)
    angles=np.zeros([daynum,2])
    trange =(ts[0],ts[-1])
    fi=np.ndarray.flatten(np.array(rs+vs))
    sol= integrate.solve_ivp(earth_mars_motion,trange,fi,t_eval=ts, max_step=3*day2sec,dense_output=True)
    return(sol.y[0:6][:,-1])
##return an array of six velocities at this time 
def estimate_errors_improved():
    daynum=1000
    ##generating np arrays for bouundary conditions
    a=np.ndarray.flatten(np.array(find_vel_past()))
    rpast=np.ndarray.flatten(np.array(rspast))
    acond=np.concatenate([rpast,a])
    bcond=np.ndarray.flatten(np.array(rs+vs))
    t=np.linspace(0,daynum*day2sec,daynum)
    y=np.zeros(([12,daynum]))
    y[:,0]=acond
    def bc(ya,yb):
        x=yb-bcond
        return np.array(x)
    sol = integrate.solve_bvp(earth_mars_motion1,bc,t,y,verbose=2)
    data1=np.transpose(sol.sol(t))
    angles=np.zeros(daynum)
    for i in range(daynum):      
        angles[i]=angle_between_planets(np.transpose(sol.sol(t)[:,0]))
    x = t/day2sec
    plt.plot(x,angles)
    plt.show()
estimate_errors_improved()

我认为我的代码无法正常工作的原因是由于正在传递的数组的形状存在一些错误。如果有人能告诉我哪里出了问题，我将不胜感激，以便我解决问题。我得到的sol.sol(t) 的输出是：

 Iteration    Max residual  Max BC residual  Total nodes    Nodes added  
Singular Jacobian encountered when solving the collocation system on iteration 1. 
Maximum relative residual: nan 
Maximum boundary residual: 2.14e+11
[[ 1.44109e+11  0.00000e+00  0.00000e+00 ...  0.00000e+00  0.00000e+00
   0.00000e+00]
 [-4.45267e+10  0.00000e+00  0.00000e+00 ...  0.00000e+00  0.00000e+00
   0.00000e+00]
 [-5.09142e+05  0.00000e+00  0.00000e+00 ...  0.00000e+00  0.00000e+00
   0.00000e+00]
 ...
 [         nan          nan          nan ...          nan          nan
           nan]
 [         nan          nan          nan ...          nan          nan
           nan]
 [         nan          nan          nan ...          nan          nan
           nan]]

【问题讨论】：

根据文档，该函数必须返回一个与y 布局相同的数组，即 (12,n) 形状。
您到底想解决什么问题？理论解完全由t=1000 处的位置和速度决定，你可以向后积分到t=0，你得到的任何差异都是数值和测量误差。你得到了什么结果？你打算用边界值求解器计算什么？假设位置准确的速度？
我想在 t=-1000,t=0 区间内计算一个新的插值解，以计算在时间 t=0,t=1000 之间地球和火星未来位置的更准确解.我用solve IVP给我一个t=-1000的试用解决方案。
比什么更准确？目标是某种最小错误计算吗？您是否知道solve_bvp 的标准容差为1e-3，并且通常比IVP 求解器更具实验性？如果您关心求解的准确性，为什么不控制求解器的公差？所涉及的比例与默认公差的合理范围相去甚远。此外，位置和速度的尺度差异几乎没有在可感知的边界上，重新调整求解器可见的数字可能会改善结果。
对于这个问题，我们被特别要求使用solve_bvp，之前已经使用solve_ivp来做同样的事情。我花了很长时间试图弄清楚一些明智的事情，我最近的问题是在求解 bvp_problem 中调用的函数的形状具有在（12,1000）和（12,999）之间波动的数组大小，正如另一个问题中所讨论的那样，我在实施上述第一个建议后发布：stackoverflow.com/questions/58799629/…

标签： python multidimensional-array scipy differential-equations

【解决方案1】：

我认为有几个问题。首先，据我所知，尝试“跑回”到 -1000 天点的唯一原因是获得一个好的 y 估计值以传递给 solve_bvp。

要做到这一点，只需反转初始速度并模拟到 +1000 天。完成此操作后，翻转生成的 sol.y 数组，它们应该可以作为 solve_bvp 的良好估计。

接下来，你实际上不需要 vel 过去，初始位置和 t=0 速度的边界条件就可以完美地完成。

这给我们带来了下一个问题，您的边界条件函数看起来是错误的。

它应该看起来像这样。

def bc(ya, yb):

return np.array([ya[0]-1.44109e11,ya[1] +4.45267e10,ya[2]+509142.,ya[3]-1.11393e11,ya[4]+1.77611e11,ya[5]-6.45385e9,yb[6]-27712.,
                 yb[7]+9730.,yb[8]+0.64148,yb[9]+20333.,yb[10]+9601.,yb[11]-300.34])

最后说明：您很可能不得不将solve_bvp问题中的节点数增加到

希望对你有帮助

【讨论】：