Решение BVP с использованием scipy.solve_bvp, где функция возвращает массив - PullRequest
Новая
       24

Решение BVP с использованием scipy.solve_bvp, где функция возвращает массив

0 голосов
/ 11 ноября 2019

Это очень общий вопрос, так как я чувствую, что мои ошибки происходят из-за неправильного понимания того, как работает scipy.solve_bvp. У меня есть функция def, которая принимает массив из 12 чисел и возвращает список системы дифференциальных уравнений для заданного времени с формой (2,6). У меня будет одномерный массив длины n для моих временных шагов, а затем массив y входных значений с формой (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]]

1 Ответ

0 голосов
/ 12 ноября 2019

несколько проблем, я думаю. Во-первых, единственная причина, по которой я пытаюсь «вернуться назад» к точке дня -1000, насколько я вижу, состоит в том, чтобы получить хорошую оценку y для перехода к solve_bvp.

, чтобы сделать это просто, чтобы полностью изменить начальные скорости. и запустить аналог до +1000 дней. как только вы это сделаете, переверните получившиеся массивы sol.y, и они должны послужить хорошей оценкой для solve_bvp.

Далее, вам на самом деле не нужно проходить мимо, граничные условия начальной позиции и скорости 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 на

надеюсь, это поможет

...