参数 space 中 solve_ivp 的问题

Question

我在创建 Matlab 中给出的示例 ODE 以使用 scipy 的 solve_ivp 时遇到问题。在 Matlab 中，函数定义为

function fixed_point_linear_center()
    clc; clf;

    stepsize=.5;
    xmin=-5;
    xmax=5;
    ymin=-5;
    ymax=5;

    [x,y] = meshgrid(xmin:stepsize:xmax,ymin:stepsize:ymax);

    A = [0 1;-1 0];

    dx = A(1,1)*x + A(1,2)*y;
    dy = A(2,1)*x + A(2,2)*y;

    % Strange scaling for nicer output, only "cosmetics"
    eunorm = ( dx.^2 + dy.^2 ).^(0.35);
    dx = dx./eunorm;
    dy = dy./eunorm;

    quiver(x,y,dx,dy);
    axis([xmin xmax ymin ymax]);
    grid on; xlabel('x'); ylabel('y');

    tspan=[0 100];

    x0stepsize=0.25;
    for x0=xmin:x0stepsize:xmax
        hold on
        ic = [x0 0];
        [~,x] = ode45(@(t,x) f(t,x,A),tspan,ic);
        plot(x(:,1),x(:,2),'r');
        hold on
        ic = [0 x0];
        [~,x] = ode45(@(t,x) f(t,x,A),tspan,ic);
        plot(x(:,1),x(:,2),'r');
    end
    hold off
end

function dx = f(~,x,A)
    dx = A*[x(1); x(2)];
end

计算如下所示的解决方案

，但是如果我像这样在 python 中重新创建函数

def fixed_point_linear_center():
    stepsize   =  0.5
    x0stepsize =  0.25
    xmin       = -5
    xmax       =  5
    ymin       = -5
    ymax       =  5
    
    x    = np.arange(xmin, xmax+stepsize, stepsize)
    xval = np.arange(xmin, xmax+x0stepsize, x0stepsize)
    y    = np.arange(ymin, ymax+stepsize, stepsize)
    yval = np.arange(ymin, ymax+stepsize*0.25, stepsize*0.25) # evaluate 4 times for smoothness
    
    [X, Y] = np.meshgrid(x, y)
    
    
    A = np.array([[0,1],[-1,0]])

    dx = A[0,0]*X + A[0,1]*Y # 21x21
    dy = A[1,0]*X + A[1,1]*Y # 21x21
    
    f = lambda t,x,A : np.dot(A,[[x[0]],[x[1]]])

    # Strange scaling for nicer output, but only "cosmetics"
    eunorm = np.float_power(( dx**2 + dy**2 ), 0.35)   #( dx**2 + dy**2 )**0.35
    eunorm[10,10] = 0.001 # center is 0 which violates division
    dx = dx/eunorm
    dy = dy/eunorm

    plt.figure(figsize = (15,12))
    plt.quiver(X, Y, dx, dy, angles = 'xy', color='#0086b3', width=0.0015)
    plt.grid() 
    plt.xlabel('x') 
    plt.ylabel('y')

    plt.axis([xmin,xmax,ymin,ymax])
    
    tspan=[0,100]
    
    for x0 in xval:
        
        ic = [x0,0]
        #[~,x] = ode45(@(t,x) f(t,x,A),tspan,ic);
        solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(A,))
        #solution = solve_ivp(f, [xmin, xmax], [x0], method='RK45', t_eval=yval, dense_output=False, args=(0,A))
        #solution = solve_ivp(f, [tmin, tmax], [ic], method='RK45', t_eval=tval, args=(A), dense_output=False)
        plt.plot(solution.y[1], solution.y[0],'r')
        
fixed_point_linear_center()

我收到类似

的错误

ValueError: shapes (2,2) and (2,1,2) not aligned: 2 (dim 1) != 1 (dim 1)

或类似的，取决于我已经尝试将 f 重写为什么。据我了解，solve_ivp 需要 x0 数组中的单个值，而我 return 是一个 2x1 向量。它也不接受向量作为 x0 数组中的值，例如 [[x0,0]]

现在我想知道 scipy.solve_ivp 是否能够像 ode45 那样对参数 space 进行计算（我该怎么做）或者我是否必须以其他方式进行计算？

（我已经检查过，所有其他矩阵和 return 值与 matlab 计算相同。）

[编辑 2] 好的，现在可以了。 x 的绘图参数当然必须是 solution.y[1]！

Answer 1

就像 Matlab 求解器一样，solve_ivp 期望状态是单个向量。变化

f = lambda t,x1,x2, A : np.dot(A,[[x1],[x2]])

至

f = lambda t, x, A : np.dot(A, x)

此外，为确保 solve_ivp 正确解释参数 A 的形状，将其传递给 solve_ivp 和 args=(A,)（注意逗号的使用).

参数 space 中 solve_ivp 的问题

Problem with solve_ivp in parameter space

python

matlab

scipy

ode