Явный метод Эйлера слишком медленный для реакции-диффузии

10

Я решаю систему реакции-диффузии Тьюринга с помощью следующего кода C ++. Это слишком медленно: для текстуры 128x128 пикселей допустимое количество итераций составляет 200, что приводит к задержке в 2,5 секунды. Мне нужно 400 итераций, чтобы получить интересное изображение, но 5 секунд ожидания - это слишком много. Кроме того, размер текстуры должен быть на самом деле 512x512 - но это приводит к огромному времени ожидания. Устройства iPad, iPod.

Есть ли шанс сделать это быстрее? Метод Эйлера сходится медленно (википедия) - более быстрый метод позволил бы уменьшить количество итераций?

РЕДАКТИРОВАТЬ: Как отметил Томас Климпел, строки: «if (m_An [i] [j] <0.0) {...}», «if (m_Bn [i] [j] <0.0) {...}» задерживают схождение: после удаления значимое изображение появляется после 75 итераций . Я закомментировал строки в коде ниже.

void TuringSystem::solve( int iterations, double CA, double CB ) {
    m_iterations = iterations;
    m_CA = CA;
    m_CB = CB;

    solveProcess();
}

void set_torus( int & x_plus1, int & x_minus1, int x, int size ) {
    // Wrap "edges"
    x_plus1 = x+1;
    x_minus1 = x-1;
    if( x == size - 1 ) { x_plus1 = 0; }
    if( x == 0 ) { x_minus1 = size - 1; }
}

void TuringSystem::solveProcess() {
    int n, i, j, i_add1, i_sub1, j_add1, j_sub1;
    double DiA, ReA, DiB, ReB;

    // uses Euler's method to solve the diff eqns
    for( n=0; n < m_iterations; ++n ) {
        for( i=0; i < m_height; ++i ) {
            set_torus(i_add1, i_sub1, i, m_height);

            for( j=0; j < m_width; ++j ) {
                set_torus(j_add1, j_sub1, j, m_width);

                // Component A
                DiA = m_CA * ( m_Ao[i_add1][j] - 2.0 * m_Ao[i][j] + m_Ao[i_sub1][j]   +   m_Ao[i][j_add1] - 2.0 * m_Ao[i][j] + m_Ao[i][j_sub1] );
                ReA = m_Ao[i][j] * m_Bo[i][j] - m_Ao[i][j] - 12.0;
                m_An[i][j] = m_Ao[i][j] + 0.01 * (ReA + DiA);
                // if( m_An[i][j] < 0.0 ) { m_An[i][j] = 0.0; }

                // Component B
                DiB = m_CB * ( m_Bo[i_add1][j] - 2.0 * m_Bo[i][j] + m_Bo[i_sub1][j]   +   m_Bo[i][j_add1] - 2.0 * m_Bo[i][j] + m_Bo[i][j_sub1] );
                ReB = 16.0 - m_Ao[i][j] * m_Bo[i][j];
                m_Bn[i][j] = m_Bo[i][j] + 0.01 * (ReB + DiB);
                // if( m_Bn[i][j] < 0.0 ) { m_Bn[i][j]=0.0; }
            }
        }

        // Swap Ao for An, Bo for Bn
        swapBuffers();
    }
}
AllCoder
источник
Кроме того, я хочу отметить, что предпочтительно, чтобы вы не пересекались с вопросами, поскольку, похоже, вы задавали очень похожие вопросы и здесь, и здесь .
Годрик Провидец
Вы случайно не видели работы Грега Турка по этому поводу ?
JM
@JM: пока нет. Я только что попытался запустить его код: для этого требуется X-сервер с PseudoColor, то есть 8-битная глубина цвета. Я думаю, что я не могу предоставить это на OSX. Я пробовал различные VNC-серверы, но не повезло.
AllCoder
Я думаю, что вы все еще должны быть в состоянии адаптировать подход Турка к данному вопросу; Реакции диффузии реакции, по-видимому, в наше время достаточно широко используются в компьютерной графике.
JM
1
Я могу ошибаться, но часть с m_An [i] [j] = 0.0; может фактически добавить элемент к этой системе, который не может быть смоделирован дифференциальным уравнением с непрерывной правой частью. Это немного затрудняет поиск более быстрого решателя.
Томас Климпел

