Решение набора краевых задач

Я пытаюсь решить набор краевых задач, задаваемых 4 дифференциальными уравнениями. Я использую bvp_solver в python и получаю сообщения об ошибках, в которых указано «недопустимое значение, обнаруженное при делении». Я предполагаю, что это означает, что в какой-то момент я делю на NaN или 0, но я не уверен, где.

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
%matplotlib inline
alpha = 1
zeta = 1
C_k = 1
sigma = 1
Q = 30
U_0 = 0.1
gamma = 5/3
theta = 3
m = 1.5
def fun(x, y):
    U, dU, B, dB, T, dT, M, dM = y;
    d2U = -2*U_0*Q**2*(1/np.cosh(Q*x))**2*np.tanh(Q*x)-((alpha)/(C_k*sigma))*dB;   
    d2B = -(1/(C_k*zeta))*dU;
    d2T = (1/gamma - 1)*(sigma*dU**2 + zeta*alpha*dB**2);
    d2M = -(dM/T)*dT + (dM/T)*theta*(m+1) - (alpha/T)*B*dB
    return dU, d2U, dB, d2B, dT, d2T, dM, d2M 

def bc(ya, yb):

    return ya[0]+U_0*np.tanh(Q*0.5), yb[0]-U_0*np.tanh(Q*0.5), ya[2]-0, yb[2]-0, ya[4] - 1, yb[4] - 4, ya[6], yb[6] - 1

x = np.linspace(-0.5, 0.5, 500)
y = np.zeros((8, x.size))


sol = solve_bvp(fun, bc, x, y)

Если я удалю два последних уравнения для M и dM, решение будет работать нормально. В прошлом у меня были проблемы с пониманием возвращаемых массивов bvp_solver, но теперь я уверен, что понимаю это. Но я продолжаю получать ошибки каждый раз, когда добавляю новые уравнения. Любая помощь приветствуется.


python numpy ode differential-equations boundary
person Patrick Lewis    schedule 28.02.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Конечно, на первом этапе это не удастся. Вы инициализируете все до нуля, а затем в функции производных делите на T, который равен нулю от инициализации.

  • Найдите более реалистичную инициализацию T, например

    x = np.linspace(-0.5, 0.5, 15)
y = np.zeros((8, x.size))
y[4] = 2.5+3*x
y[5] = 3+0*x

    or

  • десингуляризация деления, что обычно делается аналогично

    d2M = (-dM*dT + dM*theta*(m+1) - alpha*B*dB) * T/(1e-12+T*T)

Всегда имеет смысл печатать после sol = solve_bvp(...) сообщения об ошибке print(sol.message). Теперь, когда компонентов больше, чем несколько, я изменил конструкцию выходной таблицы на систематическую.

%matplotlib inline
plt.figure(figsize=(10,2.5*4))
for k in range(4):
    v,c = ['U','B','T','M'][k],['-+b','-*r','-xg','-om'][k];
    plt.subplot(4,2,2*k+1); plt.plot(sol.x,sol.y[2*k  ],c, ms=2); plt.grid(); plt.legend(["$%c$"%v]);
    plt.subplot(4,2,2*k+2); plt.plot(sol.x,sol.y[2*k+1],c, ms=2); plt.grid(); plt.legend(["$%c'$"%v]);
plt.tight_layout(); plt.savefig("/tmp/bvp3.png"); plt.show(); plt.close()
person Lutz Lehmann    schedule 28.02.2020
comment
Я решал эти уравнения, но, тем не менее, я заметил, что мое решение для M и dM меняется, когда я меняю количество точек сетки в x, которые я передаю в bvp_solve. Мне не ясно, прихожу ли я к истинному решению или что именно происходит. - person Patrick Lewis; 03.03.2020
comment
Если решающая программа вернется успешно, то, скорее всего, она близка к решению. Обратите внимание, что БВП, нелинейная в любой из своих частей, может иметь несколько решений. Первоначальная догадка влияет на то, какой из них будет достигнут. - person Lutz Lehmann; 03.03.2020
comment
Спасибо, это имеет смысл. Как мне вывести успешное сообщение из моего кода? - person Patrick Lewis; 05.03.2020
comment
У вас есть логический флаг sol.success, который вы можете использовать для остановки дальнейшей обработки в случае сбоя, и текстовое сообщение sol.message, которое вы можете поместить в оператор печати. - person Lutz Lehmann; 05.03.2020