Реализация численного решения модели Фитцхью-Нагумо с коэффициентом пространственной диффузии

Я пытаюсь придумать реализацию Python для модели Фитцхью-Нагумо.

V_t = V_xx + V(V - a)(1 - V) - W + I
W_t = eps(beta*V - W)

Используя действительно простой код для eps = 0.05, a = 0.2, beta = 5, I = .1, я могу численно решить систему (без V_xx), но я не могу понять, как реализовать пространственное рассеивание.

def func_v(v, w):
    return v * (1 - v) * (v - .2) - w + .1


def func_w(v, w):
    return .05 * (5 * v - w)


def get_yn(t0, v, w, h, t):
    while t0 < t:
        w += h * func_w(v, w)
        v += h * func_v(v, w)
        t0 += h
    return v, w

Я знаю, что формула центрированной разности для производных второго порядка

V_xx(x_i, t) = (V(x_i+1, t) - 2*V(x_i, t) + V(x_i-1, t)) / dx^2

но как мне реализовать разные значения для x_i (скажем, от x=0 до 10), чтобы волна распространялась вдоль оси x?

введите здесь описание изображения

В результате должна получиться волна, которая распространяется примерно так.


python differential-equations pde
person vlovero    schedule 16.12.2019    source источник
comment
Есть ли у вас какие-либо тестовые данные или тестовые образцы, которые мы можем использовать?   -  person wundermahn    schedule 16.12.2019
comment
@wundermahn, что ты имеешь в виду под тестовыми данными? Например, временной интервал и начальные условия, которые я использую для решения без V_xx?   -  person vlovero    schedule 16.12.2019
comment
Что вы думаете о предыдущих вопросах по теме, например stackoverflow.com/q/14915398/3088138 или stackoverflow.com/q/38116507/3088138   -  person Lutz Lehmann    schedule 17.12.2019
comment
@Dr.LutzLehmann Я видел вопрос в MATLAB, и когда я перевел код на python, у меня это не сработало.   -  person vlovero    schedule 17.12.2019

Ответы (1)


arrow_upward
1
arrow_downward

Решатель ОДУ (на самом деле любая компьютерная программа) может обрабатывать только задачи, имеющие конечномерное состояние. Состояние в этом PDE представляет собой пару функций v,w из x. В необходимой общности они не могут быть представлены в компьютере. Таким образом, вам нужно работать с конечными приближениями. Первой, которая считается достаточной во многих контекстах, является простая таблица функций. Затем производные x вычисляются с использованием формул конечных разностей.

x = np.linspace(0,L,N+1);
dx = x[1]-x[0];
v0,w0 = initial_functions(x);

def func_v(v, w):
    d2v = -2*v;
    d2v[0] += v[-1]+v[1];
    d2v[1:-1] += v[:-2] + v[2:]
    d2v[-1] += v[-2]+v[0];
    return d2v/dx**2 + v * (1 - v) * (v - .2) - w + .1

и т. д.

Для проверки концепции может быть достаточно метода Эйлера, но полученные значения будут сомнительными. Используйте метод более высокого порядка, чтобы получить полезные результаты, не используя смехотворно маленькие временные шаги.

person Lutz Lehmann    schedule 17.12.2019
comment
Чтобы использовать новый func_v с начальными v0 и w0, должен ли я повторить итерацию, как v0 += func_v(v0, w0) ; w0 += func_w(v0, w0)? Потому что, когда я его использовал, новый график просто перемещался вверх по оси v, а не по оси x. - person vlovero; 17.12.2019
comment
Вы можете вызвать свой неизмененный метод Эйлера vn, wn = get_yn(t0, v0, w0, h, t). Затем используйте plot(x,vn,x,wn). То, что вы получите, зависит от начального условия (например, петля вокруг круга радиуса 4, x=linspace(0,2*pi,N+1)[:-1]; v0=4*np.cos(x); w0=4*sin(x) или подобного, и, что наиболее важно, от размера шага h. Вам потребуется некоторое время стабилизации, когда решение сходится к предельному циклу уравнение Вандерполя. - person Lutz Lehmann; 17.12.2019
comment
Я не уверен, почему это все еще не работает. Может быть, это потому, что мой get_yn изначально должен был принимать числа с плавающей запятой вместо массивов? Сейчас происходит то, что я получаю что-то похожее на то, что хочу, но оно смещено влево, а не вправо. - person vlovero; 17.12.2019
comment
Что значит не работает? Если я правильно понял, код запускается, но дает неожиданные результаты? Обратите внимание, что ОДУ для решения бегущей волны имеет порядок 2+1, так что его система первого порядка имеет размерность 3. GIF-анимация показывает гомоклиническую орбиту в этом трехмерном пространстве, что является совершенно особой ситуацией, я не уверен если оно устойчиво при возмущении. Определить положение стационарной точки несложно, но все, что после этого, — математическая задача, здесь не по теме. - person Lutz Lehmann; 19.12.2019
comment
графики рассеиваются, но вместо того, чтобы двигаться по оси x, они ведут себя более похоже на то, что вы получили бы из уравнения головы вместо бегущей волны. - person vlovero; 19.12.2019
comment
Да, это то, что я тоже вижу, используя odeint. Вам понадобится тщательно сконструированное начальное условие и правильные параметры. Также проверьте коэффициент рассеяния, слишком сильное перемешивание приведет к выравниванию решения по пространству. - person Lutz Lehmann; 19.12.2019