在Python scipy中实现求解偏微分方程的行方法，其性能与Matlab的ode15s相当

Question

我想用线的方法来解thin-film equation。我已经用gamma=mu=0用ode15s用Matlab实现了它，它似乎工作得很好：

N = 64;
x = linspace(-1,1,N+1);
x = x(1:end-1);
dx = x(2)-x(1);

T = 1e-2;
h0 = 1+0.1*cos(pi*x);

[t,h] = ode15s(@(t,y) thinFilmEq(t,y,dx), [0,T], h0);

function dhdt = thinFilmEq(t,h,dx)
    phi = 0;
    hxx = (circshift(h,1) - 2*h + circshift(h,-1))/dx^2;
    p = phi - hxx;
    px = (circshift(p,-1)-circshift(p,1))/dx;
    flux = (h.^3).*px/3;
    dhdt = (circshift(flux,-1) - circshift(flux,1))/dx;
end

胶片在一段时间后会变平，在很长一段时间内，胶片应趋向于h(t->inf)=1。我没有做过任何严格的检查和收敛性分析，但至少在花不到5分钟的时间编写代码后，结果看起来很有希望

我想用Python做同样的事情，我尝试了以下方法：

import numpy as np
import scipy.integrate as spi

def thin_film_eq(t,h,dx):
    print(t)   # to check the current evaluation time for debugging
    phi = 0
    hxx = (np.roll(h,1) - 2*h + np.roll(h,-1))/dx**2

    p = phi - hxx   
    px = (np.roll(p,-1) - np.roll(p,1))/dx       
    flux = h**3*px/3        
    dhdt = (np.roll(flux,-1) - np.roll(flux,1))/dx

    return dhdt

N = 64
x = np.linspace(-1,1,N+1)[:-1]
dx = x[1]-x[0]

T = 1e-2
h0 = 1 + 0.1*np.cos(np.pi*x)

sol = spi.solve_ivp(lambda t,h: thin_film_eq(t,h,dx), (0,T), h0, method='BDF', vectorized=True)

我在函数中添加了一个print语句，以便检查程序的当前进度。由于某些原因，它的时间步长非常小，在等待几分钟后，它仍然停留在t=3.465e-5，dt小于1e-10。（当我把这个问题打完的时候，我还没有写完，而且可能在任何合理的时间内都不会写完）。对于Matlab程序，它是在一秒钟内完成的，只需要14个时间步（我只指定时间跨度，它输出14个时间步，其他所有内容都保持默认值）。我想问以下问题：

1. 我是否做错了什么事情，大大降低了Python代码的计算速度？我应该为solve_ivp函数调用选择什么设置？我不确定的一件事是我是否正确地进行了矢量化。我是否以正确的方式编写了函数？我知道这是一首僵硬的颂歌，但是

2. 差异真的只是ode解算器中的差异吗^根据官方{}{a2}的说法，{}是{}的推荐替代品。但是对于这个特定的例子，性能差异是1秒，而解决这个问题需要花费很多时间。差别比我想象的要大得多

3. 我可以在Python中尝试其他方法来解决类似的PDE吗？（沿着有限差分线/线的方法）我的意思是利用现有的库，最好是scipy中的库