Матрица вторых производных трехмерного волнового уравнения

Матрица вторых производных трехмерного волнового уравнения

Глава 5. Решение дифференциальных уравнений

5.8 Функции решения параболических и гиперболических уравнений

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

Уравнения в частных производных можно разделить на три типа:

1) параболические, содержащие первую производную по одной переменной и вторую – по другой, причем все производные входят в уравнение с одинаковым знаком;

2) гиперболические, содержащие первую производную по одной переменной и вторую – по другой, входящие в уравнение с разными знаками;

3) эллиптические, содержащие вторые производные, причем одного знака.

Функции решения параболических и гиперболических уравнений

MathCAD 11 включает в себя две функции для решения параболических и гиперболических уравнений, pdesolve и numol .

Функция pdesolve используется в составе вычислительного блока Given – pdesolve для решения параболических и гиперболических уравнений (или систем уравнений) в частных производных, имеющих в качестве аргументов, как правило, время t и пространственную координату x .

Обращение к этой функции:

возвращает скалярную (для одного уравнения) или векторную (для системы уравнений) функцию, являющуюся решением уравнения (или системы уравнений). Здесь u –явно заданный вектор имен функций (без указания имен аргументов), подлежащих вычислению. Эти функции, а также граничные условия должны быть определены внутри вычислительного блока Given – pdesolve ; х – пространственная координата; x range – вектор значений аргумента х для граничных условий. Он должен состоять из двух чисел, представляющих две границы расчетного интервала; t – время (имя второго аргумента неизвестной функции); t range – вектор значений аргумента t для граничных условий, состоящих из двух чисел, представляющих две границы расчетного интервала; x pts – количество пространственных точек дискретизации (может не указываться); t pts – количество временных слоев (может не указываться).

Пример использования функции pdesolve приведен на рис. 5.19 запись вычислительного блока с функцией pdesolve аналогична записи блока с функцией Odesolve . Результаты расчета показаны на рис. 5.20.

Решение одномерного волнового уравнения

Здесь w -перемещение, v -скорость перемещения

где тогда

Представим первое уравнение как систему двух

уравнений первого порядка

Given

граничные условия

Рис. 5.1 9 Пример использования функции pdesolve

Единичное решение волнового уравнения

Сетка решений волнового уравнения на временном и пространственном интервалах

Рис. 5. 20 Результаты решения волнового уравнения

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

Недостатком функции pdesolve (как и функции Odesolve ) является невозможность ее использования в составе выражения – программы для многократного решения дифференциального уравнения. При необходимости многократного решения обыкновенных дифференциальных уравнений в состав программного модуля можно включать функции Rkadapt или Bulstoer.

При необходимости многократного решения дифференциальных уравнений в частных производных в состав программного модуля можно включать функцию numol , которая, как и pdesolve , появилась в MathCAD 11.

Функция numol предназначена для решения тех же уравнений, что и функция pdesolve .

Обращение к этой функции:

возвращает матрицу решения дифференциального уравнения в частных производных в каждой точке по пространственной (по строкам) и временной (по столбцам) координате. Если решается не одно уравнение, а система уравнений, то результатом решения является составная матрица, образованная путем слияния (слева направо) со значениями каждой искомой сеточной функции. Здесь x range – вектор значений аргумента х для граничных условий. Он должен состоять из двух чисел, представляющих две границы расчетного интервала; t range – вектор значений аргумента t для граничных условий, состоящих из двух чисел, педставляющих две границы расчетного интервала; x pts – количество пространственных точек дискретизации (может не указываться); t pts – количество временных слоев (может не указываться); N pde – количество дифференциальных уравнений в частных производных в системе; N ae – количество дополнительных алгебраических уравнений, входящих в систему; rhs – вектор правых частей уравнений; init – векторная функция, определяющая начальные условия для каждой неизвестной функции; bc – функциональная матрица граничных условий.

Вектор граничных условий может иметь значения трех типов:

– rhs содержит вторые пространственные производные: граничные условия (или Дирихле « D », или Неймана « N ») требуются по одному с каждой стороны интервала интегрирования;

– rhs содржит первые пространственные производные: граничные условия Дирихле на левой или правой границе интервала, на другой стороне NA ;

– нет пространственных производных – граничные условия не требуются.

Функциональная матрица bc содержит три столбца, имеющих ледующий вид:

– ( init _ left ( t ) init _ right ( t ) » D «) – для граничных условий Дирихле;

– ( init _ left ( t ) init _ right ( t ) » N «) – для граничных условий Неймана.

Пользоваться функцией numol намного сложнее, чем функцией pdesolve .

Граничное условие называется условием Дирихле, если задано значение функции, или Неймана, если задана первая производная функции.

Трехмерное волновое уравнение

Пусть мы по-прежнему имеем дело с плоской волной. Повернем координатные оси так, чтобы направление распространения волны задавалось каким-то единичным вектором n. Решение, очевидно, имеет вид:

