使用scipy.solve_ivp复现Matlab ode45求解线性中心常微分方程报错问题

问题背景 #

尝试将Matlab中的常微分方程示例迁移到Python中使用scipy.integrate.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

Matlab实现效果 #

重写的Python实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

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数组中传入单个值，但自定义f函数返回的是2x1向量，同时solve_ivp也不接受类似[[x0,0]]这样的向量作为x0数组的输入。
需要确认scipy.solve_ivp是否可以像Matlab的ode45一样实现参数空间的求解，如果可以的话应该如何实现，还是需要采用其他计算方式？
补充说明：已检查过其他所有矩阵和返回值都和Matlab计算结果一致。

解决方法 #

solve_ivp完全可以实现和Matlab ode45一致的求解能力，报错均为用法不匹配导致，核心修改点如下：

• ODE函数返回格式修正：原f函数返回的是shape为(2,1)的二维数组，solve_ivp要求ODE函数返回和输入状态x维度一致的一维数组，修改f为f = lambda t,x,A : np.dot(A, x)即可匹配维度要求，x本身是长度为2的一维数组，点乘后返回长度为2的一维数组符合接口要求
• 求解时间区间修正：原代码中solve_ivp传入的求解区间是[xmin, xmax]，对应状态变量的范围而非时间范围，需修改为和Matlab一致的时间区间[0, 100]
• t_eval参数修正：原代码传入的yval是空间坐标采样值，实际应该传入时间采样点，比如t_eval = np.linspace(0, 100, 500)来保证轨迹绘制的光滑度
• 逻辑补全：原Python循环仅实现了初始条件[x0, 0]的求解，还需补充初始条件[0, x0]的求解和绘图，和Matlab逻辑对齐
• 经最终验证，绘图时x轴取solution.y[1]、y轴取solution.y[0]可得到和Matlab一致的效果

内容的提问来源于stack exchange，提问作者phx16