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

问题描述

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

  from datetime import datetime 
进口matplotlib.pyplot as plt 
％matplotlib内联
进口numpy as np 
从scipy进口整合
从scipy进口信号
 
 G = 6.67408e- 11＃m ^ 3 s ^ -1 kg ^ -2 
 AU = 149.597e9＃m 
土= 5.9721986e24＃kg 
 Mmars = 6.41693e23＃kg 
 Msun = 1.988435 e30＃公斤
 day2sec = 3600 * 24＃一天中的秒数
 
 rs = [[-4.8957151e10，-1.4359284e11，501896.65]，＃地球
 [-1.1742901e11 ，2.1375285e11，7.3558899e9]]＃火星（m单位）
 vs = [[27712.，-9730。-0.64148]，＃地球
 [-20333。，-9601。，300.34 ]]＃火星（m / s单位）
＃行星在（2019/6/2）-1000天的位置
 rspast = [[1.44109e11，-4.45267e10，-509142。] ，＃Earth 
 [1.11393e11，-1.77611e11，-6.45385e9]]＃火星
 def运动（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（nx.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 = GMears / xyz12-GMsun * rx2 / xyz2 
 dvy2 = GMearth * ry12 / xyz12-GMsun * ry2 / xyz2 
 dvz2 = GMearth * rz12 / xyz12-GMsun * rz2 / xyz2 
 
 $返回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]）$ ​​b $ b fi = np.ndarray.flatten（np.array（rs + vs））
 sol = integration.solve_ivp（earth_mars_motion，trange，fi，t_eval = ts，max_step = 3 * day2sec，dense_output = True）
 return（sol.y [0：6] [：，-1]）$ ​​b $ b ##此时返回六个速度数组
 def Estimate_errors_improved（）：
 daynum = 1000 
 ##为有条理的条件生成np数组
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 
返回np.array（x）
 sol =积分.solve_bvp（earth_mars_motion1，bc，t，y，verbose = 2）
 data1 = np.transpose（sol.sol（t））
 angles = np.zeros（daynum）$ b i在范围（daynum）中的$ b：
 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）的输出是：

 迭代最大残差最大BC残差总节点已添加的节点
在迭代1中求解配置系统时遇到的奇异Jacobian。
最大相对残差：nan 
最大边界残差：2.14 e + 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]] 
  

解决方案

一些问题nk。首先，据我所知，试图跑回 -1000天的唯一原因是获得一个良好的y估计值，然后传递给solve_bvp。

为此，只需反转初始速度，然后模拟到+1000天即可。完成此操作后，翻转所得的sol.y数组，即可将其作为solve_bvp的理想估计。

接下来，您实际上不需要滑行，初始位置的边界条件和t = 0速度将非常完美。

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

它应该看起来像这样。

\\

  def bc（ya，yb）：
 
返回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]）
   
 
 
 \\ 
 
 
 
 resolve_bvp问题
 
 
 
希望这会有所帮助
 
This is a very general question as I feel that my errors are resulting from some misunderstanding of how scipy.solve_bvp works. I have a function def that takes an array of 12 numbers and returns a list of the system of differential equations for a given time, with shape (2,6). I will have a one dimensional array of length n for my timesteps and then an array yof input values with shape (12,n). My code aims to take simulate the motion of earth and mars over a 1000 day period subject to boundary conditions; at t=0 positions = rpast (the corresponding velocities are returned by the function find_vel_past()), the positions and velocities at t=1000 are given by rs and vs respectively. My code is at the bottom with the two functions I'm trying to solve above:
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()
I think that the reason my code isnt working is due to some error in the shapes of arrays that are being passed. I would be very grateful if someone could tell me where I am going wrong so I can fix things.The output for sol.sol(t) I'm getting is:
 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]]
 解决方案 
a few problems I think. Firstly the only reason for trying to 'run back' to the -1000 day point as far as I can see would be to obtain a good y estimate to pass to solve_bvp.
to do this simply reverse the initial velocities and run a similation to +1000 days. once you have done this flip the resulting sol.y arrays and they should serve as a good estimate for solve_bvp.  
Next, you dont actually need vel past, boundary conditions of the initial position and the t=0 velocity will do perfectly. 
This brings us to the next problem, your Boundary condition function looks mistaken. 
it should look something like this.
\\
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])
\\ 
Final note: you will most likley have to increase the number of nodes in the solve_bvp problem to
hope this helps
                        
这篇关于使用scipy.solve_bvp解决BVP，其中函数返回一个数组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！