Реализуйте метод линий для решения PDE в Python scipy с сопоставимой производительностью с ode15s Matlab - PullRequest
Новая
       104

Реализуйте метод линий для решения PDE в Python scipy с сопоставимой производительностью с ode15s Matlab

1 голос
/ 18 февраля 2020

Я хочу использовать метод линий для решения уравнения тонкой пленки . Я реализовал его (с помощью gamma=mu=0) Matlab, используя ode15s, и, кажется, он работает нормально: 

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? В одном я не уверен, правильно ли я делаю векторизацию. Также я правильно написал функцию? Я знаю, что это жесткий ODE, но сверхмалый временной шаг, сделанный

  2. Действительно ли разница сводится только к разнице в решателе од? scipy.integrate.solve_ivp(f, method='BDF') является рекомендуемой заменой ode15s в соответствии с официальным numpy веб-сайтом . Но для этого конкретного примера разница в производительности составляет одну секунду, а для решения требуются годы. Разница намного больше, чем я думал.

  3. Есть ли другие альтернативные методы, которые я могу попробовать в Python для решения подобных PDE? (что-то вроде конечной разницы / метод линий) Я имею в виду использование существующих библиотек, предпочтительно библиотек scipy.

...