Двойная прямая интеграция - PullRequest
Новая
       91

Двойная прямая интеграция

1 голос
/ 27 февраля 2020

Я пытаюсь решить множество связанных краевых задач таким образом, чтобы 

U'' +aB'+ b*(cosh(lambda z))^{-2}tanh(lambda*z) = 0,

B'' + c*U' = 0, 

T'' = (gamma^{-1} - 1)*(d*(U')^2 + e*(B')^2)

с учетом граничных условий U(+/- 1/2) = +/-U_0*tanh(lambda/2), B(+/- 1/2) = 0 and T(-1/2) = 1, T(1/2) = 4. Я разложил этот набор уравнений в набор дифференциальных уравнений первого порядка и использовал массив производных так, чтобы [U, U', B, B', T, T']. Но BVP решить возвращает ошибку, что у меня есть один якобиан. Когда я удаляю последние два уравнения, я получаю решение для U и B, и это прекрасно работает. Однако я не уверен, почему добавление двух других уравнений приводит к этой проблеме. 

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
%matplotlib inline
alpha = 1E-7
zeta = 8E-3
C_k = 0.01 
sigma = 0.005
Q = 30
U_0 = 0.1
gamma = 5/3
theta = 3

def fun(x, y):

    return y[1], -2*U_0*Q**2*(1/np.cosh(Q*x))**2*np.tanh(Q*x)-((alpha)/(C_k*sigma))*y[3],  y[3],\
        -(1/(C_k*zeta))*y[1], y[4], (1/gamma - 1)*(sigma*(y[1])**2 + zeta*alpha*(y[3])**2)

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

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


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

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

Буду признателен за любую помощь.

1 Ответ

2 голосов
/ 27 февраля 2020

Когда выражения становятся больше, часто бывает удобнее сохранять вычисления удобочитаемыми, а не компактными. 

def fun(x, y):
    U, dU, B, dB, T, dT = 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);
    return dU, d2U, dB, d2B, dT, d2T

Это позволяет избежать пропуска ошибки индекса, поскольку в этом вычислении нет индексов, все имеют имена близко к исходным формулам.

Тогда компоненты решения (используя инициализацию только с 5 точками, что приводит к уточнению с 65 точками) отображают как

enter image description here

...