在 python 中使用 bvp 求解器时数组索引超出范围

Question

我在 Python 中使用bvp solver 来模拟d4y/dx4= F(y)/(t-y)^2 形式的四阶边值问题，其中F(y) 是变量y 的复杂函数。我为此编写的代码已附在下面：

import numpy as np
from scipy.integrate import solve_bvp
from numpy import arange
import matplotlib.pyplot as plt
from scipy.optimize import root

def xmesh(k1):
 return np.linspace(0,L,k1)

def solveit(gamma):
 def fun(x,y): # this function returns all the derivatives for the bvp solver
  m1=y[0] # m1 is an 1-D array
  array=np.empty(10,)
  def calc(): # this function calculates F(y) at all the mesh points and finally returns an array 
    for i in range(10) :
     def fun2(var):
       return gamma*(var)+var**(0.5)+var**(3)+m1[i]
       
     sol= root(fun2,x0=0.00001,method='krylov',tol=1e-12)  
     array[i]=sol.x
         
    return array   

  m2=calc() # m2 is an 1-D array and is equal to F(y)
  rhs= m2/((t-m1)**2) # rhs contains the RHS of the differential equation
  return np.vstack([y[1],y[2],y[3],rhs])
  
 def bc(ya,yb): # boundary conditions for the bvp solver
    return np.array([ya[0],ya[1],yb[0],yb[1]])
 q= 10
 yinit= np.zeros((4,len(xmesh(q))))
 sol= solve_bvp(fun,bc,xmesh(q),yinit,tol=1e-6,max_nodes=5000)
 arr= sol.sol(xmesh(q))[0] # array containing y
 return arr   
a=solveit(0.1) 

每当我尝试运行上述代码时，都会遇到错误：index 9 is out of bounds for axis 0 with size 9。我不明白为什么会出现这个错误。这意味着数组有 9 个元素，我试图引用数组的第 10 个元素，因此会引发此错误。但是上面代码中定义的所有数组每个都有 10 个元素。所以，我不明白为什么甚至会引发这个错误。我已经坚持了2天了。对上述问题的任何解决方案都将受到高度赞赏。 P.S.- 我在上面的代码中定义的函数没有物理意义，完全是随机的。因此，我不需要上述代码的正确数值解。我所需要的只是代码可以正常工作而不会出现任何错误。另外，我之前在 Python 中成功使用了 bvp 求解器。

edit：我已经添加了堆栈跟踪。

IndexError                                Traceback (most recent call last)
<ipython-input-4-8423177548fe> in <module>()
     34  arr= sol.sol(xmesh(q))[0] # array containing y
     35  return arr
---> 36 a=solveit(0.1)

12 frames
<ipython-input-4-8423177548fe> in solveit(gamma)
     31  q= 10
     32  yinit= np.zeros((4,len(xmesh(q))))
---> 33  sol= solve_bvp(fun,bc,xmesh(q),yinit,tol=1e-6,max_nodes=5000)
     34  arr= sol.sol(xmesh(q))[0] # array containing y
     35  return arr

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in solve_bvp(fun, bc, x, y, p, S, fun_jac, bc_jac, tol, max_nodes, verbose, bc_tol)
   1084                                        fun_jac_wrapped, bc_jac_wrapped, x, h)
   1085         y, p, singular = solve_newton(n, m, h, col_fun, bc_wrapped, jac_sys,
-> 1086                                       y, p, B, tol, bc_tol)
   1087         iteration += 1
   1088 

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in solve_newton(n, m, h, col_fun, bc, jac, y, p, B, bvp_tol, bc_tol)
    439     n_trial = 4
    440 
--> 441     col_res, y_middle, f, f_middle = col_fun(y, p)
    442     bc_res = bc(y[:, 0], y[:, -1], p)
    443     res = np.hstack((col_res.ravel(order='F'), bc_res))

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in col_fun(y, p)
    324 
    325     def col_fun(y, p):
