解决一组边值问题 [英] Solving set of Boundary Value Problems

问题描述

我正在尝试解决由4个微分方程给出的一组边值问题.我在python中使用bvp_solver，并且收到错误消息，指出除法中遇到无效的值".我假设这意味着我在某个时候被NaN或0除，但是我不确定在哪里.

I am trying to solve a set of boundary value problems given by 4 differential equations. I am using bvp_solver in python, and I am getting errors which state 'invalid value encountered in division'. I am assuming this means I am dividing by NaN or 0 at some point, but I am unsure where.

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的返回数组时遇到了麻烦，但是我有把握，现在我已经知道了.但是，每次添加更多方程式时，我仍然会不断出错.任何帮助，我们将不胜感激.

If I remove the last two equations for M and dM, then the solution works fine. I have had trouble in the past understanding the return arrays of bvp_solver, but I am confident I understand that now. But I continue to get errors each time I add more equations. Any help is greatly appreciated.

推荐答案

当然，这第一步将失败.将所有内容初始化为零，然后在派生函数中，除以 T ，该值从初始化开始为零.

Of course this will fail in the first step. You initialize everything to zero, and then in the derivatives function, you divide by T, which is zero from the initialization.

• 例如，找到更实际的 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

或

将除法分解为整数，通常与

desingularize the division, which usually is done similar to

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

在错误代码 print(sol.message)后，在 sol = solve_bvp(...)之后打印总是有意义的.既然有更多的组件，我将输出表格的结构更改为系统的

It makes always sense to print after sol = solve_bvp(...) the error message print(sol.message). Now that there are more than a few components, I changed the construction of the output tableau to the systematic

%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()

这篇关于解决一组边值问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！