双直接积分

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 solve 返回的错误是我只有一个雅可比行列式。当我删除最后两个方程时，我得到了 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)

但是，我得到的错误是 'setting an array with sequence'。第一个函数和边界条件求解两个耦合方程，然后我使用这些结果来评估我给出的方程。我曾尝试将我所有的方程式写在一个函数中，但这似乎返回了平凡的解决方案，即一个全为零的数组。

如有任何帮助，我们将不胜感激。

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

双直接积分

Double Direct Integration

python

integration

boundary

ode

differential-equations