Skip to content
Get App

Лукьяненко Д. В. - Численные методы - 27. Уравнения в ЧП: многомерные уравнения

Лекция о численных методах решения многомерных уравнений в частных производных с примерами и постановкой задачи.

Key Takeaways

  • Многомерные уравнения в частных производных требуют расширения постановки задачи с дополнительными пространственными координатами.
  • Для решения нелинейных уравнений аналитические методы часто невозможны, поэтому применяются численные схемы.
  • Граничные и начальные условия должны быть согласованы, чтобы избежать переопределённости задачи.
  • Метод свода к системе ОДУ и схемы Розенброка обеспечивают более эффективное решение по сравнению с простыми явными схемами.
  • Построение сеток и аппроксимация производных играют ключевую роль в точности и устойчивости численных методов.

What the video covers

  • Обсуждение базовых особенностей многомерных уравнений в частных производных с пространственной многомерностью.
  • Переход от одномерных к двумерным уравнениям по пространству с добавлением координаты Y.
  • Постановка задачи на примере уравнения типа Бюргерса с нелинейными слагаемыми.
  • Геометрическая интерпретация задачи как трёхмерного параллелепипеда в пространстве X, Y, T.
  • Необходимость граничных и начальных условий для корректной постановки задачи.
  • Обсуждение явных численных схем и метода свода к системе ОДУ с применением схем Розенброка.
  • Пояснения по построению сеток и аппроксимаций производных в двумерном случае.
  • Важность устойчивости численных схем и практические рекомендации по выбору параметров.
  • Примеры реализации численных методов и обсуждение их сложности и громоздкости.
  • Общая концепция курса — изучение численных методов на конкретных примерах для практического понимания.

Answers

Questions about this video

Что понимается под многомерностью в контексте уравнений в частных производных?

Под многомерностью понимается многомерность по пространству, то есть наличие нескольких пространственных координат, например X и Y, в отличие от одномерных задач, где учитывается только одна пространственная координата.

Почему для решения уравнения типа Бюргерса применяются численные методы?

Уравнение типа Бюргерса содержит нелинейные слагаемые, такие как U в кубе, что делает невозможным получение аналитического решения, поэтому применяются численные методы для приближённого решения.

Какие условия необходимы для корректной постановки задачи в многомерном случае?

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

Full Transcript — Download SRT & Markdown

