Double Direct Integration

Question

我正在尝试解决一系列耦合的边值问题；

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)

但是，我得到的错误是“用序列设置数组”。第一个函数和边界条件求解两个耦合方程，然后使用这些结果评估我给出的方程。我试图将所有方程式写在一个函数中，但是这似乎返回了平凡的解，即一个由零组成的数组。

任何帮助将不胜感激。

Answer 1

当表达式变大时，使计算易于人读而不是紧凑通常会更有帮助。

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个点进行细化）绘制为