“简单”微分方程系统期间solve_bvp 中的溢出

Question

我想求解以下微分方程组：

d(phi)/dz= Rcl*i
d(i)/dz = i0*exp(phi*2.3/b)

我尝试使用 scipy 中的solve_bvp 来解决这个问题，但没有成功。为了检查是否可以找到解决方案，我编写了自己的代码并使用了射击方法（颂歌系统的初始值求解器与优化相结合以找到起始值）并且它起作用了。我想建立一个更复杂的微分方程组，因此如果我可以使用内置函数那就太好了。你能看到我使用solve_bvp做错了什么吗？或者你知道改进算法的方法吗？我的代码是：

solve_bvp：

from scipy.integrate import solve_bvp
import numpy as np
import matplotlib.pyplot as plt
#define parameters
b = 50e-3
Rcl = 0.4
i0 = 1e-7

il = 0
ir = 1
#define ode system
def fun(z,y):
    phi = y[0]
    i = y[1]
    didz = i0*np.exp(phi*2.3/b)
    dphidz = Rcl*i
    dydz = np.array([dphidz,didz])
    return dydz
#Define boundary conditions
def boundary_conditions(ya, yb):

    return np.array([ya[1]-il, yb[1]-ir])

#define numerical aprameters
a = 0
b = 1
num_points = 100

#Define the initial mesh and initial guess for the solution
x_mesh = np.linspace(a, b, num_points)  # Define the interval [a, b] and the number of points
y_initial_guess = np.zeros((2, x_mesh.size))  # Initial guess for the solution
solution = solve_bvp(fun, boundary_conditions, x_mesh, y_initial_guess)
#plot
#Evaluate the solution on a finer mesh for plotting
x_fine = np.linspace(a, b, 1000)
y_fine = solution.sol(x_fine)
#%%
#Plot the solution
fig, ax = plt.subplots()
plt.plot(x_fine, y_fine[0], label='Phi')
plt.legend()
fig, ax = plt.subplots()
plt.plot(x_fine, y_fine[1]/Rcl, label='J')
plt.legend()

自己的代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
#define parameters
l_cl = 1       # [cm]
R_cl = 0.8     # Ohm/cm²
b = 50e-3      # 1/V
i0 = 1e-7      # A/cm²
il
ir = 1         # Target current density

num_points = 1000
step_width = l_cl / num_points
#define function to solve ode system with a given starting value
def fun(phi_start):
    state_vars = np.zeros((3, num_points))  # Single array for state variables (phi, dphidx, i)
    state_vars[0, 0] = phi_start     #define start value for phi
    state_vars[2,0] = il             #define start value for i
#calculate solution of the ode system with phi_start and give out difference between i at right boundary and ir
    for ii in range(num_points - 1):
        didx = i0 * np.exp((state_vars[0, ii]) / b * np.log(10))                        #didx
        state_vars[2, ii + 1] = state_vars[2, ii] + didx * step_width                   #i
        state_vars[1, ii + 1] = state_vars[2, ii + 1] * R_cl                            #dphidx
        state_vars[0, ii + 1] = state_vars[0, ii] + state_vars[1, ii + 1] * step_width  #phi
    return abs(state_vars[2, -1] - ir)
#find phi_start which aims to i at right boundary equals ir
x0 = 0.31
res = minimize(fun, x0, method='Nelder-Mead', tol=1e-6)
phi_start_optimized = res.x
print(res.x)

state_vars = np.zeros((3, num_points))  # Single array for state variables (phi, dphidx, i)
state_vars[0, 0] = res.x
mesh = np.linspace(0,1,num_points)
#solve everything with phi_start to show solution
for ii in range(num_points - 1):
    didx = J0 * np.exp((state_vars[0, ii]) / b * np.log(10))
    state_vars[2, ii + 1] = state_vars[2, ii] + didx * step_width
    state_vars[1, ii + 1] = state_vars[2, ii + 1] * R_cl
    state_vars[0, ii + 1] = state_vars[0, ii] + state_vars[1, ii + 1] * step_width

phi = state_vars[0, :]
dphidx = state_vars[1, :]
i = state_vars[2, :]
#plot
plt.plot(mesh,phi,label='over potential')
plt.plot(mesh,i,label='i_sol')
plt.legend()

非常感谢

Answer 1

感谢您的评论。完整的错误消息是：“c:... RuntimeWarning: 在 exp didz = J0np.exp(phi2.3/b) 中遇到溢出 C:... RuntimeWarning: 在减法中遇到无效值 0.125 * h * (f [:, 1:] - f[:, :-1])) C:... RuntimeWarning: 减法中遇到无效值 df_dy[:, i, :] = (f_new - f0) / hi C:... RuntimeWarning：除法中遇到无效值 r_middle /= 1 + np.abs(f_middle)" 但是，我知道颂歌的解决方案，并且两个变量都在 0 和 1 之间的范围内 - 所以不应该溢出... 切割指数（按照 Lutz Lehmann 的建议）不起作用，我收到以下错误消息：“ValueError：具有多个元素的数组的真值不明确。使用 a.any() 或 a.all() “求解器似乎可以同时处理多个解决方案，但无法处理这个问题......