00:19
Speaker A
Коллеги, сегодня мы с вами поговорим про некоторые базовые особенности решения многомерных уравнений в частных производных. И под многомерностью мы будем подразумевать многомерность по пространству, потому что формально до этого мы тоже рассматривали двухмерные уравнения в частных производных. Одно измерение — это X, другое измерение — это время. Но мы всё равно называли такие задачи одномерными, подразумевая, что по пространству задача была поставлена только по иксу, то есть по одному измерению. Сейчас же мы рассмотрим пример, где по пространству будет присутствовать и координата по иксу, и координата по игреку. И в
00:43
Speaker A
качестве примера давайте рассмотрим следующий пример. Вспомним девятнадцатую лекцию, где мы рассматривали одномерное по пространству уравнение типа Бюргерса или уравнение параболического типа.
01:06
Speaker A
Давайте я сейчас выпишу соответствующую постановку, напомню эту постановку из девятнадцатой лекции и обобщу её на двухмерный по пространству случай. Итак, на девятнадцатой лекции мы рассматривали следующую задачу для уравнения типа Бюргерса.
01:21
Speaker A
Эпсилон D²U по DX² минус DU по DT. Сразу обращаю ваше внимание, что я оставляю место на слагаемые, которые я попозже допишу.
01:43
Speaker A
U dU по dx и - U в кубе. Это уравнение определяло неизвестную функцию в следующей области: X принадлежит от A до B.
02:00
Speaker A
И T принадлежит от T нулевого до T большого. Дифференциальное уравнение параболического типа. Для того, чтобы поставить задачу корректно, чтобы она имела одно единственное решение, мы доопределяли неизвестную функцию двумя дополнительными условиями на границах отрезка от A до B и в начальный момент времени.
02:19
Speaker A
Поэтому у нас было следующее: U от A T — тоже место оставляю про запас для дополнительной переменной. Это у нас было U left — условие Дирихле на левой границе.
02:54
Speaker A
На правой границе у нас функция была доопределена правым граничным условием. И всё это было, опять же, при T от T ноль до T большого. И ещё было начальное условие: U от X T ноль — это было U initial от
03:13
Speaker A
X. X принадлежала от A до B. Есть вот такую вот задачу мы с вами рассматривали на девятнадцатой лекции. И напомню вам, на девятнадцатой лекции я приводил вот такой вот рисунок.
03:44
Speaker A
То есть наша цель заключалась в чём? Восстановить, то есть найти функцию U двух переменных X и T вот в этой прямоугольной области.
04:04
Speaker A
Теперь же, теперь мы рассматриваем обобщение этого уравнения на двухмерный по пространству случай. В этом случае уравнение усложняется следующим образом.
04:18
Speaker A
Добавляется вторая производная по игреку. Естественно, теперь предполагается, что функция U — это не функция двух переменных X и T, а функция уже трёх переменных X, Y и T.
04:31
Speaker A
Вот здесь добавляется частная производная по игреку. Соответственно, к чему это приводит? А это приводит к тому, что у нас появляется ещё одно дополнительное измерение.
04:49
Speaker A
И, если мы конкретизируем, что по игреку задача ставится на отрезке от C до D, то как раз и получается, что уже вот это уравнение определяет неизвестную функцию на вот таком вот открытом множестве.
05:07
Speaker A
То есть фактически, что это из себя представляет? Как раз прямоугольничек от A до B по иксу, от C до D по игреку, открытый, не включая свои границы.
05:27
Speaker A
Соответственно, что у нас ещё вот здесь меняется? Теперь понятно, почему я оставил вот здесь место, здесь добавляю Y, Y, Y и вот здесь дописываю, что Y принадлежит от C до D.
05:38
Speaker A
Смотрите, какая у нас геометрическая интерпретация. Геометрическая интерпретация теперь гораздо более сложная. Естественно, мы решаем задачу уже не в каком-то прямоугольнике, а уже в трёхмерном параллелепипеде в пространстве. Поэтому я вот здесь нарисовал вот такую картинку, чтобы было понятно, по какому принципу строится
06:01
Speaker A
соответствующий параллелепипед. Естественно, чтобы на картинке не было вот такой мешанины, часть проекций я убираю и убираю, в принципе, всё лишнее.
06:27
Speaker A
То есть вот наша трёхмерная область, наш трёхмерный параллелепипед, в котором мы решаем нашу задачу. Но понятно, предыдущие рисунки нужны были для чего?
06:38
Speaker A
Чтобы понимать, по какому принципу был построен вот этот параллелепипед. Итак, и также, как в девятнадцатой лекции, для постановки задачи нам нужны граничные и начальные условия.
06:51
Speaker A
И вот здесь как мы действуем? На нижней границе по времени у нас условие ставится в начальный момент времени. Но теперь ставится не на отрезке по иксу, а уже на соответствующей грани параллелепипеда. Поэтому я вот здесь дописываю вот такую конструкцию и тоже
07:08
Speaker A
вот здесь дописываю, что теперь не X принадлежит AB, а это множество точек, принадлежащих соответствующему прямоугольнику, который совпадает с нижней гранью.
07:30
Speaker A
Нам нужны граничные условия. Вот, видите, левое граничное условие, правое граничное условие тоже поставлены на соответствующих гранях трёхмерного параллелепипеда. Поэтому вот здесь и добавилась ещё одна переменная Y, и здесь, кроме T, что T принадлежит от T0 до T большого, мы ещё и дописали, что Y
07:46
Speaker A
принадлежит от C до D, чтобы как раз вот эти граничные условия определяли нашу функцию именно на соответствующей грани параллелепипеда. И ещё вы догадываетесь, что нужно дописать? Нужно граничное условие на грани, вот задней, ну, если мы так на рисунок смотрим. Это что у нас будет?
08:11
Speaker A
Это будет U от X C T. Я назову это U bottom. А на передней грани будет U top.
08:35
Speaker A
Здесь соответственно будет X T. Почему вот это нижняя, нижней гранью я считаю? Ну, понятно, относительно направления оси Y. То есть относительно направления оси Y вот это условие является нижним. А вот здесь у нас U от X D T, это U
08:49
Speaker A
top от X T. Здесь X принадлежит от A до B. T по-прежнему от T0 до T большого. Коллеги, на всякий случай напомню вам, что вот здесь у нас есть произвол некоторый, где мы точки вкалываем, а где выкалываем.
09:15
Speaker A
С чем это связано? То есть можете, например, вот напоминаю, например, почему вот здесь точка T0 выколота? А потому что при T0 функция доопределена начальным условием. И аналогично вы можете, например, какой-нибудь источник здесь вколоть, но нужно будет где-нибудь вот здесь выколоть. То есть нужно
09:40
Speaker A
следить просто затем, чтобы у нас соответствующие множества не пересекались, потому что если будет пересечение, это будет эквивалентно тому, что у нас задача переопределённая, что в какой-то точке поставлено не одно условие на функцию, а два. Зачем же нам решать переопределённую задачу, лишние
09:59
Speaker A
заморочки себе придумывать? Поэтому вот, вот постановка уравнения в двухмерном по пространству случае. Опять же хочу акцентировать внимание на терминологию. Иногда могут называть эту постановку и трёхмерной, ну потому что формально у нас функция U зависит от трёх переменных X, Y, T. Либо называют
10:18
Speaker A
2+1-мерной, подразумевая под двойкой двухмерную по пространству, а под единичкой, что ещё есть дополнительное измерение — это время. Так, теперь какая наша цель? Наша цель, естественно, с вами, как и была, найти функцию U от X, Y, T в соответствующем
10:43
Speaker A
трёхмерном параллелепипеде. Естественно, если мы можем найти её аналитически, мы с вами никакие численные методы не применяем. Мы, естественно, исходим из того, что аналитически решение не находится. А я специально как раз и выписал достаточно сильно нелинейное уравнение. Во-первых, оно нелинейное за счёт вот этого
11:11
Speaker A
слагаемого и вот ещё дополнительная сильная нелинейность — U в кубе. То есть аналитическое решение вот этого уравнения найти нельзя. Поэтому ничего нам не остаётся, как применять численные методы. И давайте для однообразия, структуру построения ну и рассмотрения соответствующих численных схем, мы
11:33
Speaker A
построим в том же стиле, что было и на девятнадцатой лекции. А на девятнадцатой лекции мы с вами что делали? Сначала мы рассмотрели простейший подход, основанный на построении численной схемы на основе чисто явной схемы, получаемой с помощью простейшего чисто явного
11:54
Speaker A
шаблона. То есть ссылались на курс теория разностных схем и оттуда применяли соответствующие знания. И второй подход мы рассматривали, это сводили задачу для уравнения в частных производных с помощью метода прямых к большой системе обыкновенных дифференциальных уравнений, которую уже решали с помощью одной из
12:17
Speaker A
схем из одностадийного семейства схем Розенбро.
12:42
Speaker A
схем из одностадийного семейства схем Розенброка. Мы с вами уже много раз говорили, что второй подход гораздо более эффективный. Но первый подход, когда мы строим простейшую схему чисто явную на основе какого-то шаблона, он очень прост для реализации. Поэтому я всё
13:09
Speaker A
равно повторю вот этот первый подход сначала, потому что это будет всё достаточно быстро. А уже потом вы увидите на втором подходе, насколько всё усложняется, насколько все выкладки становятся громоздкими. А это приводит соответственно к чему? К тому, что иногда получается, что вы потратите
13:30
Speaker A
огромное количество времени на расписывание схемы. Поэтому возникает вопрос: а может использовать менее надёжную схему, например, чисто явную, но просто взять шашки расчёта гораздо меньше. Поэтому вот именно в связи с этим, так как на ещё раз напоминаю, вот этот курс рассчитан в первую очередь на
13:49
Speaker A
студентов второго курса. И если вы в своей научной работе достаточно встречаетесь с необходимостью прикинуть вообще, как выглядит решение какого-то сложного уравнения, ну, задачи для какого-то сложного нелинейного уравнения, можно иногда, да, реализовать явную схему, примерно посмотреть, как всё выглядит.
14:14
Speaker A
Если вас результат устраивает, уже принимать решение, использовать ли более надёжные схемы, чтобы получить более детальное решение.
14:23
Speaker A
Итак коллеги напоминаю в чём заключается первый подход. А первый подход заключается в следующем. Мы естественно функцию аналитически найти не можем, поэтому что мы с вами делаем? Мы вводим на отрезке AB по иксу, на отрезке CD по игреку и на отрезке
14:49
Speaker A
от Т0 до Т большого по Т сетки. Для простоты введём равномерные сетки. Вот давайте запишу. Например, определю шажок по иксу, это b - a делить на число интервалов nx по иксу, и узлы сетки, опять же, рассматриваем простейший случай, равномерную
15:17
Speaker A
сетку, это будет, как с кучу раз записывали, a + h x умножить на n, n от нуля до n большое с индексом x.
15:30
Speaker A
Сеточку по y вводим. Шажок равномерной сетки вычисляется как d - c делённое на n y.
15:41
Speaker A
y n-тое определяется как c + h y умножить на n, n от нуля до n y. Поэтому, смотрите, сразу замечаем, у нас здесь n и здесь n. Будет у нас некоторая путаница. Поэтому, чтобы соответствующей путаницы не было, давайте по иксу будем использовать
16:11
Speaker A
индекс i, x i-тое, а по игреку будем использовать индекс j. А уже много раз используемый индекс m маленькое оставим для сетки по времени. Поэтому, давайте формально запишу и сетку по времени.
16:40
Speaker A
Определяем величину шажка тау, как делали кучу раз, делить на число интервалов сетки по времени и t m-тое определяется как t ноль + тау на m. m опять же от нуля до m большое. Вот мы ввели сетки и к чему это приводит?
17:08
Speaker A
Смотрите, напоминаю, что мы делали на девятнадцатой лекции. На девятнадцатой лекции мы говорили, когда мы ввели по иксу и по t сетки говорили "Проведём-ка через каждый узел прямую, параллельную соответствующей оси координат". Здесь же мы что делаем?
17:29
Speaker A
Здесь мы делаем следующее. Мы уже через каждый узел проводим не прямую, а плоскость, параллельную соответствующей координатной плоскости. То есть через все узлы сетки по иксу мы проводим плоскости, параллельные координатной плоскости yt.
17:50
Speaker A
Через каждый узел сетки по игреку проводим плоскости, параллельные координатной плоскости xt. И пересечение вот этих всех плоскостей как раз дадут в пространстве сетку, в узлах которой мы и будем искать соответствующие сеточные значения. Поэтому мы получим сетку, и наша цель в узлах этой сетки уже будет
18:23
Speaker A
найти соответствующие сеточные значения, которые для простоты будем определять вот так: u и, запятая, j.
18:33
Speaker A
То есть индексы по координатам будем писать внизу. Упорядочены сначала индекс по иксу, потом и индекс по игреку, и индекс по времени будем писать наверху. Тоже на всякий случай вам напоминаю.
18:50
Speaker A
Что? Естественно, мы с вами договорились о чём? Что когда мы пишем u от x, y, t, это точное значение, которое мы не знаем. Вот это мы будем находить приближённые значения. Поэтому мы держим в уме, что вот это некоторое, напоминаю,
19:08
Speaker A
условность. То есть если мы считаем u u это точная функция, то если мы её вычислим в узлах соответствующей пространственной сетки, а вот здесь забыл индекс M ещё поставить, то мы с вами что получим?
19:26
Speaker A
Мы получим точное значение функции в узле этой сетки, а вот это приближённое, поэтому я напоминаю, что вот это всего лишь обозначение, а вот это равенство неточное. И смотрите, что мы делаем с вами дальше?
19:40
Speaker A
Естественно, как я сказал, если через все узлы провели соответствующие плоскости параллельные соответствующим координатным плоскостям, здесь на графике получится жуткая мешанина. Поэтому я делаю что? Поэтому давайте, что я вот здесь отметил? Я здесь отметил линии пересечения плоскостей, проходящих через узлы по
20:05
Speaker A
иксу, через узлы по игреку с плоскостью координатной X Y. Но меня, естественно, узлы вот этой сетки не интересуют, потому что мне нужно получить соответствующую сеточку для T0, T1, T2 и так далее T большое. Поэтому, чтобы было понятно, как определяется
20:31
Speaker A
соответствующая сеточка, смотрите, что я делаю. Я на вот это расстояние вот эту конструкцию поднимаю.
20:38
Speaker A
Поднял. Получается у меня, конечно, мешанина. Нижнюю плоскость убираю. И вот получается соответствующая сеточка, соответствующая моменту времени T0.
20:56
Speaker A
И, соответственно, моя цель уже будет определить неизвестную функцию во всех узлах вот этой сетки для каждого временного слоя. Понятно? Напоминаю, кстати, терминологию, что набор сеточных значений для фиксированного момента времени называется слоем или временным слоем.
21:23
Speaker A
Соответственно, вот здесь чтобы для полноты картины, естественно, мне что нужно сделать, чтобы вот всю сетку записать, всю сетку? Мне нужно, естественно, вот такие вот конструкции для каждого временного слоя нарисовать. Но будет мешанина, поэтому я для наглядности останавливаюсь только на
21:45
Speaker A
нулевом временном слое, соответствующему моменту времени Т с индексом ноль, остальное не рисую, исходя из предположения, что всем понятно, как всё получилось. Поэтому давайте тоже немножко формализм какой.
21:58
Speaker A
Мне нужно определить вот эти все сеточные значения. При каких и, ж и м? Формально я сначала пишу, что это при и от нуля до nx, при ж от нуля до ny, при м от нуля до м большого.
22:18
Speaker A
Но дальше я замечаю следующее. Смотрите. В момент времени t0 вообще-то я все эти сеточные значения уже знаю, потому что они определяются начальным условием. Поэтому тоже напоминаю вот эти обозначения, что красненьким я рисую либо узлы сетки крестиками, либо красными кружками
22:43
Speaker A
найденные в этих узлах уже приближённые значения, но приближённые. Зелёненьким я рисую то, что мне известно точно.
22:53
Speaker A
А с учётом вот этого начального условия на нулевом временном слое я могу вычислить все сеточные значения точно, используя вот это начальное условие. Это что будет? U с индексом и ж верхний индекс, соответствующий моменту времени с индексом m, это ноль,
23:15
Speaker A
определяется как U initial от x итого y житого. И вот это всё я в программке проделаю в цикле для и соответственно, у меня здесь x от a до b включая, поэтому это всё будет для и от нуля до nx, включая граничные индексы. И ещё
23:43
Speaker A
для всех ж от нуля до ny. есть я сразу в начале программы смогу определить все сеточные значения на нулевом временном слое.
24:02
Speaker A
В принципе, сделать это смогу точно. Всё, теперь исхожу из того, что мне сеточные значения в нулевой момент времени на нулевом временном слое мне известны.
24:15
Speaker A
Теперь моя цель какая? Смотрите, что произошло на картинке. Вот даже давайте мигну. Предыдущая картинка, новая, предыдущая, новая. Что? Я здесь сдвинулся на шаг тау по оси t. То есть до этого у меня все сеточные значения были заданы в момент
24:35
Speaker A
времени t0, теперь я вот отметил сеточку соответствующую временному слою для t1. И на что мы обращаем ваше внимание? На то, что моя цель найти приближённые значения вот в этих узлах, кроме граничных. Почему кроме граничных?
24:57
Speaker A
А потому что на границе мне функция известна из граничных условий. Поэтому я сразу тоже в начале программки что могу сделать? Я могу сделать следующее. Вот смотрите, вот это что у меня такое?
25:13
Speaker A
У Ну давайте я буду конкретно для первого временного слоя писать. А по иксу соответствует индексу ноль.
25:26
Speaker A
Здесь жи. Это будет у left от игрек житого Т один. Это всё будет при жи. Смотрите, а вот здесь жи будет не от нуля до NY, а от единицы до NY - 1. Почему? Потому что вот здесь, как мы с вами написали, игрек
25:53
Speaker A
принадлежит от С до Д, не включая. Не включаются граничные точки этого отрезка. Значит, не включаются граничные узлы ноль и N с индексом игрек.
26:05
Speaker A
Дальше. На правой границе. Правая граница здесь соответствует индекс NX. Это у right от игрек житого Т один. То есть я сейчас всё пишу, как вычислить именно, вот обращаю ваше внимание, именно граничные значения в соответствующих узлах и именно для
26:33
Speaker A
первого временного слоя. Поэтому я здесь пишу один и один, ну конкретизирую, жи в таких-то пределах, от единицы до NY - 1. И давайте зададим вот на этих границах, на этой границе и на вот этой границе. Как это выписываться?
26:54
Speaker A
Вот отсюда. Это будет U, опять же в момент времени один, здесь будет и, ну, давайте для границы, которая ближе всего к нам. У неё индекс n y, это будет U топ от x и t t один. Здесь и
27:21
Speaker A
здесь и будут от нуля, включая до nx включая, потому что здесь x принадлежит от ab, от a до b включая соответствующие граничные точки, значит, включая соответствующие граничные узлы.
27:36
Speaker A
U и ноль, один, это будет U ботом x и t t один и равно от нуля до nx. То есть вот эти действия, обращаю ваше внимание, я специально использую вот эти квадратные скобки чтобы обратить внимание на то, что вот эти действия
28:09
Speaker A
делаются для вот таких вот индексов. Ну, теперь, да, я, конечно, вычислил на вот этих границах для соответствующих узлов точные значения для первого временного слоя. Ну, почему бы мне в программке это сразу не сделать и для всех остальных временных слоёв
28:33
Speaker A
для m от единицы, то есть нолик не включая, до m большого, потому что t вколотое, t выколотое, поэтому нолик, левый граничный узел по времени не включается, правый граничный узел по времени включается. То есть таким образом, я, в принципе, вот вот таким вот действием
28:56
Speaker A
задаю сеточные значения на нижней границе и всех боковых границах. Теперь остаётся вопрос следующий. А как теперь найти приближённые значения вот в этих узлах, отмеченных красненьких?
29:17
Speaker A
И напоминаю, как мы с вами действовали на девятнадцатой лекции. Мы с вами вводили понятие шаблона, а шаблон определяет фактически конфигурацию узлов сетки, которую мы используем для того, чтобы на выписать конечно-разностную аппроксимацию вот этого уравнения. И соответствующая конечно-разностная аппроксимация связывает между собой узлы
29:46
Speaker A
соответствующего шаблона. И мы с вами на девятнадцатой лекции использовали вот такой вот чисто явный шаблон, очень простой для реализации. То есть мы с вами писали, что вот здесь, давайте, индексы у нас был индекс n и индекс m.
30:22
Speaker A
Теперь у нас два индекса. Здесь будет i + 1 j и - 1 j, третий индекс это m.
30:37
Speaker A
И вот здесь у нас даже, давайте, индексы я буду писать абсолютно так же, как и вот здесь. Индексы, давайте, наверху буду писать, иначе три индекса достаточно плотно выписаны.
30:49
Speaker A
Вот это предполагается текущий временной слой, а это следующий временной слой. Кстати заметьте вот здесь я не сделал, вот когда я записал m + 1, я обратил внимание, что я что не уточнил? Смотрите, вот это сейчас, вот эта вся конструкция была
31:03
Speaker A
записана только для чего? Для m равного единице. Поэтому для произвольного m как обобщается? Вот здесь единичку, которая соответствовала первому временному слою, меняю на m. Вот теперь действительно я на всех временных слоях вот это всё вычисляю.
31:25
Speaker A
Смотрите, так вот, возвращаемся к шаблончику. У нас был вот такой шаблон. То есть фактически вто- второго индекса у нас j не было. Теперь аналогичный для двумерного случая шаблон будет вот такой вот.
31:40
Speaker A
Соответственно, здесь индексация будет, ну, если оси направлены так же, как на картинке, вот здесь будет i j + 1 m, здесь i j - 1 m.
31:55
Speaker A
И что этот шаблон позволяет сделать? Смотрите, возвращаемся к нулевому временному слою. Если мы его прикладываем к какой-то к какой-то, мы говорили, ячейке, ну, здесь под ячейкой будет вот такая конфигурация по xy трансуваться.
32:17
Speaker A
Ну, либо к узлу прикладываем, то он накладывает, когда мы распишем конечно-разностную аппроксимацию соответствующего уравнения, он накладывает связь между сеточными значениями вот в этих сколько? Шести узлах.
32:31
Speaker A
Пять узлов лежат на текущем временном слое, они известны, один узел лежит на следующем временном слое.
32:40
Speaker A
Понятно, в чём преимущество здесь как раз видно чисто явной схемы. Этот шаблон наложит связь между шестью сеточными значениями нашей неизвестной функции, при этом пять из которых известна, одна неизвестна. Понятно, что эту неизвестную будет легко выразить.
33:02
Speaker A
Поэтому давайте вот здесь запишем конечно-разностную аппроксимацию соответствующего уравнения. Это что будет? Это будет эпсилон d2u по dx в квадрате.
33:19
Speaker A
Это что будет? Это будет с учётом вот этой конфигурации u i + 1 j в m-тый момент времени - 2 u i j-тое + u i 1 j делить всё на hx в квадрате.
33:42
Speaker A
Плюс производную по y расписываю. Это будет с учётом этого шаблона u i j + 1 m - 2 u i j m + u i j - 1 m делить всё на hy в квадрате.
34:12
Speaker A
Аппроксимацию вот этой производной записываю. Это будет что такое? Это будет u i j-тое на m + 1 временном слое u i j-тое на m-том временном слое делить на тау. Ну, то есть, давайте подытожу, понятно.
34:35
Speaker A
Производную по иксу записал на основе вот этих трёх узлов, которые лежат вдоль икса. Конечно-разностную аппроксимацию по игреку записал на основе вот этих трёх узлов, которые лежат вдоль игрека.
34:49
Speaker A
Конечно-разностную аппроксимацию производной по Т написал на основе вот этих двух узлов, которые лежат вдоль оси Т.
34:58
Speaker A
Продолжаем. Здесь у нас что будет? У нас здесь будет минус у и ж м. На всякий случай напоминаю, про что я говорил.
35:13
Speaker A
Вот здесь есть произвол. В каком узле этого шаблона записывать значение у. И этот произвол может повлиять на свойства схемы.
35:26
Speaker A
Что, в зависимости от того, в каком узле берём, схема может обладать либо более высоким порядком точности, либо менее, либо лучшей устойчивостью, либо худшей устойчивостью, либо монотонностью и так далее. Этому как раз и посвящён, напоминаю, курс теория разностных схем, который предлагает вам
35:43
Speaker A
для изучения методологию, которая позволяет все эти свойства исследовать. Сейчас мы просто наивным образом выписываем более простую конечно-разностную аппроксимацию, чтоб нам легче просто было выписывать. Ну, потому что вот здесь, как я сейчас записываю, если всё практически на одном временном слое, у меня будет только одна
36:01
Speaker A
неизвестная. Так. Также, раз уж речь про теорию разностных схем зашла, смотрите. Я, например, вот здесь вам предлагаю записать конечно-разностную аппроксимацию первой производной со вторым порядком точности на основе вот этих двух узлов.
36:22
Speaker A
Но идея в чём? Опять же, в теории разностных схем обсуждается, что вот эта аппроксимация производной с полусуммой, как иногда говорят, она обладает свойством немонотонности. И возможно, иногда имеет смысл записать не со вторым порядком точности, а с первым.
36:40
Speaker A
Ну, например, вот здесь -1 убрать, и внизу будет просто hx. Будет конечно-разностная аппроксимация производной с первым порядком точности, но будет обладать свойствами монотонности. В этом базовом курсе мы вот это не обсуждаем, но на всякий случай я ваше внимание на этот момент
37:02
Speaker A
обращаю. есть, возможно, что это не самая лучшая аппроксимация для вот этой задачи. Но моя цель просто показать общий принцип, как всё расписывается. Ну, и вот здесь тоже аналогично возьму, давайте, вот в этом узле. Но также вопрос дискуссионный, в каком узле брать лучше.
37:30
Speaker A
И смотрите, вот это вот это сейчас было записано, вот это уравнение записывается для каких иксов и игреков, то есть для каких узлов сетки. То есть, мы вот это можем записать, давайте, я вот так сделаю, для всех и от единицы до
37:52
Speaker A
nx -1, потому что граничные точки a и b не включаются, и значит, не включаются граничные узлы 0 и nx. И для всех g от единицы до ny -1.
38:08
Speaker A
Теперь, что мы с вами можем сделать? Естественно, мы с вами отсюда смотрим на вот это уравнение. То есть, это уравнение, как я уже вам говорил, выписано за счёт конечно-разностной аппроксимации вот этого уравнения с помощью конфигурации узлов, определяемых вот этим чисто явным
38:29
Speaker A
шаблоном. И здесь, вот, кстати, давайте здесь записать, да забыл m дописать. И здесь видим что?
38:36
Speaker A
С учётом того, что построено на рисунке, мы видим, что у нас, если мы исходим из предположения, что в текущий момент времени m сеточные значения нам все известны в текущем временном слое, и наша цель найти на следующем временном слое с индексом m +
38:55
Speaker A
1, из этих соображений понятно что? Что это уравнение, несмотря на то, что здесь шесть разных сеточных значений написано, что формально содержит только одну неизвестную, вот это. Всё остальное мы считаем известными.
39:13
Speaker A
Поэтому из вот этого уравнения мы совершенно спокойно, давайте я вот это делать не буду, мы выражаем u и житое с индексом m + 1. Ну, почему я буду не буду делать? Потому что здесь соответствующие преобразования чисто алгебраические, а соответствующая
39:32
Speaker A
формула будет выписана в программке правильно. И что мы с вами делаем? То есть, представьте, вот отсюда, из этого уравнения выразили явно u и житое с индексом m + 1 и в двойном цикле для вот таких значений индексов мы
39:53
Speaker A
можем найти приближённые значения во всех вот этих узлах. То есть, видите, что нарисовано. Вот я приложил шаблон к нулевому временному слою, вот к какому-то конкретному узлу, и с помощью него на следующем временном слое я нахожу, вот красненьким выделено,
40:16
Speaker A
приближённое значение по вот этой формуле. Дальше я могу приложить к другому узлу шаблон и найти соответствующее другое приближённое значение. Теперь очень важный момент. Почему?
40:29
Speaker A
Почему? Я как раз прикладываю вот здесь не последовательно узл шаблоны, а к произвольным точкам, ну удовлетворяющих естественно вот этим индексам. Очень важное замечание какое? Смотрите, вот мы да, в программе, естественно будем формулу для вычисления U и G этого с
40:50
Speaker A
индексом M + 1 записывать в двойном цикле и у нас всё это будет выделяться последовательно.
40:56
Speaker A
Сначала будет вот это значение вычисляться, потом вот это, вот это, вот это, вот это и и так далее.
41:04
Speaker A
Но при вычислении для разных индексов, вот смотрите, какая нравится мне красивая фраза, нету никаких информационных зависимостей. То есть для каждого каждой пары индексов I и G соответствующие вычисления могут производиться независимо друг от друга. А это позволяет сделать что? Например, в
41:27
Speaker A
частности реализовать параллельный алгоритм реализации вот вот этого большого двойного цикла. То есть запускается на многопроцессорной системе программа и каждый процессор обрабатывает свои пары индексов I и G независимо от других процессов. Потом все обмениваются информацией, на у всех есть
41:50
Speaker A
в результате, ну, мы получаем, например, на каком-нибудь управляющем процессоре набор сеточных значений на соответствующем временном слое. То есть, если мы вот этот двойной цикл, ну, в нашей последовательной программе, выполним, мы в результате находим сеточные значения на очередном временном
42:06
Speaker A
слое. Конечно, обращаю ваше внимание, часть сеточных значений вычислена точно из граничных условий, часть сеточных значений найдена приближённо с помощью вот этого шаблона, который даёт так называемую явную схему, которая будет эквивалентна, когда мы будем применять метод прямых.
42:26
Speaker A
Вот эта схема будет эквивалентна использованию одностадийной схемы Рунге-Кутты. Для про неё мы знаем, что первый порядок точности, нету устойчивости, а нету устойчивости, это что значит? Это значит, что на оптимальное соотношение шагов hx, hy и тау, когда схема не разваливается, не
42:48
Speaker A
получается разболка, уводящая решение на бесконечность, нужно накладывать условия. В случае линейных задач, если бы нелинейности не было, это условие Куранта, которое вводится в теории разностных схем. В случае нелинейных задач, явным образом, конечно, условия не выписать, но всё равно условия, накладывающие ограничения
43:10
Speaker A
на взаимные ограничения между шагами hx, hy от и тау, или то же самое между числом интервалов, оно всё равно есть. Явно выписать нельзя, но оно есть.
43:23
Speaker A
А теперь обращаю ваше внимание. Дальше мы что делаем? Дальше процедура та же самая. Нужно найти на очередном временном слое. Поэтому мы прикладываем шаблон, опять же, в общем случае, к произвольным парам узлов, определяющих соответствующей сеткой, и находим сеточные значения на следующем
43:47
Speaker A
временном слое. И так далее. То есть, в результате, то есть, смотрите, если мы в двойном цикле оформляем вот эту формулу, ну, когда выражена U и G Т с индексом m + 1, то да, вот этот двойной цикл по
44:05
Speaker A
приближённым значениям на текущем временном слое с индексом m позволяет найти набор сеточных значений на следующем временном слое с индексом m + 1. Поэтому, если уже вот это всё, вот этот двойной цикл мы оформим ещё третьим внешним циклом с
44:25
Speaker A
индексами пробегающими от нуля до m большое - 1, мы найдём набор сеточных значений на всех временных слоях. Ну, понятно, ещё раз вот обращаю ваше внимание, что есть когда мы делаем для m равного нулю все сеточные значения с индексом ноль,
44:48
Speaker A
на нулевом временном слое, у нас уже известны из начального условия. Потом, когда делаем для m маленького равного m большое - 1, финально мы с вами что находим? Найдём решение в момент времени m большое - 1 + 1 есть
45:04
Speaker A
в момент времени t с индексом m большое, есть в финальный момент времени. А теперь, коллеги, давайте посмотрим на соответствующую программную реализацию. Итак, коллеги, посмотрим на программную реализацию.
45:21
Speaker A
Опять, как всегда, напоминаю, что я демонстрирую интеллектуальную часть кода, преамбулу, подключение пакетов, вывод графики я в основном, как всегда, опускаю. Полная версия программки выложена на портале Тючин, в разделе материалы лекций, и называется пример 27.1. 27 - наша сегодняшняя последняя
45:42
Speaker A
лекция, и один номер примера внутри этой лекции. Итак, обращаю внимание, как я программку писал. Я взял пример 19.1 из девятнадцатой лекции.
45:56
Speaker A
Всё-таки я же обобщал одномерную по пространству задачу на двухмерную, поэтому почему бы уже не использовать шаблон, который. И в начале что сделано?
46:08
Speaker A
Функция, определяющая начальные условия, была зависимость только от x, добавилась y. Функция, определяющая левые и правые граничные условия, была зависимость от t, добавилась y. И добавились функции u top и u bottom, определяющих верхнее и нижнее граничные условия. Под верхним и нижним я
46:29
Speaker A
подразумеваю в том случае, если мы отдельно бы рисовали плоскость xoy. Поэтому на плоскости xoy вот эти как раз функции определяли бы верхнее значение и нижнее в прямоугольнике по иксу и по игреку. Дальше, основная часть программы.
46:49
Speaker A
Задал границы соответствующего прямоугольника по пространству, в котором решается задача, начальный момент времени, конечный момент времени и эпсилон, который, напоминаю, определяет крутизну вот этого переходного слоя, который возникает в задаче. И чем больше порядок малости эпсилон, то есть, если вот здесь будет -2, -3, то есть,
47:12
Speaker A
если эпсилон будет не 1/10, 1/100, 1/1000, 1/10000, крутизна растёт, это будет приводить к необходимости выбирать для решения гораздо более густые сетки.
47:26
Speaker A
Теперь, определил сетки по пространству, ну, число интервалов nx, ny, число интервалов по времени. Определил соответствующую пространственно-временную сетку, hx, шажок по иксу вычислен.
47:40
Speaker A
Задана сетка по иксу, равномерная, как мы делали уже с вами много раз. Аналогично, шажок по сетке по игреку, равномерная сетка по игреку. И то же самое по времени.
47:53
Speaker A
Теперь, какие будут ещё отличия? Теперь массив U, если я напоминаю, у нас структура массива U была какая? Это был двумерный массив.
48:05
Speaker A
Каждая строка содержала набор сеточных значений в соответствующий момент времени. И мы этот массив двумерный заполняли последовательно. Сначала строку с индексом ноль, который соответствует моменту времени Т с индексом ноль, мы заполняли начальным, ну, сеточными значениями, вычисленными с помощью начального
48:26
Speaker A
условия. А потом с помощью схемы заполняли последовательно последующие строчки этого массива. Теперь массив будет уже трёхмерный. Потому что добавилось ещё одно измерение. Ну и маленький нюанс, обращаю ваше внимание.
48:43
Speaker A
Я обычно всегда использовал функцию zeros, которая и выделяет место в памяти под массив и заполняет его сразу нулями.
48:52
Speaker A
Empty делает что? Просто выделила место в памяти под этот массив, ничем не заполняя. Почему я только сейчас решил empty использовать? А потому что здесь массивы достаточно большие и плюс и то есть я не хочу тратить время на их
49:07
Speaker A
заполнение. Ну, хоть хоть какая-то мелкая оптимизация пусть будет. И второе, что если часто в предыдущих программках бывало, что я задавал какой-то двумерный массив, а потом определял только не нулевые элементы, предполагая что нулевые элементы сразу же уже являются нулевыми, когда я массив инициализировал
49:28
Speaker A
с помощью функции zeros, то здесь такого не будет, поэтому я использую вот эту команду empty, которая просто выделяет место в памяти. Теперь, как мы вы здесь писали, в двойном цикле по i и j, которые в пределах от нуля до nx включая
49:44
Speaker A
и от нуля до ny включая, я заполняю массив U для индекса m равного нулю.
49:56
Speaker A
Это первый индекс временной, я его оставил первым. Соответствующими сеточными значениями функции, определяющей начальные условия. И в двойных циклах определяю граничные условия. То есть изменения минимальные, все, думаю, очевидные. И вот основная часть кода.
50:17
Speaker A
Вот я вот здесь выразил u с индексом m + 1 в простра- для пространственных узлов i и j через всё остальное, что мне известно. Вот в двойном цикле по i и j осуществляется пересчёт пересчёт соответствующих сеточных значений в момент времени t с
50:41
Speaker A
индексом m + 1 через сеточные значения в момент времени t с индексом m. И потом во внешнем цикле что реализуется?
50:50
Speaker A
Реализуется пе- нахождение сеточных значений на каждом очередном временном слое. Поэтому, когда я говорил про информационные зависимости, вот на что я обращаю ваше внимание. Смотрите.
51:03
Speaker A
Для разных индексов i и j соответствующие вычисления делаются независимо друг от друга. То есть не нужно вычислять предварительно какое-нибудь значение с индексом друг- ну, с другими индексами i и j, которые используются вот в этом текущем вычислении. А вот для вычислений для
51:23
Speaker A
разных m, понятно, информационные зависимости у нас есть. Мы не можем вычислить, например, u с первым индексом 10, если у нас нету сеточных значений для u с индексом первым, например, девять. Поэтому для различных m вот эти вычисления независимо мы проводить не
51:43
Speaker A
можем. То есть по m у нас информационные зависимости есть. По i и ж, информационных зависимостей нету.
51:50
Speaker A
Теперь, я как раз обращаю ваше внимание на это, почему? Потому что мы видим, что именно явная схема реализовалась очень просто, очень просто. То есть обобщение, вот, я не побоюсь этого слова, элементарные, элементарнейшие. Это с чем связано? Это будет связано с тем, что вы увидите,
52:07
Speaker A
когда мы будем применять метод прямых, чтобы построить неявную схему, которая обладает большим порядком точности, большей устойчивостью. Всё будет гораздо более громоздко и кропотливо. Здесь всё действительно очень просто. Обобщение элементарное, но минус, что явная схема с первым порядком точности и с
52:29
Speaker A
отсутствием безусловной устойчивости. Устойчивость будет в каком случае? Если nx, ny вы берёте очень большими. Поэтому обращаю ваше внимание, что на практике часто используется следующий подход.
52:44
Speaker A
Но при э-э реализации некоторых реальных задач выписывается явная схема и потом, за счёт отсутствия вот этих информационных зависимостей, про которые я вам сказал, реализуется её распараллеливание. Потому что явная схема распараллеливается относительно просто.
53:04
Speaker A
Потому что нужно организовать, если у вас есть много вычислительных устройств, нужно распределить между ними пары индексов и и ж, для которых будут проводиться соответствующие вычисления. Кстати, вот именно вот это распараллеливание, вот этой программки, мы проделываем в курсе параллельные вычисления. Там в качестве
53:27
Speaker A
примера, на относительно несложного, как раз вот эта программка и рассматривается, на которой демонстрируются некоторые особенности процедуры распараллеливания, ну и алгоритмов распараллеливания. Вот. А вот неявные схемы, они сами по себе громоздкие и гораздо более времязатратные для распараллеливания.
53:48
Speaker A
Поэтому вы поняли идею, почему я вот сейчас обращаю ваше внимание на вот эту вроде бы примитивную и не сверхнадёжную схему, потому что когда вы решаете именно реальные задачи, одно дело, понимаете, на что хочу тоже обратить внимание ваше, на очень важный момент.
54:04
Speaker A
Одно дело есть теоретический результат, который говорит: "Эта схема не особо, эта схема гораздо лучше". Это одно дело, а другое дело это практическая реализация. На практике бывает какую-нибудь замечательнейшую схему реализовать очень сложно либо невыгодно.
54:22
Speaker A
Поэтому здесь всегда вопрос стоит об эффективности. Какой результат вы получите, сколько времени вы затратите. И понима- и нужно уметь находить баланс. Здесь баланс какой?
54:34
Speaker A
Явная схема не особо надёжная, нужно брать очень большие числа nx, ny и m, чтобы всё хорошо считалось, но при этом расписывать схему легко. Понятно, дальше за вас схему расписываете вы ручками на бумажке, а считает за вас компьютер.
54:57
Speaker A
А во втором случае, когда мы будем говорить про неявные схемы, да, компьютер будет считать быстро, но вы потратите огромное число времени на расписывание схемки. То есть для сравнения просто скажу, вот эту схемку расписать у меня ушло, понятно, аккуратно вот эти все выкладки дома
55:13
Speaker A
минут 10-15, ту схему, которую мы будем чуть позже смотреть, у меня ушло часа 4-5 аккуратно все выкладки сделать. Понятно, если речь будет идти про трёхмерные задачи или про что-то более сложное, про оптимизацию вычислений, уже количество ручных действий возрастает. А помним общий
55:32
Speaker A
принцип какой тоже? Чем больше действий вы совершаете сами вручную, тем больше вероятность совершить ошибку, Ну потому что она примерно прямо пропорциональна количеству выполняемых вами действий. А вот для больших n y n x m всё-таки будет делать компьютер. Компьютер, если схема
55:51
Speaker A
протестирована, он будет считать, понятно правильно. Вот. А сейчас скоро, коллеги, мы посмотрим на то на видео, как выглядит вообще решение, что выдаёт эта программа. Итак, коллеги, то есть если я обычно, когда демонстрировал код, куски кода, соответствующие графическому выводу, опускал, они
56:14
Speaker A
присутствовали только в полных версиях программы, выложенных в разделе материалы лекции портала ТЕЧИН, то здесь всё-таки на соответствующие кусочки кода я хочу обратить внимание, потому что вот давайте, у нас два варианта. Один вариант вот такой, всё-таки я вам его
56:31
Speaker A
демонстрирую. И давайте, кстати, посмотрим на видео. Итак, видео. Он реализует построение вот такого вот графика трёхмерного. То есть, понятно, смысл этого видео как какой? Что каждый кадр, каждый кадр соответствует поверхности с координатами x y, определяющей значение функции в соответствующий
56:55
Speaker A
момент времени. И вот для всех m мы строим соответствующую картинку. Что-то здесь происходит, какая-то там волна бежит, но, например, лично мне вот такие графики не особо нравятся, мне трудно что-либо рассмотреть. Поэтому сейчас будет ещё чуть позже одно видео,
57:13
Speaker A
которое, давайте с этим видео закончим, с этим закончим, и вот другой вариант построения видео.
57:20
Speaker A
Итак, видео посмотрим. Здесь мне нравится это видео чем? То есть здесь плоскость x y, а каждой сеточке, ну Ну узлу сеточки соответствует цвет. Цвет соответствует конкретному значению решения. Здесь, несмотря на то, что решение, как мы видели на трёхмерном графике, меняется
57:43
Speaker A
мало, достаточно, вот здесь более явно видно, что происходит. Вопрос, вопрос, конечно, как строить соответствующий график, как интерпретировать результаты счёта, в том виде график строить или в этом, это абсолютно дело вкуса. Также, кстати, коллеги, ещё про дело вкуса маленький комментарий. Обычно, если вы
58:06
Speaker A
оформляете какой-то результат, который представляет какую-то значимость, естественно, вот эта картинка бессмысленна без чего? Без колорбара.
58:13
Speaker A
Под колорбаром что подразумевается? Функция, которая выводит цвета, которые соответствуют, написано, каким конкретным значениям. Почему я колорбар не вывожу? Потому что здесь у нас просто это какой-то модельный пример, конкретные значения не имеют для нас никакого ну, сейчас с вами значения.
58:32
Speaker A
А так, если вы, естественно, получаете какой-то научный результат, нужно ещё добавлять и колорбар. Ну, колорбар я не добавил, чтобы не перегружать картинку.
58:41
Speaker A
Коллеги, переходим теперь ко второму подходу к решению. Также, как в девятнадцатой лекции, такая же структура. Рассмотрели простейший, интуитивно понятный подход, который к тому же и очень просто реализуется. А теперь возникает вопрос: а всё-таки давайте построим надёжные схемы решения, которые являются
59:02
Speaker A
безусловно устойчивыми. У нас уже отработанный подход в этом курсе, который в чём заключается? Есть система для уравнения в частных производных, и что мы дальше делаем? Мы вводим в сетку сейчас только полупространственный переменный. То есть, обращаю ваше внимание, что по Т
59:25
Speaker A
уже сетка отсутствует. И после этого что мы делаем? Мы запишем конечно-разностные аппроксимации только пространственных производных.
59:40
Speaker A
После того, как мы это сделаем, у нас получится свести вот эту задачу для уравнения в частных производных к очень большой системе обыкновенных дифференциальных уравнений относительно функции U с фиксированными координатами X Y, которые соответствуют узлам этой сетки.
60:06
Speaker A
И в результате это будет приводить к тому, что нам нужно восстановить все эти неизвестные функции вдоль вот этих вот прямых. Поэтому давайте, то есть напоминаю вам ещё раз формально подход. Опять же, наша цель найти функцию U от X
60:28
Speaker A
Y T. Желательно аналитически. Этого, к сожалению, мы сделать не можем. Поэтому опять, как мы и делали, ввели сетку только по пространству, то есть по иксу и по игреку.
60:44
Speaker A
Ну, давайте формально напишу, вдруг кто-то отмотал лекцию сразу к вот этому кусочку. То есть X и т определяем как A + H X на и, и от нуля до N X, H X это, соответственно, B - A / N X.
61:06
Speaker A
Сеточку по игреку, Y житое это C + H Y на G, G от нуля до N Y, H Y это D - C / на число интервалов вот этой сетки. И теперь наша цель найти функцию не x, y, t,
61:27
Speaker A
а функции, вот обозначения помним какие, функции с индексами i и j только аргумента t для всех возможных пар узлов i, j, которые образуют вот это пространственная сетка.
61:48
Speaker A
И как мы постоянно с вами делаем, используем какую договорённость? Что если у нас аргумент t, то мы его опускаем.
61:59
Speaker A
То есть если аргумент будет отличен от t, мы запишем конкре- в скобочках конкретное значение аргумента.
62:07
Speaker A
Вот это это подразумевается, что вот эта функция зависит от t. Соответственно, что мы с вами дальше делаем? Дальше мы записываем конечно-разностную аппроксимацию вот этого уравнения в соответствующем пространственном узле сетки с индексами i и j. Что у нас получается? Ну, по аналогии, что мы
62:34
Speaker A
сейчас писали. Эпсилон u i плюс один j-тое, то есть предполагается, что зависит от t, минус два u i j-тое плюс u i минус один j-тое делить на hx в квадрате плюс, вот здесь я плюс случа- минус случайно записал, здесь должен быть
62:59
Speaker A
плюс, плюс u i j плюс первое минус два u i j-тое плюс u i j минус первое делить на hy в квадрате минус D U и житое по Dt то, что кучу раз обсуждали.
63:28
Speaker A
Здесь U функция нескольких переменных. Теперь вот здесь, когда записываем в конкретном и житом узле, это функция одного, вот эта функция одного переменного t, поэтому производные стоят полные.
63:43
Speaker A
Равно минус U и житое здесь U и плюс один, записываем конечно-разностную аппроксимацию первой производной по иксу, U и плюс один житое минус U и минус один житое делить на два h x плюс аппроксимация по игреку, U и жи плюс
64:08
Speaker A
один минус U и жи минус один делить на два h y и минус U и житое в кубе.
64:23
Speaker A
Естественно, вот это уравнение оно не одно, оно записывается для каких и и жи? Оно записывается для и от единицы до nx минус один, то есть узлы ноль nx не включаются, потому что здесь выколоты a и b.
64:42
Speaker A
Жи аналогично от единицы до ny минус один t остаётся от t0 до t большого.
64:55
Speaker A
Теперь учитываем граничные условия записываем аппроксимации, это точка a соответствует по иксу узлу с индексом ноль то есть U ноль житое это U left вычисленная в точке Y житое от T.
65:17
Speaker A
Аналогично здесь, U с индексом N X с индексом I, это U right от Y житое T.
65:38
Speaker A
G вот здесь от единицы до N Y - 1, потому что выколота. T от T 0 до T большого.
65:52
Speaker A
И вот эти граничные условия. U I 0 C соответствует по Y узлу с индексом 0.
66:06
Speaker A
U bottom от Y Опечатался. Здесь уже, естественно, не Y, а здесь у нас X X и T T.
66:25
Speaker A
И U с индексом I здесь N Y, это U top с индексом X I T T.
66:42
Speaker A
И вот здесь I, смотрите, уже от нуля до N X граничные узлы включаются, потому что здесь включаются граничные точки этого отрезка.
66:53
Speaker A
T принадлежит от T 0 до T большого. И дописываем начальные условия. U с индексом I G вот здесь аргумент конкретно Т0, потому что он отличается от Т, поэтому я его явно пишу. Здесь будет У initial от X итого
67:19
Speaker A
Y житого и и изменяется в пределах от нуля до NX, G изменяется в пределах от нуля до NY.
67:33
Speaker A
Всё. Методом прямых свели задачу к большой системе обыкновенных дифференциальных уравнений. Количество дифференциальных уравнений у нас здесь сколько?
67:50
Speaker A
Единиц NX - 1 умножить на NY - 1. Неизвестных функций у нас во всех узлах это NX + 1 умножить на NY + 1.
68:03
Speaker A
есть вот именно только вот эти уравнения дают переопределённую систему. Но с учётом граничных условий у нас количество неизвестных функций и количество уравнений совпадает.
68:16
Speaker A
Напоминаю формальную терминологию, какую? Формально мы получили не чисто дифференциальную систему, а дифференциально-алгебраическую, потому что вот это дифференциальные уравнения, потому что вот у нас присутствуют производные вот этих функций по Т, но есть ещё и алгебраические выражения.
68:38
Speaker A
Что мы с вами делаем? А мы с вами делаем обычно следующее. Мы подставляем вот эти известные значения функций, заданные граничными условиями, в соответствующие уравнения и в результате сводим систему к чисто дифференциальной. Давайте частично с вами это проделаем. Почему
69:05
Speaker A
частично? Потому что, смотрите, здесь какой момент интересный. Вроде бы пока мы ничего такого особенного с вами особо сложного не сделали. Но вот дальше вдруг всё начнётся очень громоздко, очень громоздко.
69:18
Speaker A
И смотрите, что я дальше сделаю. Смотрите, что я сделаю. Давайте-ка я запишу выражу d u с индексом и же по d t, оставлю в одной стороне, всё остальное в другую сторону.
69:35
Speaker A
В результате что я получу. Я получу следующую конструкцию. d u и же по d t.
69:45
Speaker A
Так, давайте прицелюсь. Я представляю, что перенёс с плюсом, значит, в другой стороне вот это с плюсом, вот это тоже с плюсом идёт, это с плюсом. В результате будет вот такая вот конструкция.
69:57
Speaker A
Давайте ещё также для удобства, чтобы не писать эпсилон hx в квадрате, я всё за скобочки выношу, ну, чтобы дроби более компактные у меня были.
70:06
Speaker A
Значит, здесь будет u и плюс один житое минус два u и житое плюс u и минус один житое плюс вот это слагаемое.
70:20
Speaker A
Ну, если скобочки раскрыл, эпсилон делённое на h y в квадрате. Ну, вопрос раскрывать скобки или как переписывать, это понятно, абсолютно вопрос вкуса, как вам нравится. Здесь будет u и же плюс один минус два u и житое плюс u и
70:43
Speaker A
же минус один. Плюс вот эта вот конструкция. Давайте, ну, единица напишу, единица на два hx вынесу.
70:55
Speaker A
Здесь u и житое. У и плюс один житое минус у и минус один житое плюс единица делённая на два h y у и житое у и жи плюс один минус у Я вот стараюсь записать, мне нужно ещё, представьте, будет немало места
71:33
Speaker A
справа. У и жи минус один и минус у и житое, точнее, плюс. Я же перенёс всё с плюсом плюс с плюсом. Теперь смотрите, я вот это записываю вот это.
71:54
Speaker A
Делаю последовательно. Я это записываю не для и от единицы до nx минус один, а вот смотрите, внимание, от двух до nx минус двух. А почему?
72:11
Speaker A
А потому, что именно при и равном единице и и равном nx минус один, я перепишу вот это выражение, подставив в него, соответственно, левое и правое граничные условия. Потому что, представьте, если у меня индекс и будет равен единице, то вот здесь возникает у с первым
72:38
Speaker A
индексом ноль, то есть ноль жи, заменён на вот это. И аналогично при и равном здесь возникнет у с индексом nx, запятая, жи. Учту вот это граничное условие. Поэтому давайте я вот это сделаю и запишу конкретно.
72:59
Speaker A
d u, конкретно я для i равного единице пишу, поэтому смотрите внимание, здесь везде u i точечка, а здесь я вот палочкой ставлю, чтоб единичка. u 1 j по dt равняется эпсилон делённое на hx в квадрате u 2 j -2 u 1 j
73:23
Speaker A
плюс вот здесь возникает u с индексом 0 j меняю на u left от y житого t плюс эпсилон делённое на hy в квадрате. Вот это перепишется без изменений, потому что с соответствующими индексами у нас ничего вот именно здесь не появится,
73:50
Speaker A
поэтому будет u 1 j + 1 - 2 u 1 j плюс u 1 j - 1 плюс 1 делённое на 2 hx u 1 j здесь u 2 j, а здесь u 0 j заменяю на u left
74:27
Speaker A
от y житого t плюс h и вот это слагаемое без изменения останется 2hy u 1 j u 1 j + 1 - U1 G - 1 и + U1 от 1, а то я здесь на случайно сначала И
74:58
Speaker A
написал, в кубе. То есть поняли, вот это мы выписали, вот здесь давайте я опущу.
75:06
Speaker A
То есть вот это выписано для вот этих И, вот это выписано для И, для И равного единице. И отдельно пишу для И равного Nx - 1, для И равного Nx - 1. DU с индексом Nx - 1 G
75:32
Speaker A
по DT равно эпсилон делённое на Hx в квадрате. Здесь смотрите, когда И равно Nx - 1, будет стоять здесь U с индексом NxG, меняю на правое граничное условие U right от Y житого T -2 U Nx - 1 G
76:08
Speaker A
+ U с индексом Nx -2 G + эпсилон на формально вот эти страшные выражения переписывать.
76:23
Speaker A
Здесь Nx - 1 G + 1 - 2 U Nx 1 G + U Nx - 1 G - 1. Конечно, возможно, я сейчас где-нибудь могу опечатку допустить, но я надеюсь, что вы её тогда заметите, а сама в самой программе всё, естественно,
76:48
Speaker A
написано правильно, потому что когда я дома всё делал на листочке, я перепроверял кучу раз и причём находил у себя, само собой, ошибки.
76:57
Speaker A
Потому что делая много таких операций, ошибку допустить, ну, точнее, опечатку, не мудрено. Вот здесь, когда пишу U с индексом nx - 1 и ещё + 1, это будет U, опять же, вот это с индексом nxj. Подставляю опять вот это правое граничное условие.
77:21
Speaker A
Y житое T U n x - 2. Кстати, коллеги, сейчас проверю, nx, да, да, да, да, здесь не забыл. Я просто подумал, что вдруг я nx - 2 в предыдущем выражении забыл.
77:41
Speaker A
g + 1 / 2 h y U с индексом nx - 1 g U с индексом n x - 1 g + 1 - U с индексом nx - 1 g - 1 и + U с индексом nx - 1
78:10
Speaker A
g в кубе. Всё, написал. А теперь смотрите интересный момент. Громоздко, да, громоздко. Но почему я вот до сих пор не дописал, это не просто так. Почему я не написал вообще-то вот вот вот эти все уравнения, они причём не только при i вот в таких пределах. То
78:30
Speaker A
есть, сейчас смотрите, вот вот вот здесь что переписано? Здесь переписана вот эта конструкция для всех i просто в другой форме записи. Просто для и равного единице мы переписали в таком виде, учтя U left, левое граничное условие. Для и
78:44
Speaker A
равного nx - 1 тоже пере- переписали, учтя правое граничное условие. А теперь вот это всё вот это всё у нас записывается для вот да, у меня сейчас не уместится для чего?
79:06
Speaker A
Для G от единицы до ny 1. Кстати, обращаю ваше внимание, что и вот здесь, когда мы подставляли граничные условия, у нас тоже G от единицы до ny - 1, поэтому мы ничего не потеряли, ничего лишнего не записали. Но
79:28
Speaker A
мы сделаем с вами что? Мы это запишем для G от двух до ny - двух. А с чем это связано? А это связано с тем что для G равного единице уже вот этот весь блок вот здесь переписывается учитывая вот эти два гра- вот это
79:53
Speaker A
граничное условие. Ну, то есть, поняли, когда G будет равно единице вот этот весь блок переписывается и вот здесь, например, будет содержаться индекс ноль. Ноль. Ноль.
80:06
Speaker A
Можем заменить соответствующие выражения будем на вот это нижнее граничное условие. И вот этот весь блок ещё записываем для G равного ny - 1. В этом случае вот здесь тоже будет содержаться индекс второй ny и используем верхнее граничное условие.
80:28
Speaker A
Поэтому уже понимаете, громоздко вот это я с вами сейчас выписывать не буду. Я думаю, принцип понятен, принцип понятен, и вы сможете при желании сами соответствующие формулы выписать. Они, естественно, в программный код у нас перенесены. Поэтому давайте ещё раз
80:43
Speaker A
подытожим. Что мы с вами делаем? А мы с вами делаем следующее. После вот этого. Вот этот блок из трёх уравнений записываем сюда для G равного единице. Везде, где второй индекс ноль встречаем, заменяем на U bottom. И вот этот весь блок из трёх
81:03
Speaker A
уравнений записываем для G равного Ny - 2. Везде, встречая второй индекс Ny, используем верхнее граничное условие.
81:15
Speaker A
Уже получается очень громоздко, уже очень громоздко. То у нас что у нас получается девять различных вариантов записей.
81:27
Speaker A
Поэтому я вот здесь, поняли, специально систему не пишу, потому что сейчас у меня её никакого эквивалентного нету.
81:31
Speaker A
Эквивалентный переход в каком случае был бы? Блок из трёх уравнений здесь, блок из трёх уравнений здесь, но под тремя уравнениями вы понимаете, я подразумеваю не само уравнение, а форму его записи, потому что уравнений у нас само собой гораздо
81:45
Speaker A
больше. И ещё плюс начальные условия. В результате получается очень большая система обыкновенных дифференциальных уравнений.
81:59
Speaker A
То есть таким образом мы всегда с вами что сделаем? Мы с вами вводили, вот сейчас я допущу некоторую неточность и поправлю её, вводили вектор Y и показывали, что вот эту огромную систему обыкновенных дифференциальных уравнений мы можем записать вот в таком вот виде.
82:28
Speaker A
dy по dt равно f от y t t принадлежит от t0 до t большого y от t0 равно y0.
82:49
Speaker A
И дальше предлагали решать вот эту систему обыкновенных дифференциальных уравнений, полученную, с помощью одной схем из семейства схем Розенброка. Рекомендуют рекомендовали либо обратную схему Эйлера либо схему Розенброка с комплексным коэффициентом как наилучший вариант. Ну, только смотрите, какие я здесь
83:12
Speaker A
неточности допустил. У меня всё-таки уже y зарезервирован за второй пространственной переменной, поэтому давайте именно вот здесь, в этой программе, ну, в этой реализации, я вот здесь y давайте на какую-нибудь другую переменную заменю, а именно p, например.
83:28
Speaker A
Ну, понятно, здесь почему именно p? А потому что первое, что в голову пришло, ну, из того, что не зарезервировано в программках.
83:39
Speaker A
Вот. И смотрите, какая у нас структура вектора p. А структура вектора p, ну, смотрите, у нас же слева стоит dp по dt, а в какой последовательности здесь у нас записаны u в соответствующих узлах сетки. Ну, как мы с вами уже обсудили, сначала у нас
84:01
Speaker A
вот здесь будет стоять блок для g равного единице, поэтому в этом блоке для g равного единице здесь будут стоять единички.
84:10
Speaker A
А потом будет стоять, ну, представьте, если этот блок, вот этот блок для g равного единице один, ну, здесь везде один стоит, а здесь один потом два, три, четыре, пять и так далее.
84:22
Speaker A
Поэтому структура будет следующая. У 1 1 У 2 1 и так далее У с индексом nx -1 1.
84:39
Speaker A
Потом следующий блок для G равного двум. Для G равного двум всё это повторяется. Будет U 1 2 U 2 2 и так далее U с индексом nx -1 2.
84:54
Speaker A
Потом для третьего блока для четвёртого и так далее. Всё закончится n x -1, а здесь n y -1.
85:09
Speaker A
То есть будет фактически набор вот этих сеточных значений, которые нам интуитивно легко представлять в виде двухмерной матрицы, он у нас разворачивается в виде вот такого вот вектора. Поняли?
85:25
Speaker A
Сначала зафиксирована второй индекс G и изменяется первый индекс от единицы до nx -1. Потом второй индекс стал равен двум, изменяется индекс от единицы до nx -1 и так далее.
85:42
Speaker A
Да. Если бы мы записывали, понятно, структуру в друг ну могли бы последовательность менять, например, не как я делал вот здесь по i, а там по G, у нас, понятно, структура бы поменялась. Если меняются местами уравнения, это бы привело к тому, что и
85:57
Speaker A
структура вектора P поменялась. Коллеги, смотрите. Вообще общая концепция всего вот этого курса какая была? Что мы не просто с вами рассматривали численные схемы, а мы их в том числе рассматривали на конкретных примерах, чтобы вы могли всё без проблем
86:14
Speaker A
применить. Теперь же вот именно вот в этой лекции всё получается настолько громоздко, что я многие выкладки вынужден пропускать, а объяснять в общих чертах, что происходит.
86:26
Speaker A
Но я исхожу из того, что большинство выкладок, которые я делаю, достаточно понятные. А если что-то будет непонятно, те выкладки, которые я пропущу, в комментариях к лекциям я всё время без проблем отвечу. То есть понятно, я пропустил блок, как задаётся
86:45
Speaker A
вот дополнительные три уравнения. Дальше, что я пропускаю? Предполагаю, что вы это сможете сделать, так как мы вот этими действиями занимались, начиная с восемнадцатой лекции. Если есть структура вектора P а у вектора P, ну помним, в памяти будет храниться как? Вот это будет храниться в
87:07
Speaker A
ячейке с индексом ноль массива P в Питоне. Это в ячейке с индексом один и так далее. Ну, например, вот это будет храниться в ячейке с индексом nx - 2.
87:20
Speaker A
Вот это будет храниться в ячейке с индексом nx - 1 и так далее. Всё закончится на вот это будет храниться в ячейке с индексом nx - 1 умноженное на ny 1 и ещё, понятно 1. То есть с учётом нумерации с нуля и
87:44
Speaker A
заканчивая вот этим, у нас всего элементов nx - 1 на ny - 1. Я, естественно, предполагаю что? Что вы сможете выписать. Да, это кропотливо, сразу скажу, у меня дома на вот это ушло минут 40, минут 40, чтобы аккуратно
88:03
Speaker A
выписать компоненты вектора F через P, как мы делали с вами много раз, начиная с восемнадцатой лекции, то есть с по двадцать пятую. Представьте, 40 минут у меня на это ушло. Да, я делал кропотливо. Сейчас бы я бы с этим
88:17
Speaker A
справился бы полчаса, но это всё равно чисто техническая работа. То есть, выписали компоненты вектора F. В программке, естественно, все формулы в полной версии будут выписаны. В версии, которую я сейчас скоро буду демонстрировать, будет выписано не всё, потому что там настолько много формул,
88:35
Speaker A
что они просто не помещаются в виде нормальной презентации. И смотрите, какие мы применим ещё приёмы. Особенность какая? Особенность заключается в следующем. Как мы решаем вот эту систему? Мы её решаем с помощью одной из схем из одностадийного семейства схем Розенброка.
88:57
Speaker A
Сначала находим W1 из вот такой вот системы. Здесь у нас вектор W 1. Справа у нас F от P в m-тый момент времени Tm плюс тау пополам. Нашли W1 и после этого находим вектор P в следующий момент времени. А помним, у
89:39
Speaker A
нас P зависит от T. Почему? Потому что все компоненты зависят от T. В момент времени, ну, в следующий момент времени, находим как в текущий момент времени плюс тау на действительная часть W1. И делаем это для всех m от нуля до m большое минус 1.
90:04
Speaker A
Нам нужно, чтобы это реализовать, нужно будет матрица Якоби. Нужна будет матрица Якоби. А представьте, какая структура у матрицы Якоби будет? Матрица Якоби, за счёт вот этой вот специфики, будет состоять из n y одного блока, размерность каждого из которых
90:37
Speaker A
n x - 1. И то же самое, понятно, матрица квадратная, будет то же самое и с количеством строк.
90:48
Speaker A
n y - 1 блок, каждого размерность n x - 1. Если вы выпишите элементы матрицы Якоби, то окажется, что она сильно разреженная. Ну, разреженная матрица. Большинство элементов нулевые.
91:07
Speaker A
Не нулевые элементы будут стоять на главной на диагонали, на двух ко диагоналях, причём внутри вот этих матричек, то есть, видите, обращаю ваше внимание, я не прочертил вот так от начала до конца, то есть, внутри вот этих блоков только будут три диагональных элемента.
91:27
Speaker A
И ещё, за счёт того, смотрите, ну, например, вот здесь у нас относительно вектора u и житого, почему вот эти три диагональки будут? Потому что есть вот соседние элементы с индексами i + 1 и - 1, которые вот здесь располагаются рядом, и поэтому в
91:44
Speaker A
матрице Якоби они будут располагаться рядом. Но вот здесь относительно i и житого есть элемент, у которого индекс j + 1, а он с индексом j + 1, ну, например вот u 1 1, а с индексом u 1 2 лежит дальше
91:59
Speaker A
на n x - 1. Поэтому за счёт вот этих элементов, где второй индекс сдвинут на единичку, добавится вот такая вот диагональ. И где на единичку вниз второй индекс сдвинут, вот такая вот диагональ.
92:16
Speaker A
То есть матрица системы имеет вот такой вот специфический вид. Она пятидиагональная. И как решать вот эту систему? То есть на что я обращу ваше внимание как раз в программной реализации.
92:33
Speaker A
Я буду обращать ваше внимание на следующее. Что метод Гаусса применять очень плохо, почему? Потому что если у нас размеры матрицы количество строк и количество столбцов равно чему? nx - 1, ny - 1. Ну, например, по порядку. По порядку nx, ny
92:51
Speaker A
у нас 100. Значит, количество строк, ну столбцов 10 в четвёртой. Количество строк 10 в четвёртой. Элементов матрицы 10 в восьмой. Матрица очень большая, в память может не влезть, если, например, вы уже возьмёте коли- количество интервалов 10 в третьей. Тогда и количество строк, и
93:12
Speaker A
количество столбцов будет 10 в шестой, а ра- количество элементов матрицы уже будет 10 в двенадцатой. Уже просто может ну не не поместиться в оперативную память вашего компьютера. Поэтому мы будем в программной реализации применять приёмы, связанные с оперированием матриц, которые мы будем хранить в
93:36
Speaker A
разреженном формате данных. То есть матрицу мы будем хранить не в виде двухмерного массива, а мы разреженный формат данных тоже вкратце, так это как технический вопрос к численным методам ну меньшее имеет отношение.
93:52
Speaker A
Мы будем хранить как? То есть три вектора. Один вектор содержит координаты, тоесть номера строк, в которых находятся эти ненулевые элементы. Второй вектор содержит координаты столбцов, в которых находятся эти ненулевые элементы. Третий вектор нахо- содержит значения этих ненулевых элементов. И поэтому уже видим
94:20
Speaker A
в разреженном формате хранения данных нам уже будет нужно не nx - 1 на ny - 1 в квадрате элементов хранить, а давайте я, ну, округлю, а примерно вот здесь пять на вот это число, потому что это количество строк. И в каждой строке,
94:41
Speaker A
почему я написал примерно пять? Ну, потому что в первой строке три элемента, во второй четыре, дальше пять идёт, вот здесь опять там встретится четыре элемента, ну, и так далее. Тоесть, количество элементов гораздо меньше. И дальше операции с разреженными матрицами
94:56
Speaker A
позволяют ещё более быстрая более быструю реализацию, например, произведения матрицы на вектор. Когда вы перемножаете, помните, матрицу на вектор, вы же фактически реализуете скалярное произведение строки на соответствующий столбец. А вам, если вы знаете, в каких местах содержатся соответствующие ненулевые элементы, у
95:18
Speaker A
вас вычисление соответствующего элемента, ну, результат произведения матрицы на вектор будет требовать скалярного произведения, которое требует, ну, если я вот это обозначу как n, не порядка n операций, а будет только порядка пять, потому что вам нужно будет пять элементов перемножить, ну, их сложить.
95:40
Speaker A
Ну, типа, пять и умножение, четыре сложения. Ну, де- девять операций. И так для всех или строк n будет девять девять n операций, то есть перемножение матрицы на вектор, если матрица хранится в разреженном формате данных, выполняется гораздо быстрее, не порядка
96:04
Speaker A
n в квадрате, а за порядка n операций. Итак, коллеги, давайте вкратце посмотрим на программную реализацию.
96:13
Speaker A
Программа, полная версия, приложена в разделе материалы лекций, пример 27.2. Смотрите, вас сразу должны немножко смутить эти многоточия. А многоточия с чем связаны? А связаны с тем, что, смотрите, вот это, вот это написано, это как раз ровно вон тот блок для g
96:37
Speaker A
от было написано от двух до n минус n y минус двух. Здесь, с учётом структуры вектора b, я был вынужден, чтобы было всё удобней, сдвинуть индексы, чтобы здесь индексация не поменялась. И в результате вот такой вот страшный вид
96:56
Speaker A
соответствующего блока. В полной версии программы, естественно, присутствует не только преамбула, когда подключаем пакеты и вывод графики, но также присутствует и полная версия вот этой функции f, которая в три раза больше, то есть ещё вот три вот таких огромных
97:13
Speaker A
блока добавлены вот здесь и вот здесь. Теперь вы понимаете, почему я не показываю всё целиком. И так здесь уже даже разобраться сложно, вот так глядя, ну естественно.
97:24
Speaker A
Но аккуратно всё сделано, всё соответствует тем формулам, которые выписаны. Единственный нюанс, и да, нужно было сместить индексы.
97:34
Speaker A
Так как мы это делали с вами не раз, я сейчас подробно это не разбираю. Если у кого-то будут вопросы в комментариях к лекциям, спокойно отвечу на всё.
97:45
Speaker A
Дальше. Основная часть программы. От примера 27.1 отличается очень мало. Ну, естественно, также, так как мы используем метод прямых, например, здесь вводится параметр альфа, который определяет схему Розенброка с комплексным коэффициентом. Можете положить единицу, это будет обратная схема Эйлера, 1/2 - схема
98:06
Speaker A
Кранка-Николсона или ноль - одностадийная схема Рунге-Кутты либо схема Эйлера. Соответственно, массив U, который будет содержать сеточные значения, и вот этот вспомогательный массив P других размеров.
98:21
Speaker A
А теперь смотрите, какой очень важный момент. Глядя на вот эту структуру вектора F, на которую я потратил на аккуратное расписывание там 40-50 минут, вы понимаете, выкладок сделано очень много.
98:35
Speaker A
Чем больше выкладок мы делаем, есть вероятность ошибиться просто при банальном Ctrl Copy, Ctrl Paste, который я использовал. Поэтому я предлагаю следующий подход. Я предлагаю следующий подход: не спешить с вычислением якобиана, а сделать что сначала, например, вам? Вычислить вот здесь реализовать просто схему
98:51
Speaker A
Эйлера, которая не требует вычисления якобиана. Если вы на тестовом примере, ну, на достаточно густых сетках, для маленького отрезочка времени убеждаетесь, что у вас программа работает, значит, вектор F вы сформировали правильно. И значит, можно приступать к вычислению якобиана.
99:11
Speaker A
Поэтому вот более полная версия программы, когда я убедился в том, что вектор F сформирован правильно, вот здесь используется из пакета дополнительный пакет из а-а, пакета SciPy функция spsolve. То есть, она решает систему с матрицей, заданной в разреженном формате данных.
99:37
Speaker A
То есть, гораздо более быстро, чем метод Гаусса бы работал бы. Я использовал вот эту функцию, поэтому здесь используется функция F и ещё формируемая матрица А, единичная матрица минус альфа тау на якобиан. И вот эта конструкция выглядит следующим образом.
99:55
Speaker A
Я вот здесь формирую три массива, один массив вот такой вот размерности по количеству ненулевых элементов, то есть все эти три массива содержат, они размера количество ненулевых элементов.
100:08
Speaker A
Это можно посчитать. А дальше я эти массивы заполняю следующим образом. Например, К равно нулю.
100:16
Speaker A
Номер строки, в котором элемент находится, номер столбца и соответствующее значение ненулевого элемента для якобиана. Здесь я привёл только три строки. Таких строк, таких конструкций в полной версии программы 33 штуки. На это у меня ушло где-то 2,5 часа.
100:37
Speaker A
Поэтому, естественно, эти формулы я здесь расписывать не буду, иначе получится жуткий перегруз. Но моя цель просто сформировать, чтобы вы имели представление про сложность уже решения двухмерных задач. Дальше, с помощью команды из пакета спарс, ну, пакет, точнее, сайпай, а в
100:56
Speaker A
нём подпакет, модуль спарс, я сформировал матрицу в так называемом координатном формате, а потом её перевёл в так называемый формат ЦСР, компрест спарс роу, то есть удобный формат, в котором наиболее быстро реализуются вычисления произведения матрицы на вектор.
101:24
Speaker A
И дальше сформировал матрицу А. Это единичная матрица тоже, команда из этого пакета в том же самом формате, сжатый строковый формат, как иногда говорят, - альфа на тау на вот этот вот якобиан. В результате получается матрица А, которая хранится в
101:48
Speaker A
разреженном формате данных. То есть не nx - 1 на ny - 1 и всё в квадрате элементов, а гораздо меньше в памяти тратится места на хранение этой матрицы.
101:59
Speaker A
И быстрее реализуются все операции, в том числе и решение системы линейных алгебраических уравнений, учитывая то, что не делаются операции с нулевыми элементами, с которыми они должны делаться. А теперь давайте сейчас посмотрим на видео.
102:17
Speaker A
А, первое видео. Здесь я запустил для гораздо, для более густой сетки 100 на 100 и а, эпсилон был более высокого порядка малости. Вот один вариант этого видео и сейчас давайте посмотрим на другой вариант видео, также по аналогии, как с
102:34
Speaker A
первой частью, где у нас двухмерная картинка движется. Вот, в принципе, всё получается достаточно красиво. Поэтому, коллеги, на сегодня закругляемся и 2 минуты подытожу, что мы сегодня с вами рассмотрели. А мы с вами рассмотрели базовые подходы к решению многомерных
102:53
Speaker A
уравнений. Мы взяли задачу одномерную из лекции 19, обобщили её на двухмерный случай, решили её сначала с помощью явной схемы, которая и выписывается и реализуется жутко просто по сравнению с вот этим подходом.
103:11
Speaker A
И после этого мы что сделали? Мы уже сделали всё честно. Свели нашу задачу с помощью метода прямых к большой системе обыкновенных дифференциальных уравнений.
103:25
Speaker A
Эта система имеет специфическую структуру. Если мы её решаем неявными методами из, с помощью какого-то метода из семейства схем одностадийных Розенброка у нас возникает необходимость в якобиане. Якобиан получается разреженной матрицей очень специфической структуры, которую выгодно хранить в виде матрицы в разреженном формате
103:51
Speaker A
данных, что существенным образом ускоряет все и вычисления, и плюс даёт возможность оперировать на больших сетках, потому что хранимые матрицы имеют достаточно приемлемый размер. И по коллеги, финальное. В связи с тем, что соответствующие выкладки в этой лекции у нас были очень громоздкие, и что не то,
104:14
Speaker A
что выкладки я не могу все сделать даже с вами, как при во всех предыдущих лекциях я делал, а ещё до кучи я даже не могу нормальные презентации показать, потому что формул чрезвычайно много.
104:25
Speaker A
Напоминаю, в разделе материала лекции соответствующая программа приложена, и если будут вопросы, я в комментариях без проблем на все вопросы отвечу. А на сегодня коллеги всё.
Topics:численные методыуравнения в частных производныхмногомерные уравненияуравнение Бюргерсаграничные условияначальные условиячисленные схемысхемы Розенброкааппроксимация производныхустойчивость численных методов

Get More with the SozAI App

Transcribe recordings, audio files, and YouTube videos — with AI summaries and speaker detection. 30 minutes free.

Or transcribe another YouTube video here →