参数空间中solve_ivp的问题

Posted 2023-02-18

【中文标题】参数空间中solve_ivp的问题

【英文标题】：Problem with solve_ivp in parameter space

【发布时间】：2021-12-31 22:43:19

【问题描述】：

我在创建 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：形状 (2,2) 和 (2,1,2) 未对齐：2 (dim 1) != 1 (dim 1)

或类似的，取决于我已经尝试将f 重写为的内容。据我了解，solve_ivp 期望 x0 数组中有一个值，而我返回一个 2x1 向量。它也不接受向量作为 x0 数组中的值，例如 [[x0,0]]

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

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

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

【参考方案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,)（注意逗号的使用）。

【讨论】：

好吧，我还得调整剧情。现在它起作用了（似乎）。我必须做一些测试，看看它是否也适用于其他例子。谢谢