Уравнение нелинейной диффузии Колмогорова-Петрова-Пискунова - Mathematica - Страница 2

Зосима · 18.10.2014, 15:11

Студворк — интернет-сервис помощи студентам

У меня получается 200с + 53с
последовательным интегрированием против
630с двойным.

В Matlab наблюдается такая же тенденция

@Mysterious Light · 19.10.2014, 23:16

Выкладываю сюда все свои результаты, поскольку дальше разбираться не хочется.

Интегрирование

Пусть нам дано

Haskell
1
2
3
4
5
6
7
8
n = 3;
a = 0;
b = 1;
k = 7/10;
u0 = Sin[x] + x;
\[CurlyPhi][i_] = x^i (1 - x);
\[Rho][j_] = t^j (1 - t);
Integrant = Expand[R \[CurlyPhi][i] \[Rho][j]]

Обращаем внимание на то, что Head[Integrant] === Plus. Посмотрим, как Вольфрам берёт два последовательных интеграла от каждого слагаемого в общем, численном случае:

Haskell
1
2
Monitor[Map[(ie = #; ii = Integrate[ie, {x, 0, 1}]; 
     Integrate[ii, {t, 0, 1}]) &, Integrant] // Timing, {ie, ii}]

Можно поступить так:

Haskell
1
2
3
Timing[Ix = Integrate[Integrant, {x, a, b}]] // First
Timing[Ixt = Integrate[Ix, {t, a, b}]] // First
II[i_, j_] = Ixt;

Довольно дорогая операция. У меня при n=3 для указанного u0 считает треть часа.
Составляем уравнения:

Haskell
1
2
3
eqns = Flatten@Array[II[##] == 0 &, {n, n}];
eqns = Simplify[eqns]; (* очень долго работает *)
vars = Flatten@Array[c, {n, n}];

далее можно писать Solve[eqns, vars], но я пойду другим путём, который короче.

Нахождение корней системы уравнений

Предлагаю решать через FindRoot.

Haskell
1
2
3
N[c[args___]] := c[args];
fs = N@Flatten@Array[II, {n, n}];
U = FindRoot[fs, Transpose[{vars, initPoint}]]

Проверим стахостически, что решение единственно:

Haskell
1
2
3
4
Rationalize[
  Table[Array[c, {n, n}] /. 
    FindRoot[fs, Transpose[{vars, RandomReal[{ -100, 100}, n^2]}]], {200}], 
  10^-10] // DeleteDuplicates

Визуализация

Haskell
1
2
3
U1 = First[U];
ComplexExpand[Im[U1]] (* Должен быть нуль *)
DensityPlot[Re[U1], {x, 0, 1}, {t, 0, 1}, PlotRange -> All]

PS. Вопрос знатокам WM: как переписать N[c[__]], чтобы N не переводило аргументы c из натуральных индексов в вещественные, вроде c[1., 2.]?
Я сейчас использую N[c[args___] = c[args], но поэтому N впадает в бесконечную рекурсию. Ещё можно сделать N без ничего, а потом Rationalize, чтобы обратно вернуть индексы в целочисленный вид. Но это я тоже рассматриваю как костыль.

@Mysterious Light · 20.10.2014, 00:00

Прилагаю файл.

Считает 1000c интегралы (символьно, в прошлом посте я ошибочно написал «численно»), после чего составляются уравнений за ~мин, остальное считается быстро.

Зосима · 20.10.2014, 12:10

Исправленный вариант, который что-то таки выдал

Haskell
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
a = 0;
b = 1;
k = 0.75;
 
n = 3;
ut0 = (Sin[Pi*x] + x)*(1-t) + x*t; (*выполняем граничные условия*)
Fun[u1_] = k*u1*(1 - u1); 
fi[i_] = x^i*(1 - x);
ro[i_] = t^i*(1 - t);
 
CC = Array[c, n^2];
u = ut0 + Sum[  c[(i - 1)*n + j]*fi[i]*ro[j] , {i, n}, {j, n}];
dU2dx2 = D[ u, {x, 2}];
dUdt = D[ u, t];
Rr = dU2dx2 - dUdt + Fun[u] (*Невязка*)
 
nn = 1;
Z = Array[z, n^2];
For[i = 1, i <= n, i++,
 For[j = 1, j <= n, j++,
  fx = Expand[ Rr*fi[i]*ro[j] ];
  ft = Integrate[ fx, {t, a, b}];
  z[nn] = Integrate[ ft, {x, a, b}]  ;
  Print[{i, j, nn++}]; (*отображаем текущий прогресс вычислений*)
  ]
 ]
Print[Z];
 
S = FindRoot[Z == 0, Transpose[{CC, Array[1 &, n^2]}] ]; (*Решаем систему*)
U = Simplify[u /. S] (*подставляем решение в начальное выражение*)
Plot3D[ U, {x, a, b}, {t, a, b}]

Результаты: Galerkin.pdf

Уравнение нелинейной диффузии Колмогорова-Петрова-Пискунова

И код для проверки:

Haskell
1
2
3
4
5
6
7
(*Точное решение для проверки*)
NDSolve[{D[f[x, t], t] == 
    D[f[x, t], {x, 2}] + k*f[x, t]*(1 - f[x, t]), 
   f[x, 0] == Sin[Pi*x] + x, f[0, t] == 0, f[1, t] == 1 }, 
  f, {x, 0, 1}, {t, 0, 1}];
Plot3D[Evaluate[f[x, t] /. %], {x, 0, 1}, {t, 0, 1}, PlotRange -> All,
  AxesLabel -> Automatic]

Новые блоги и статьи Все статьи Все блоги /
Про "Windows must die!" kumehtar 04.11.2026 Знаете, чет накатило. . . . Когда-то это было почти боевым кличем: мол, нормальный человек сидит на DOS, OS/ 2 или Linux, а вот это вот - Must Die. А теперь человек произносит «мастдай» с усмешкой - и. . .	Красочные поля с номером страницы svg+background-image russiannick 04.11.2026 Красочные поля с номером страницы не менее важны чем в эпоху бумажных книг. Посему спешу представить почтеннейшей публике этот важный инструмент который эти поля и формирует. В процессе разработки. . .	Вступаем в новую фазу генерации favicons, векторную. russiannick 02.11.2026 Теперь нет нужды скрывать что используется простой svg. Исходный код изображения помещаем в поле и жмем кнопку. Исходник преобразуется в вид корректный для помещения в favicon или в url(). Нажатие. . .	Генератор паролей (вариант №4) Maks 01.11.2026 После публикации 3-го варианта генератора паролей, разослал его своим коллегам (с целью оценки), в итоге были получены следующие рекомендации: 1. Вернуть возможность использования знаков \| 1 l I;. . .
Генератор паролей (вариант №3) Maks 30.10.2026 Финальная версия генератора паролей 1С. Доработка решения в этот раз коснулась только визуальной части, а именно, добавлен реквизит ДопПараметры (тип Булево), который скрывает/ отображает реквизиты:. . .	Магическое мышление и магия, современная наука не отличима от магии (говорит ИИ Джемени) zorxor 30.10.2026 Джемени, расскажи, пожалуйста, что такое магическое мышление, что такое магия, почему современная наука чем-то напоминает магию? Фраза «Любая достаточно развитая технология неотличима от магии». . .	За какое время ИИ может подобрать пароль по хешу (рассказывает Джемени ИИ) zorxor 30.10.2026 Расскажи, Джемени, если известен хэш МД5, например, за сколько ИИ сможет подобрать соответствующий пароль? Искусственный интеллект (ИИ) может подобрать простой или средний по сложности пароль по. . .	Кто убил Suno? Почему наше солнышко погасло? (диалог с ИИ Джемени) zorxor 30.10.2026 Как убивали Suno. Хронология. «Убийство» музыкальной нейросети Suno — это затяжной процесс, в котором технологический хайп стартапа разбился о жесткую реальность судебных исков от гигантов. . .