--> 326         return collocation_fun(fun, y, p, x, h)
    327 
    328     def sys_jac(y, p, y_middle, f, f_middle, bc0):

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in collocation_fun(fun, y, p, x, h)
    311     y_middle = (0.5 * (y[:, 1:] + y[:, :-1]) -
    312                 0.125 * h * (f[:, 1:] - f[:, :-1]))
--> 313     f_middle = fun(x[:-1] + 0.5 * h, y_middle, p)
    314     col_res = y[:, 1:] - y[:, :-1] - h / 6 * (f[:, :-1] + f[:, 1:] +
    315                                               4 * f_middle)

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in fun_p(x, y, _)
    648     if k == 0:
    649         def fun_p(x, y, _):
--> 650             return np.asarray(fun(x, y), dtype)
    651 
    652         def bc_wrapped(ya, yb, _):

<ipython-input-4-8423177548fe> in fun(x, y)
     23     return array
     24 
---> 25   m2=calc() # m2 is an 1-D array and is equal to F(y)
     26   rhs= m2/((t-m1)**2) # rhs contains the RHS of the differential equation
     27   return np.vstack([y[1],y[2],y[3],rhs])

<ipython-input-4-8423177548fe> in calc()
     18        return gamma*(var)+var**(0.5)+var**(3)+m1[i]
     19 
---> 20      sol= root(fun2,x0=0.00001,method='krylov',tol=1e-12)
     21      array[i]=sol.x
     22 

/usr/local/lib/python3.6/dist-packages/scipy/optimize/_root.py in root(fun, x0, args, method, jac, tol, callback, options)
    197         sol = _root_nonlin_solve(fun, x0, args=args, jac=jac,
    198                                  _method=meth, _callback=callback,
--> 199                                  **options)
    200     else:
    201         raise ValueError('Unknown solver %s' % method)

/usr/local/lib/python3.6/dist-packages/scipy/optimize/_root.py in _root_nonlin_solve(func, x0, args, jac, _callback, _method, nit, disp, maxiter, ftol, fatol, xtol, xatol, tol_norm, line_search, jac_options, **unknown_options)
    300                                   line_search=line_search,
    301                                   callback=_callback, full_output=True,
--> 302                                   raise_exception=False)
    303     sol = OptimizeResult(x=x)
    304     sol.update(info)

/usr/local/lib/python3.6/dist-packages/scipy/optimize/nonlin.py in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)
    277 
    278     dx = np.inf
--> 279     Fx = func(x)
    280     Fx_norm = norm(Fx)
    281 

/usr/local/lib/python3.6/dist-packages/scipy/optimize/nonlin.py in <lambda>(z)
    273 
    274     x0 = _as_inexact(x0)
--> 275     func = lambda z: _as_inexact(F(_array_like(z, x0))).flatten()
    276     x = x0.flatten()
    277 

<ipython-input-4-8423177548fe> in fun2(var)
     16     for i in range(10) :
     17      def fun2(var):
---> 18        return gamma*(var)+var**(0.5)+var**(3)+m1[i]
     19 
     20      sol= root(fun2,x0=0.00001,method='krylov',tol=1e-12)

IndexError: index 9 is out of bounds for axis 0 with size 9

Answer 1

在函数中

  def calc(): # this function calculates F(y) at all the mesh points and finally returns an array 
    for i in range(10) :
        def fun2(var):
            return gamma*(var)+var**(0.5)+var**(3)+m1[i]

您正在使用i 在嵌套函数内迭代m1，但在上面的代码中

m1=y[0] # m1 is an 1-D array

您将此值定义为无法迭代的标量（或一维数组）。

这会给你你正在观察的错误；对于您定义的一维数组m1，任何大于0 的索引都超出了此一维数组的范围。

这应该为您提供更改代码以使其正常工作所需的内容。

我不会为你完成这个，因为我们不知道你需要做什么，但我希望在 for 循环中一次又一次地定义一个函数是不必要的。