(3.65)

Соотношения между w, kи v остаются прежними.

Волновой вектор — это вектор, модуль которого равен волновому числу, а направление совпадает с направлением распространения волны:

Фронт волны – плоскость, ортогональная волновому вектору k, – движется со скоростью v, оставаясь параллельным самому себе.

Найдем уравнение, которому удовлетворяет решение (3.65). Дважды дифференцируем выражение (3.65)по координатам х, у, z:

(3.66)

Складывая эти три уравнения, находим:

(3.67)

Вторая производная решения по времени дает уравнение:

(3.68)

Учитывая соотношение

получаем из (3.67), (3.68):

(3.69)

Выражение в скобках в левой части уравнения является дифференциальным оператором, который называется лапласианом и имеет специальное обозначение D

Записываем волновое уравнение для волн в трехмерном пространстве в окончательной форме:

(3.70)

Если волновая функция и зависит только от одной координаты (скажем, х), то лапласиан превращается во вторую производную по x и мы возвращаемся к прежней форме волнового уравнения.

Подчеркнем, что D не есть греческая буква D («дельта»), а Du не есть приращение величины u, но сумма вторых ее производных по координатам.

Но волновое уравнение (3.70) имеет и другие решения, нежели плоские волны. Простым дифференцированием можно убедиться, что сферическаяволна

(3.71)

удовлетворяет волновому уравнению. Фронт волны является сферой с центром в месте расположения источника колебаний (r=0), причем радиус сферы увеличивается со скоростью v.

Действительно, поверхность постоянной фазы дается уравнением

дифференцируя которое, находим

Амплитуда сферической волны

убывает с увеличением расстояния до точки наблюдения. Интенсивность волны

убывает по закону обратных квадратов. Это, как и закон Кулона, также связано с трехмерностью нашего пространства. Если среда не поглощает излучение, то поток энергии через поверхность сферы одинаков для сфер любых радиусов, окружающих источник излучения. Поскольку площадь сферы равна 4pr 2 , то энергия, проходящая через единицу площади, обратно пропорциональна r 2 .

Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. Студалл.Орг (0.004 сек.)

Матрица вторых производных трехмерного волнового уравнения

11.3.1. Параболические и гиперболические уравнения

Разработчики впервые применили дополнительные встроенные функции для решения параболических и гиперболических уравнений в частных производных в версии Mathcad 11, отлично осознавая значимость этих задач для современного исследователя и инженера. Предусмотрены два варианта решения: при помощи вычислительного блока Given/pdeso l ve , а также при помощи встроенной функции numo l . Первый путь проще в применении и нагляднее, зато второй позволяет автоматизировать процесс решения уравнений в частных производных, например, если нужно включить его в качестве составного шага в более сложную Mathcad-программу.

Вычислительный блок Given /pdesolve

Встроенная функция pdesolve применяется в рамках вычислительного блока, начинающегося ключевым словом Given , и пригодна для решения различных гиперболических и параболических уравнений. Она предназначена для решения одномерного уравнения (или системы уравнений) в частных производных (того, которое определит пользователь в рамках вычислительного блока Given ), зависящего от времени t и пространственной координаты х , имеет целый набор различных аргументов и работает следующим образом:

• pdesolve(u,x,xrange,t,trange,[xpts],[tpts])) — Возвращает скалярную (для единственного исходного уравнения) или векторную (для системы уравнений) функцию двух аргументов (x,t) , являющуюся решением дифференциального уравнения (или системы уравнений) в частных производных. Результирующая функция получается интерполяцией сеточной функции, вычисляемой согласно разностной схеме:

• u — явно заданный вектор имен функций (без указания имен аргументов), подлежащих вычислению. Эти функции, а также граничные условия (в форме Дирихле или Неймана) должны быть определены пользователем перед применением функции pdeso l ve в вычислительном блоке после ключевого слова Given . Если решается не система уравнений в частных производных, а единственное уравнение, то, соответственно, вектор и должен содержать только одно имя функции и вырождается в скаляр;
• х — пространственная координата (имя аргумента неизвестной функции);
• xrange — пространственный интервал, т. е. вектор значений аргумента х для граничных условий. Этот вектор должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала);
• t — время (имя аргумента неизвестной функции);
• t range — расчетная временная область: вектор значений аргумента t, который должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала по времени);
• xpts — количество пространственных точек дискретизации (может не указываться явно, в таком случае будет подобрано программой автоматически);
• tpts — количество временных слоев, т. е. интервалов дискретизации по времени (также может не указываться пользователем явно).

В качестве примера использования функции pdeso l ve (листинг 11.4) используем то же самое одномерное уравнение теплопроводности (11.5) с граничными и начальными условиями (11.6) и (11.7).

Листинг 11.4. Решение одномерного уравнения теплопроводности