Ответы:

9

Вы, кажется, ограничены стабильностью, которая ожидается, так как диффузия жесткая, когда вы улучшаете сетку. Хорошие методы для жестких систем по крайней мере частично неявны. Это займет некоторое усилие, но вы можете реализовать простой многосеточный алгоритм (или использовать библиотеку) для решения этой системы со стоимостью менее десяти «рабочих единиц» (по сути, стоимость одного из ваших временных шагов). Когда вы уточняете сетку, количество итераций не будет увеличиваться.

Джед браун
источник
Если бы здесь была только диффузия, он мог бы использовать метод ADI, такой как Дуглас-Ганн, и все было бы хорошо. Однако, по моему собственному опыту, часть реакции часто намного хуже в отношении жесткости, в дополнение к тому, что она сильно нелинейная.
Томас Климпел
1
У ADI, к сожалению, ужасная память. Также обратите внимание, что реакция может быть обработана неявно независимо от того, является ли диффузия. При уточнении сетки диффузия в конечном итоге станет доминирующей, но мы не можем сказать, где находится порог, не зная постоянных.
Джед Браун
Пример кода, реализующего обратный Эйлер для этого (в Python), можно найти
Дэвид Кетчесон,
@DavidKetcheson: использование неявных методов требует решения уравнения? Вот почему в коде есть linalg.spsolve ()?
AllCoder
1
@AllCoder Да, это требует решения, но решение может быть выполнено намного быстрее, чем все временные шаги, необходимые для стабильности явного метода.
Джед Браун
2

С практической точки зрения: процессор A5 не такой мощный, поэтому вы можете подождать несколько итераций HW или, если ваш ipod / ipad будет подключен к Интернету, решить вашу проблему удаленно или в облаке.

fcruz
источник
Я удивлен, как малая мощность предлагает A5. Как Pages, Safari и другие большие приложения могут работать так хорошо? Мне нужно генерировать случайные, абстрактные изображения, думал, что морфогенез будет достаточно простым ..
AllCoder
Что ж, A5 - это энергоэффективный процессор, оптимизированный для Интернета и видео (Pages, Safari и т. Д.). Напротив, большинство числовых рабочих нагрузок выполняет тонны операций с плавающей запятой и перемещений данных, эти функции не являются целью мобильных процессоров с низким энергопотреблением.
fcruz
0

Эйлер медленно сходится по отношению к другим методам, однако я не думаю, что это то, что вас интересует. Если вы просто ищете «интересные» изображения, увеличьте размер своего временного шага и выполняйте меньше итераций. Проблема, как указывает Джед, состоит в том, что явный метод Эйлера имеет проблемы со стабильностью с большими временными шагами по отношению к размеру сетки. чем меньше ваша сетка (т. е. чем выше разрешение вашего изображения), тем меньше должен быть ваш временной шаг для ее учета.

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

Годрик Провидец
источник
Эта проблема ограничена стабильностью, поэтому простое увеличение временного шага не будет работать.
Джед Браун
Если я изменяю 0,01, например, 0,015, то я получаю «концентрацию хим. Зр. Около нуля» во всех точках - т.е. серый квадрат. Вот происхождение моего кода: drdobbs.com/article/print?articleId=184410024
AllCoder
Да, это будет результатом проблем со стабильностью, упомянутых Джедом. Как он упоминает в своем ответе, использование неявного метода, характеризующегося лучшей стабильностью работы, решит эту проблему для вас. Я обновлю свой ответ, чтобы удалить ненужную информацию.
Годрик Провидец