Решая ODE с помощью scipy.integrate.solve_bvp, возврат моей функции не может транслироваться в соответствующую форму массива для resolve_bvp - PullRequest
Новая
       63

Решая ODE с помощью scipy.integrate.solve_bvp, возврат моей функции не может транслироваться в соответствующую форму массива для resolve_bvp

0 голосов
/ 13 июля 2020

Я пытаюсь решить ODE второго порядка с помощью resolve_bvp. Я разделил ОДУ второго порядка на систему буксировочных ОДУ первого порядка. У меня есть изменяющийся набор констант в зависимости от значения x (me sh). Итак, я передаю их как массив формы (N,) в мою функцию numdens. При попытке запустить resolve_bvp я получаю сообщение об ошибке, что возвращаемые значения имеют разные формы, а именно (N,) и (N-1,) и, следовательно, не могут быть переданы в один массив. Но когда я проверяю каждый возврат вручную вне функции, он имеет форму (N,). Если я запустил решающую программу без изменения констант, я получу решение, похожее на правильное. 

import numpy as np
from scipy.integrate import solve_bvp,odeint
import matplotlib.pyplot as plt

E_0 = 1 * 0.0000016021773 #erg: gcm^2/s^2
m_H = 1.6*10**(-24) #g
c = 3e11 #cm
sigma_c = 2*10**(-23)
n_0 = 1*10**(20)   #1/cm^3
v_0 = (2*E_0/m_H)**(0.5)  #cm/s
T = 10**7
b = 20.3
n_eq = b*T**3
n_s = 2.03*10**(19)
Q = 1 

def velocity(v,x):
    dvdx = -sigma_c*n_0*v_0*((8*v_0*v-7*v**2-v_0**2)/(2*v*c))
    return dvdx


n_num = 100
x_num = np.linspace(-1*10**(6),3*10**(6), n_num)

sol_velo = odeint(velocity,0.999999999999*v_0,x_num)

sol_new = np.reshape(sol_velo,n_num)

def constants(v):
    D1 = (c*v/(3*n_0*v_0*sigma_c))
    D2 = ((v**2-8*v_0*v+v_0**2)/(6*v))
    D3 = sigma_c*n_0*v_0*((8*v_0*v-7*v**2-v_0**2)/(2*v*c))
    return D1,D2,D3

def numdens(x,y):
    v = sol_new
    D1,D2,D3 = constants(v)
    return np.vstack((y[1],(-D2*y[1]-D3*y[0]+Q*((1-y[0])/n_eq))/(D1)))

def bc_num(ya, yb):
    return np.array([ya[0]-n_s,yb[0]-n_eq])


y_num = np.array([np.linspace(n_s, n_eq, n_num),np.linspace(n_s, n_eq, n_num)])


sol_num = solve_bvp(numdens, bc_num, x_num, y_num)




plt.plot(sol_num.x, sol_num.y[0], label='$n(x)$')
plt.plot(x_num, sol_velo-v_0/7, label='$v(x)$')
plt.yscale('log')
plt.grid(alpha=0.5)
plt.legend(framealpha=1)
plt.show()

1 Ответ

0 голосов
/ 13 июля 2020

Необходимо принять во внимание, что решатель BVP использует адаптивный меня sh. То есть, после уточнения начального предположения на исходной сетке решатель определяет области с чрезмерно большими ошибками и создает там новые узлы me sh. Насколько я видел, обратное не реализовано, даже если в некоторых приложениях может быть разумным уменьшить количество узлов me sh на особенно «хороших» сегментах.

Итак, что вы делаете функция numdens непонятна, она должна работать точно так же, как любая другая функция, которую вы передали бы решателю ODE. Если бы мне пришлось предложить какое-то быстрое исправление и не зная, в чем заключается основная проблема, которую вы хотите решить, я бы изменил присвоение v на 

v = np.interp(x,x_num,sol_velo)

, поскольку это должно по крайней мере создать массив правильного формата.

...