Для корректного использования функции pdeso l ve предварительно, после ключевого слова Given , следует записать само уравнение и граничные условия при помощи логических операторов (для их ввода в Mathcad существует специальная панель). Обратите внимание, что уравнение должно содержать имя неизвестной функции u(x,t) вместе с именами аргументов (а не так, как она записывается в пределах встроенной функции pdeso l ve ). Для идентификации частных производных в пределах вычислительного блока следует использовать нижние индексы, например, uxx(,t) , для обозначения второй производной функции и по пространственной координате х .

Как видно из рис. 11.14, на котором изображены результаты расчетов по листингу 11.4, встроенная функция с успехом справляется с уравнением диффузии, отыскивая уже хорошо знакомое нам решение. Заметим, что использование встроенной функции pdeso l ve связано с довольно громоздкими вычислениями, которые могут отнимать существенное время.

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

Рис. 11.14. Решение уравнения диффузии тепла при помощи встроенной функции pdesoдve (листинг 11.4)

Пример: волновое уравнение

Приведем еще один пример применения функции pdesoive для решения уравнений в частных производных. Рассмотрим одномерное волновое уравнение, которое описывает, например, свободные колебания струны музыкального инструмента:

Здесь неизвестная функция u(x,t) описывает динамику смещения профиля струны относительно невозмущенного (прямолинейного) положения, а параметр с характеризует материал, из которого изготовлена струна.

Как вы видите, уравнение (11.11) содержит производные второго порядка, как по пространственной координате, так и по времени. Для того чтобы можно было использовать встроенную функцию pdesolve , необходимо переписать волновое уравнение в виде системы двух уравнений в частных производных, введя вторую неизвестную функцию v=ut . Программа для решения волнового уравнения приведена в листинге 11.5, а результат— на рис. 11.15.

Листинг 11.5. Решение волнового уравнения

Рис. 11.15. Решение волнового уравнения (продолжение листинга 11.5)

Встроенная функция numol

Альтернативный вариант решения дифференциальных уравнений в частных производных заключается в применении еще одной встроенной функции numo |, которая реализует тот же самый алгоритм сеток, позволяя вручную задать большинство его параметров:

• numol(xrange,xpts,trange,tpts,Npde,Nae,rhs,init,bc) — Возвращает матрицу решения дифференциального уравнения в частных производных, представляющую искомую сеточную функцию в каждой точке по пространственной (по строкам) и временной координате (по столбцам). Если решается не одно уравнение, а система уравнений, то результатом является составная матрица, образованная путем слияния (слева-направо) матриц со значениями каждой искомой сеточной функции:

• Npde — общее количество дифференциальных уравнений в частных производных в системе;
• Nae — общее количество дополнительных алгебраических уравнений, которые также могут входить в систему;
• rhs — векторная функция, определяющая систему дифференциальных и алгебраических уравнений (формат этого и двух следующих матричных параметров объяснен в листинге 11.9);
• init — векторная функция, определяющая начальные условия для каждой неизвестной функции;
• be — функциональная матрица граничных условий;
• xrange — пространственный интервал, т. е. вектор значений аргумента х для граничных условий. Этот вектор должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала);
• xpts — количество пространственных точек дискретизации (может не указываться явно, в таком случае будет подобрано программой автоматически);
• trange — расчетная временная область: вектор значений аргумента t, который должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала по времени);
• tpts — количество временных слоев, т. е. интервалов дискретизации по времени (также может не указываться пользователем явно);

Пример решения волнового уравнения при помощи функции numol приведен в листинге 11.6, особое внимание в котором мы призываем уделить формату представления векторов rhs , init и be , а также принципу извлечения отдельных сеточных решений из матрицы-результата. График решения, показанный на рис. 11.16, полезно сравнить с результатом применения вычислительного блока из предыдущего раздела (см. листинг 11.5 и рис. 11.15).

Листинг 11.6 . Решение волнового уравнения при помощи функции numol

Как вы видите, функция numol имеет еще большее число аргументов, нежели pdesolve , и позволяет автоматизировать применение метода сеток. Однако пользоваться ею намного сложнее, чем вычислительным блоком, поскольку и уравнения, и начальные и граничные условия должны быть записаны в специальном формате. Применение функции numol оправданно, когда необходимо включить решение уравнений в частных производных в более сложные вычисления в качестве подпрограммы, организовать серию расчетов с меняющимся параметром, подготовить анимацию графиков решения и т. п.

Рис. 11.16. Решение волнового уравнения (продолжение листинга 11.6)

Именно в целях визуализации решения параболических и гиперболических уравнений в частных производных использование функции numol наиболее полезно. График решения динамических уравнений (зависящих от времени t) выглядит намного эффектнее и воспринимается несравненно лучше, если он оформлен в виде анимации. Для создания анимационных роликов расчетное время следует выразить через константу FRAME и затем применить команду View / Animate (Вид / Анимация) (см. разд. 13.3.2).