带有自适应解决方案的 scipy solve_ivp

如何解决带有自适应解决方案的 scipy solve_ivp

我正在努力理解 scipy.solve_ivp() 如何处理 ODE 系统中的错误。假设我有以下单个 ODE 的简单代码，我想我可能在某些方面做错了。假设我的 rhs 看起来像：

from scipy.integrate import solve_ivp
def rhs_func(t,y):  
    z = 1.0/( x - y + 1j )
    return z

假设我们使用以下签名调用solve_ivp：

Z_solution  = ivp_adaptive.solve_ivp(    fun    = rhs_func,t_span = [100,0],y0     = y0,#some initial value of 0 for example  
                                         method ='RK45',t_eval = None,args   = some_additional_arguments_to_rhs_func,dense_output = False,rtol   = 1e-8
                                         atol   = 1e-10
                                    )

现在，应该使用绝对和相对容差来修正计算错误。在这种情况下，我遇到的问题与“t_eval=None”有关。显然，这个选择让积分器（在这种情况下是 RK45 类型）根据上面指定的容差选择时间步长是否被超过，即步长不是固定的，但不知何故在 t 中采取更大的步长意味着已经找到了低于上述容差的解决方案 (atol=1e-10,rtol=1e-8)。这在时间尺度变化很大的问题中特别有用，在这种情况下，t 的均匀离散化非常低效。

我的大问题与 scipy.integrate._ivp.solve_ivp() 中第 575 行附近的以下代码段有关，其中包含“t_eval == None”情况：

while status is None:
    message = solver.step()

    if solver.status == 'finished':
        status = 0
    elif solver.status == 'Failed':
        status = -1
        break

    t_old = solver.t_old
    t = solver.t
    y = solver.y

    if dense_output:
        sol = solver.dense_output()
        interpolants.append(sol)
    else:
        sol = None

    if events is not None:
        g_new = [event(t,y) for event in events]
        active_events = find_active_events(g,g_new,event_dir)
        if active_events.size > 0:
            if sol is None:
                sol = solver.dense_output()

            root_indices,roots,terminate = handle_events(
                sol,events,active_events,is_terminal,t_old,t)

            for e,te in zip(root_indices,roots):
                t_events[e].append(te)
                y_events[e].append(sol(te))

            if terminate:
                status = 1
                t = roots[-1]
                y = sol(t)

        g = g_new
    
    # HERE I HAVE MODIFIED THE FILE BY CALLING AN INTERPOLATION FUNCTION FOR THE SOLUTION
    if t_eval is None:
        ts.append(t)
        #ys.append(y)
        # this calls to adapt the solution to a new set of values x over which y(x,t) is 
        # defined
        interp_solution(t,y,solver,args)
        
        y = solver.y
        ys.append(y)

我在其中定义了一个函数：

def interp_solution( t,args ):
    import numpy as np
    from scipy import interpolate

    x_old          =   args.get_old_grid() # this call just returns an array of the style of 
                                           # x_new,and is where y is defined 
    x_new          =   np.linspace( -t,t,dim ) # the new array where components of y are 
                                                 # defined 
    y_interp       =   interpolate.interp1d( x_old,y )
    y_new          =   y_interp( x_new )
    
    solver.y       =   y_new # update the solver y

    # finally,we change the maximum allowed step of the integrator if t is below
    # some threshold value

    if ( t < args.get_threshold() ):
        solver.max_step = #some number
    

    return y_new

当我查看结果时，似乎这对容差和集成步骤的执行方式非常敏感，但不知何故，我看不出这种方法中的错误可能来自哪里——谁能解释一下方法以某种方式影响解决方案和相关错误？如何以这种方式实施类似的方法？非常感谢任何帮助。