Разностном методе

10.4.1.0 разностном методе

Разберем идею разностного метода решения краевых задач на примере взаимодействия световых пучков (см. Рисунок 10.1), переобозначив в системе (10.1) интенсивность излучения вправо на Y, а интенсивность излучения влево на у (просто в целях удобства, чтобы не писать в дальнейшем индекс). Суть метода заключается в покрытии расчетного интервала сеткой из N точек. Тем самым определяются (м-l) шагов (Рисунок 10.7). Затем надо заменить дифференциальные уравнения исходной краевой задачи аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого 1-го шага. В нашем случае достаточно просто заменить первые производные из (10.1) их разностными аналогами (такой метод называется еще методом Эйлера):

Алгоритм стрельбы

10.2.1. Алгоритм стрельбы

Суть метода стрельбы заключается в пробном задании недостающих граничных условий на левой границе интервала и решении затем полученной задачи Коши хорошо известными методами (см. главу 11). В нашем примере не хватает начального условия для y1(0), поэтому сначала зададим ему произвольное значение, например, у1(0)=10. Конечно, такой выбор не совсем случаен, поскольку из физических соображений ясно, что, во-первых, интенсивность излучения — величина заведомо положительная, и, во-вторых, отраженное излучение должно быть намного меньше падающего. Решение задачи Коши с помощью функции rkfixed приведено в листинге 10.1, причем в его последней строке выведены расчетные значения у0 и уг на правой крайней точке интервала х=1. Разумеется, они не совпадают, т. е. правое краевое условие не выполнено, из чего следует, что полученный результат не является решением поставленной краевой задачи.

О постановке задачи

10.5.1. О постановке задачи

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

a(x,Y,y)=a0(x)-?(x)(Y+y), (10.8)

где а0(х) — постоянная (не зависящая от интенсивности света) составляющая коэффициента поглощения, а е(х) — малый коэффициент пропорциональности. В этом случае модель (10.1) перепишется в виде:

Таким образом, физическая подоплека квадратичной нелинейности задачи (10.9) связана с зависимостью коэффициента поглощения среды от Y+y. Иными словами, чем мощнее свет в некоторой точке среды, тем меньше локальный коэффициент поглощения среды в этой точке.

Механизм этого явления будет хорошо понятен, если предположить, что средой является воздух, насыщенный парами жидкости (попросту говоря, туман). Чем сильнее в некоторой точке туман, тем больше локальное поглощение. Если свет, проходящий через туман, имеет малую интенсивность, то он не оказывает на среду никакого (или почти никакого) воздействия. Однако если пучок света будет очень мощным, за счет его поглощения пары воды будут разогреваться и, следовательно, испаряться. По мере испарения капель воды туман станет рассеиваться, и среда просветлеет (т. е. ее коэффициент поглощения уменьшится). Таким образом, чем мощнее свет, тем активнее пойдет процесс испарения, тем существеннее уменьшится коэффициент поглощения. Именно такой эффект, совершенно понятный с физической точки зрения, и описывает формула (10.8) и, соответственно, модель (10.9).

Подчеркнем еще раз и физический смысл коэффициентов а0(х) и ?(х). Первый описывает постоянный (не зависящий от интенсивности света) коэффициент поглощения среды, а второй — его зависимость от эффекта разогрева среды светом. Чем больше коэффициент а0(х), тем более жесткой является краевая задача; а чем больше ?(х), тем сильнее ее нелинейность. В предельном случае ?(х)->0 задача (10.9) переходит в задачу (10.1), становясь линейной.

Двухточечные краевые задачи

10.2.2. Двухточечные краевые задачи

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

sbval(z,x0,x1, D, load, score) — поиск вектора недостающих L начальных условий для двухточечной краевой задачи для системы N ОДУ:

 z — вектор размера Lx1, присваивающий недостающим начальным условиям (на левой границе интервала) начальные значения;  х0 — левая граница расчетного интервала;  xi — правая граница расчетного интервала;  load(x0,z) — векторная функция размера Nx1 левых граничных условий, причем недостающие начальные условия поименовываются соответствующими компонентами векторного аргумента z;  score (x1, у) — векторная функция размера Lx1, выражающая L правых граничных условий для векторной функции у в точке x1;  D(x,y) — векторная функция, описывающая систему N ОДУ, размера Nx1 и двух аргументов — скалярного х и векторного у. При этом у — это неизвестная векторная функция аргумента х того же размера Nx1.

ВНИМАНИЕ!

Решение краевых задач в Mathcad, в отличие от большинства остальных операций, реализовано не совсем очевидным образом. В частности, помните, что число элементов векторов о и load равно количеству уравнений N, а векторов z, score и результата действия функции sbval— количеству правых граничных условий L. Соответственно, левых граничных условий в задаче должно быть (N-L) .

Как видно, функция sbval предназначена не для поиска собственно решения, т. е. неизвестных функций yi(x), а для определения недостающих начальных условий в первой точке интервала, т. е. yi(x0). Чтобы вычислить yi(х) на всем интервале, требуется дополнительно решить задачу Коши. Разберем особенности использования функции sbval на конкретном примере (листинг 10.2), описанном выше (см. разд. 10.1.1). Краевая задача состоит из системы двух уравнений (N=2), одного левого (L=1) и одного правого (N-L=2-1=1) граничного условия.

Метод стрельбы

10.5.2. Метод стрельбы

Встроенная функция sbval, реализующая в Mathcad алгоритм стрельбы (см. разд. 10.2), может справляться и с нелинейными задачами. Приведем пример решения краевой задачи (10.9) (листинг 10.8 и Рисунок 10.11) с теми же граничными условиями, что были поставлены для модели (10.1). Интерес представляет решение у, описывающее рост интенсивности отраженного пучка по мере его распространения справа-налево (обратите внимание, что на Рисунок 10.11 оно отложено с десятикратным увеличением).

Жесткие краевые задачи

10.4.2. Жесткие краевые задачи

Один из случаев, когда применение разностных схем может быть очень полезным, связан с решением жестких краевых задач (подробнее о жестких ОДУ читайте в главе 11). В частности, рассматриваемая задача о встречных световых пучках становится жесткой при увеличении коэффициента ослабления а(х) в несколько десятков раз. Например, при попытке решить ее с а(х) :=100 с помощью листинга 10.2 вместо ответа выдается сообщение об ошибке "Can't converge to a solution. Encountered too many integration steps" (He сходится к решению. Слишком много шагов интегрирования). Это и неудивительно, поскольку жесткие системы характерны тем, что требуют исключительно малого значения шага в стандартных алгоритмах.

Для жестких задач неприменимы и явные разностные схемы, о которых рассказывалось в предыдущем разделе (см. разд. 10.3.1). Результат расчетов по программе листинга 10.6, например с а(х) :=20 (Рисунок 10.9), дает характерную для неустойчивых разностных схем "разболтку" — колебания нарастающей амплитуды, не имеющие ничего общего с реальным решением.

Краевые задачи с условием во внутренней точке

10.2.3. Краевые задачи с условием во внутренней точке

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

bvaifit(z1,z2,x0,x1,xf,D,load1,load2,score) — поиск вектора недостающих граничных условий для краевой задачи с дополнительным условием в промежуточной точке для системы N ОДУ.

 z1 — вектор, присваивающий недостающим начальным условиям на левой границе интервала начальные значения.  z2 — вектор того же размера, присваивающий недостающим начальным условиям на правой границе интервала начальные значения.  х0 — левая граница расчетного интервала.  x1 — правая граница расчетного интервала.  xf — точка внутри интервала. Dо(х,у) — векторная функция, описывающая систему N ОДУ размера Nx1 и двух аргументов — скалярного х и векторного у. При этом у — это неизвестная векторная функция аргумента х того же размера Nx1.  load1(x0,z) — векторная функция размера Nx1 левых граничных условий, причем недостающие начальные условия поименовываются соответствующими компонентами векторного аргумента z.  load2(x1,z) — векторная функция размера Nx1 левых граничных условий, причем недостающие начальные условия поименовываются соответствующими компонентами векторного аргумента z.  score (xf, у) — векторная функция размера Nx1, выражающая внутреннее условие для векторной функции у в точке xf.

Обычно функция bvaifit применяется для задач, в которых производная решения имеет разрыв во внутренней точке xf. Некоторые из таких задач не удается решить обычным методом пристрелки, поэтому приходится вести пристрелку сразу из обеих граничных точек. В этом случае дополнительное внутреннее условие в точке xf является просто условием сшивки в ней левого и правого решений. Поэтому для данной постановки задачи оно записывается В форме score (xf, у) :=у.

Рассмотрим действие функции bvaifit на знакомом примере модели взаимодействия пучков света (см. Рисунок 10.1), предположив, что в промежутке между xf=0.5 и x1=1.0 находится другая, оптически более плотная среда с другим коэффициентом ослабления излучения а(х)=3. Соответствующая краевая задача решена в листинге 10.3, причем разрывный показатель ослабления определяется в его второй строке.

Разностные схемы

10.5.3. Разностные схемы

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

Подробная разработка алгоритмов и соответствующих программ Mathcad выходит далеко за пределы данной книги, поэтому мы конспективно представим один из приемов решения нелинейных краевых задач, сводящийся к их линеаризации. В общих чертах подход заключается в следующем. Предположим, что некоторое приближение (обозначим его Y(х) и у(х)) к решению нелинейной задачи (10.9) нам известно, и можно считать, что Y->Y+Z и y->y+z, где z и z — близкие к нулю функции х. Тогда, пользуясь их малостью (по сравнению с Y0 и у0), можно разложить нелинейные слагаемые в уравнениях (10.9) в ряд Тейлора по z и z. Получим:

Y'+Z' = -aY + ry + ?(Y2 + Yy) - aZ + rz + 2?Y(Z + z); (10.10)

y'+z' = ay - rY - ?(y2 + Yy) + az - rz - 2?y(Z + z) .

Теперь, поскольку Y(x) и у(х) являются приближением к решению исходной задачи, то можно считать, что они (приблизительно) удовлетворяют и уравнению, и граничным условиям. Тогда, вычитая (10.9) из (10.10), получим краевую задачу для z (х) и z (х):

Z' = -aZ + rz + 2?Y(Z + z);

z' = az - rZ - 2?y(Z + z) ; (10.11)

Z(0) = z(0) = Z(l) = z(l) = 0.

Это и есть та самая новая задача для "маленьких" функций z (х) и z (х), которую надо решить, и которая, благодаря малости z и z, является линейной. Вся беда в том, что мы не знаем "больших" функций Y(x) и у(х), а они, увы, входят в задачу (10.11). Тем не менее рецепт получения этих функций напрашивается сам собой: если нелинейность исходных ОДУ не слишком сильная, то в качестве Y и у можно взять решение линейной краевой задачи, т. е. задачи (10.9) с ?(х)=0 (см. разд. 10.4.1).

Сказанное иллюстрируют листинги 10.9 и 10.10, первый из которых решает линейную краевую задачу (являясь, фактически, небольшой модификацией листинга 10.7 из разд. 10.4.2), а второй решает линеаризованную задачу (10.11), учитывая результат листинга 10.9. В листинге 10.9 матрица о и вектор правых частей в являются разностной аппроксимацией ОДУ (ее первые N строк аппроксимируют первое уравнение, а оставшиеся N строк — второе). Такой же смысл и точно такую же структуру имеют матрица с и вектор правых частей с для второй (линеаризованной) задачи (10.11). Для решения сеточных уравнений Dу=B и Cz=G используется (конечно, весьма неэкономично) встроенная функция isolve, реализующая алгоритм Гаусса.

Важно привлечь внимание читателя к последним строкам листинга 10.9. В них осуществляется интерполяция полученного решения системы сеточных уравнений для того, чтобы в нелинейной задаче (в листинге 10.10) можно было использовать непрерывные "большие" функции из линейной задачи. В последней строке листинга 10.10 осуществляется сложение "больших" и "маленьких" функций (результатов листингов 10.9 и 10.10) для получения полного решения нелинейной задачи (10.9), которое изображено на Рисунок 10.12.

Не будем давать дополнительных комментариев к Mathcad-программам, надеясь, что читатель, заинтересовавшийся нелинейным примером со световыми пучками и эффектом разогрева светом среды, сам разберется в листингах, тем более что техника разработки разностных схем была нами детально разобрана раньше (см. разд. 10.4).

Обыкновенные дифференциальные уравнения краевые задачи

Глава 10. Обыкновенные дифференциальные уравнения: краевые задачи

Обыкновенные дифференциальные уравнения: краевые задачи 10.1. О постановке задач 10.2. Решение краевых задач средствами Mathcad 10.2.1. Алгоритм стрельбы 10.2.2. Двухточечные краевые задачи 10.2.3. Краевые задачи с условием во внутренней точке 10.3. Задачи на собственные значения для ОДУ 10.4. Разностные схемы для ОДУ 10.4.1. О разностном методе 10.4.2. Жесткие краевые задачи 10.5. Нелинейные краевые задачи 10.5.1. О постановке задачи 10.5.2. Метод стрельбы 10.5.3. Разностные схемы

Содержание

Иллюстрация метода стрельбы (продолжение листинга 10 1)

Рисунок 10.2. Иллюстрация метода стрельбы (продолжение листинга 10.1)

Решение пробной задачи Коши для модели (10 1)

Листинг 10.1. Решение пробной задачи Коши для модели (10.1)

График полученных решений показан на Рисунок 10.2 (слева). Из него также видно, что взятое наугад второе начальное условие не обеспечило выполнение граничного условия при х=1. В целях лучшего выполнения этого граничного условия следует взять большее значение y1(0), например, y1(0)=15, и вновь решить задачу Коши. Соответствующий результат показан на том же Рисунок 10.2 (в центре). Граничное условие выполняется с лучшей точностью, но опять-таки оказалось недостаточным. Для еще одного значения y1(0)=20 получается решение, показанное на Рисунок 10.2 (справа). Из сравнения двух правых графиков легко заключить, что недостающее начальное условие больше 15, но меньше 20. Продолжая подобным образом "пристрелку" по недостающему начальному условию, возможно отыскать правильное решение краевой задачи.

В этом и состоит принцип алгоритма стрельбы. Выбирая пробные начальные условия (проводя пристрелку) и решая соответствующую серию задач Коши, можно найти то решение системы ОДУ, которое (с заданной точностью) удовлетворит граничному условию (или, в общем случае, условиям) на другой границе расчетного интервала.

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

Решение краевой задачи

Листинг 10.2. Решение краевой задачи

Первые три строки листинга задают необходимые параметры задачи и саму систему ОДУ. В четвертой строке определяется вектор z. Поскольку правое граничное условие всего одно, то недостающее начальное условие тоже одно, соответственно, и вектор z имеет только один элемент z0. Ему необходимо присвоить начальное значение (мы приняли z0=10, как в листинге 10.1), чтобы запустить алгоритм стрельбы (см. разд. 10.1.2).

Примечание 1
Примечание 1

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

В следующей строке листинга векторной функции load(x,z) присваиваются левые граничные условия. Эта функция аналогична векторной переменной, определяющей начальные условия для встроенных функций, решающих задачи Коши. Отличие заключается в записи недостающих условий. Вместо конкретных чисел на соответствующих местах пишутся имена искомых элементов вектора z. В нашем случае вместо второго начального условия стоит аргумент z0 функции load. Первый аргумент функции load — это точка, в которой ставится левое граничное условие. Ее конкретное значение определяется непосредственно в списке аргументов функции sbval. Следующая строка листинга определяет правое граничное условие, для введения которого используется функция score(х,у). Оно записывается точно так же, как система уравнений в функции о. Аргумент х функции score аналогичен функции load и нужен для тех случаев, когда граничное условие явно зависит от координаты х. Вектор score должен состоять из такого же числа элементов, что и вектор z.

Реализованный в функции sbval алгоритм стрельбы ищет недостающие начальные условия таким образом, чтобы решение полученной задачи Коши делало функцию score (х, у) как можно ближе к нулю. Как видно из листинга, результат применения sbval для интервала (0,1) присваивается векторной переменной и. Этот вектор похож на вектор z, только в нем содержатся искомые начальные условия вместо приближенных начальных значений, заданных в z. Вектор и содержит, как и z, всего один элемент и0. С его помощью можно определить решение краевой задачи у (х) (последняя строка листинга). Тем самым функция sbval сводит решение краевых задач к задачам Коши. График решения краевой задачи показан на Рисунок 10.3.

На Рисунок 10.4 показано решение той же самой краевой задачи, но с другим правым граничным условием, соответствующим R=0, т. е. без зеркала на правой границе. В этом случае слабый обратный пучок света образуется исключительно за счет обратного рассеяния излучения от лазера. Конечно, многие из читателей уже обратили внимание, что реальная физическая среда не может создавать такого большого рассеяния назад. Иными словами, более реальны значения r(х)<<а (х). Однако когда коэффициенты в системе ОДУ при разных yi очень сильно (на порядки) различаются, система ОДУ становится жесткой, и функция sbval не может найти решения, выдавая вместо него сообщение об ошибке "Could not find a solution" (Невозможно найти решение).

ВНИМАНИЕ!

Метод стрельбы не годится для решения жестких краевых задач. Поэтому алгоритмы решения жестких ОДУ в Mathcad приходится программировать самому (см. разд. 10.4).

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

Листинг 10.3. Краевая задача с граничным условием в промежуточной точке

Система уравнений и левое краевое условие вводится так же, как и в предыдущем листинге для функции sbval. Обратите внимание, что таким же образом записано и правое краевое условие. Для того чтобы ввести условие отражения на правой границе, пришлось определить еще один неизвестный пристрелочный параметр z20. Строка листинга, в которой определена функция score, задает условие стрельбы — сшивку двух решений в точке xf. В самой последней строке листинга выдан ответ — определенные численным методом значения обоих пристрелочных параметров, которые объединены в вектор а (мы применили в предпоследней строке операцию транспонирования, чтобы результат получился в форме вектора, а не матрицы-строки). Для корректного построения графика решения лучше составить его из двух частей — решения задачи Коши на интервале (x0,xf) и другой задачи Коши для интервала (xf,x1). Реализация этого способа приведена в листинге 10.4, который является продолжением листинга 10.3. В последней строке листинга 10.4 выведено значение второй искомой функции на правой границе интервала. Всегда полезно проконтролировать, что оно совпадает с соответствующим пристрелочным параметром (выведенным в последней строке листинга 10.3).

Решение краевой задачи (продолжение листинга 10 3)

Листинг 10.4. Решение краевой задачи (продолжение листинга 10.3)

Решение краевой задачи приведено на Рисунок 10.5. С физической точки зрения естественно, что интенсивность света уменьшается быстрее по мере распространения в более плотной среде в правой половине расчетного интервала. В средней точке xf=0.5, как и ожидалось, производные обоих решений имеют разрыв.

Примечание 1
Примечание 1

Еще один пример решения краевых задач с разрывными коэффициентами ОДУ приведен в справочной системе Mathcad.

Решение задачи о собственных колебаниях струны

Листинг 10.5. Решение задачи о собственных колебаниях струны

В первых двух строках листинга определяются функции, входящие в задачу, в том числе р' (х) :=0, и границы расчетного интервала (0,1). В третьей строке дается начальное приближение к собственному значению ?0, в четвертой вводится система ОДУ. Обратите внимание, что она состоит не из двух, а из трех уравнений. Первые два из них определяют эквивалентную (10.3) систему ОДУ первого порядка, а третье необходимо для задания собственного значения в виде еще одного компонента у2 искомого вектора у. Поскольку по определению собственное значение постоянно при всех х, то его производная должна быть приравнена нулю, что отражено в последнем уравнении. Важно также, что во втором из уравнений собственное значение записано как у2, поскольку является одним из неизвестных.

В следующих двух строках листинга задается левое граничное условие, включающее и недостающее условие на собственное значение для третьего уравнения, и правое граничное условие у0=0. В предпоследней строке листинга обычным образом применяется функция sbval, а в последней выводится результат ее работы вместе с известным аналитически собственным значением п2-я2. Как легко убедиться, мы нашли первое собственное значение для n=1, а чтобы найти другие собственные значения, необходимо задать другие начальные приближения к ним (в третьей строке листинга 10.5). Например, выбор ?0=50 приводит ко второму собственному значению 22?2, а ?0=100 — к третьему 32?2.

Чтобы построить график соответствующей собственной функции, надо добавить в листинг строку, программирующую решение задачи Коши, например, такую: u:=rkfixed(load(a,A) ,а,b, 100,D) . Полученные кривые показаны на Рисунок 10.6 в виде коллажа трех графиков, рассчитанных для трех собственных значений.

Примечание 1
Примечание 1

Примеры решения нескольких задач на собственные значения можно найти в различных Ресурсах Mathcad.

Реализация явной разностной схемы

Листинг 10.6. Реализация явной разностной схемы

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

В первой строке листинга определяются функции и константы, входящие в модель, во второй задается число точек сетки N=5 и ее равномерный шаг. Следующие две строки определяют матричные коэффициенты, аппроксимирующие уравнения для Yi; а пятая и шестая — для уi. Седьмая и восьмая строки листинга задают соответственно левое и правое граничное условие, а строки с девятой по одиннадцатую — правые части системы (10.6). В следующей строке завершается построение матрицы А вырезанием из нее левого нулевого столбца. В предпоследней строке листинга применена встроенная функция isolve для решения системы (10.6), а в последней выведены рассчитанные ею неизвестные граничные значения. Графики решения приведены на Рисунок 10.8, причем первые N элементов итогового вектора есть вычисленное излучение вперед, а последние N элементов — излучение назад.

Реализация неявной разностной схемы для жесткой краевой задачи

Листинг 10.7. Реализация неявной разностной схемы для жесткой краевой задачи

Не будем специально останавливаться на обсуждении листинга 10.7, поскольку он почти в точности повторяет предыдущий листинг. Отличие заключается лишь в формировании матрицы А другим способом, согласно неявной схеме. Решение, показанное на Рисунок 10.10, демонстрирует, что произошло "небольшое чудо": "разболтка" исчезла, а распределение интенсивностей стало физически предсказуемым. Обратите внимание, что (из-за взятого нами слишком большого коэффициента ослабления излучения) отраженный пучок света имеет очень маленькую интенсивность, и ее пришлось построить на графике с увеличением в тысячу раз.

Решение нелинейной краевой задачи методом стрельбы

Листинг 10.8. Решение нелинейной краевой задачи методом стрельбы

Следует помнить, что чем сильнее нелинейность и чем жестче краевая задача, тем более сильные требования предъявляются к начальному значению алгоритма, ввод которого осуществляется посредством функции load. Попробуйте повторить расчеты листинга 10.8 с иным сочетанием параметров а0(х) и ?(х), и вы увидите, что с немного более жесткими задачами алгоритм стрельбы уже не справляется.

Примечание 1
Примечание 1

На компакт-диске, прилагаемом к книге, вы найдете еще более интересное решение для другой задачи, с нелинейностью типа (y+Y)2.

Решение линейной (приближенной) краевой задачи

Листинг 10.9.Решение линейной (приближенной) краевой задачи

Решение линеаризованной задачи (продолжение листинга 10 9)

Листинг 10.10. Решение линеаризованной задачи (продолжение листинга 10.9)

="27.gif">

Последний важный момент, который следует обозначить, связан с решением задач, обладающих значительной нелинейностью. Решение, приведенное в листингах 10.9, 10.10 и на Рисунок 10.12, согласно самой постановке, должно не сильно отличаться от решения линейной краевой задачи, поскольку функции z(x) и z(x) малы по сравнению с Y(x) и у(х). Если же нелинейность сильная, то решение может сильно отличаться от Y и у, и линеаризация (10.11) будет просто неправильной. В этом случае следует слегка усложнить алгоритм решения нелинейной краевой задачи.

Модель краевой задачи

Рисунок 10.1. Модель краевой задачи

Полученную задачу называют краевой (boundary value problem), поскольку условия поставлены не на одной, а на обеих границах интервала (0,1). И, в связи с этим, их не решить методами предыдущей главы, предназначенными для задач с начальными условиями. Далее для показа возможностей Mathcad будем использовать этот пример с R=1 и конкретным видом a(x)=const=1 и r(x)=const=0.1, описывающим случай изотропного (не зависящего от координаты х) рассеяния.

Примечание 2
Примечание 2

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

Примечание 3
Примечание 3

Модель встречных световых пучков привела нас к системе уравнений (10.1), в которые входят производные только по одной переменной х. Если бы мы стали рассматривать более сложные эффекты рассеяния в стороны (а не только вперед и назад), то в уравнениях появились бы частные производные по другим пространственным переменным у и z. В этом случае получилась бы краевая задача для уравнений в частных производных, решение которой во много раз сложнее ОДУ.

Нелинейные краевые задачи

10.5. Нелинейные краевые задачи

Все примеры, приведенные пока в этом разделе, связаны с краевыми задачами для линейных ОДУ. Между тем на практике часто приходится исследовать именно нелинейные задачи, которые несравненно сложнее с вычислительной точки зрения. Рассмотрим оба подхода к решению нелинейных задач (алгоритм стрельбы и разностный метод), не забывая о том, что нелинейные ОДУ также могут быть в той или иной степени жесткими.

Неверное решение жесткой краевой задачи по неустойчивой явной разностной схеме

Рисунок 10.9. Неверное решение жесткой краевой задачи по неустойчивой явной разностной схеме

Выходом из положения будет использование неявных разностных схем. Применительно к нашей задаче достаточно заменить правые части уравнений (10.1) значениями не на левой, а на правой границе каждого шага:

Граничные условия, конечно, можно оставить в том же виде (10.5). Поскольку мы имеем дело с линейными дифференциальными уравнениями, то и схему (10.7) легко будет записать в виде матричного равенства (10.6), перегруппировывая соответствующим образом выражение (10.7) и приводя подобные слагаемые. Разумеется, полученная матрица А будет иной, нежели матрица А для явной схемы (10.4). Поэтому и решение (реализация неявной схемы) может отличаться от изображенного на Рисунок 10.9 результата расчетов по явной схеме. Программа, составленная для решения системы (10.7), приведена в листинге 10.7.

О постановке задач

10.1.О постановке задач

Постановка краевых задач для ОДУ отличается от задач Коши, рассмотренных в главе 9, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеется система N обыкновенных дифференциальных уравнений первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия — на противоположной границе.

Примечание 1
Примечание 1

Дифференциальные уравнения высших порядков можно свести к эквивалентной системе ОДУ первого порядка (см. главу 9).

Чтобы лучше понять, что из себя представляют краевые задачи, рассмотрим их постановочную часть на конкретном физическом примере модели взаимодействия встречных световых пучков. Предположим, что надо определить распределение интенсивности оптического излучения в пространстве между источником (лазером) и зеркалом, заполненном некоторой средой (Рисунок 10.1). Будем считать, что от зеркала отражается R-Я часть падающего излучения (т. е. его коэффициент отражения равен R), а среда как поглощает излучение с коэффициентом ослабления а(х), так и рассеивает его. Причем коэффициент рассеяния назад равен г(х). В этом случае закон изменения интенсивности у0(х) излучения, распространяющегося вправо, и интенсивности y1 (х) излучения влево определяется системой двух ОДУ первого порядка:

Для правильной постановки задачи требуется, помимо уравнений, задать такое же количество граничных условий. Одно из них будет выражать известную интенсивность излучения 10, падающего с левой границы х=0, а второе — закон отражения на его правой границе x=1:

y0(0) = 10; (10.2)

y1(l) = Rxy0(l).

Обыкновенные дифференциальные уравнения краевые задачи

Обыкновенные дифференциальные уравнения: краевые задачи

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

Краевые задачи во множестве практических приложений часто зависят от некоторого числового параметра. При этом решение существует не для всех его значений, а лишь для счетного их числа. Такие задачи называют задачами на собственные значения (см. разд. 10.3).

Несмотря на то, что, в отличие от задач Коши для ОДУ, в Mathcad не предусмотрены встроенные функции для решения жестких краевых задач, с ними все-таки можно справиться, применив программирование разностных схем, подходящих для решения задач этого класса (см. разд. 104). О подходе к решению нелинейных краевых задач написано в конце главы (см. разд. 10.5).

Первые три собственные функции

Рисунок 10.6. Первые три собственные функции задачи колебаний струны (коллаж трех графиков, рассчитанных листингом 10.5 для разных ?0)

Разностные схемы для ОДУ

10.4. Разностные схемы для ОДУ

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

Решение краевой задачи для R=0 (продолжение листинга 10 2 при R=0)

Рисунок 10.4. Решение краевой задачи для R=0 (продолжение листинга 10.2 при R=0)

Решение краевой задачи для R=l (продолжение листинга 10 2)

Рисунок 10.3. Решение краевой задачи для R=l (продолжение листинга 10.2)

Решение краевой задачи по неявной разностной схеме (продолжение листинга 10 7)

Рисунок 10.10. Решение краевой задачи по неявной разностной схеме (продолжение листинга 10.7)

Решение краевой задачи разностным методом (продолжение листинга 10 6)

Рисунок 10.8. Решение краевой задачи разностным методом (продолжение листинга 10.6)

Как мы увидели, реализация в Mathcad разностных схем вполне возможна и не слишком трудоемка — предложенная программа состоит всего из двух десятков математических выражений. Конечно, для их написания требуется и время, и часто кропотливые расчеты, но, собственно, в этом и состоит работа математика. Кстати говоря, при небольшом числе шагов расчеты по разностным схемам не требуют существенного времени (программа, приведенная в листинге 10.6, работает быстрее, чем метод стрельбы, встроенный в функцию sbval). Существуют, кроме того, весьма очевидные для многих читателей пути ускорения расчетов, связанные с применением более подходящих методов решения систем линейных уравнений с разреженной матрицей.

Решение краевой задачи с разрывом в xf=0 5 (продолжение листингов 10 3—10 4)

Рисунок 10.5. Решение краевой задачи с разрывом в xf=0.5 (продолжение листингов 10.3—10.4)

Ради справедливости необходимо заметить, что разобранную краевую задачу легко решить и с помощью функции sbval, заменив в листинге 10.2 зависимость а (х) на третью строку листинга 10.3. В этом случае листинг 10.2 даст в точности тот же ответ, что показан на Рисунок 10.5. Однако в определенных случаях (в том числе из соображений быстроты расчетов) удобнее использовать функцию bvaifit, т. е. вести пристрелку с обеих границ интервала.

СОВЕТ

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

Решение краевых задач средствами Mathcad

10.2. Решение краевых задач средствами Mathcad

Для решения краевых задач в Mathcad реализован наиболее популярный алгоритм, называемый методом стрельбы или пристрелки (shooting method). Он, по сути, сводит решение краевой задачи к решению серии задач Коши с различными начальными условиями. Рассмотрим сначала основной принцип алгоритма стрельбы на примере модели световых пучков (см. Рисунок 10.1), а встроенные функции, реализующие этот алгоритм, приведем в следующем разделе.

Решение нелинейной краевой задачи (продолжение листинга 10 8)

Рисунок 10.11. Решение нелинейной краевой задачи (продолжение листинга 10.8)

Решение нелинейной краевой задачи разностным методом (продолжение листингов 10 9 и 10 10)

Рисунок 10.12. Решение нелинейной краевой задачи разностным методом (продолжение листингов 10.9 и 10.10)

Обозначим полученное в результате решение, как и в листинге 10.10, вектором J(?), подчеркивая тем самым его зависимость от параметра нелинейности. Очевидно, что J(0) есть решение линейной задачи. Для того чтобы решить задачу с сильной нелинейностью, т. е. довольно большим ?=?1, можно организовать продолжение по ? как по параметру. Иными словами, используя в качестве начального приближения J (0), можно решить задачу для другого, малого ?=Д?, получив J(A?), затем, взяв это J(?) в качестве приближенного решения, получить J(2A?) и т. д. малыми шагами добраться до желаемого ?1.

Примечание 1
Примечание 1

Упрощенную реализацию этого алгоритма вы найдете на компакт-диске, прилагаемом к книге. Она связана с выводом во внешний файл данных результата задачи из листинга 10.10 и считыванием из него же этих данных в качестве входной информации для следующей итерации. В качестве нулевой итерации используется решение линейной задачи, выводимое предварительно в файл из усовершенствованного листинга 10.9.

Сетка покрывающая расчетный интервал

Рисунок 10.7. Сетка, покрывающая расчетный интервал

Примечание 1
Примечание 1

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

Получилась система (по числу шагов) 2(N-1) разностных линейных алгебраических уравнений с 2N неизвестными Yi и уi. Для того чтобы она имела единственное решение, надо дополнить число уравнений до 2N.

Это можно сделать, записав в разностном виде оба граничных условия: \

Y0=10, YN=RYN. (10.5)

Сформированная полная система алгебраических уравнений называется разностной схемой, аппроксимирующей исходную краевую задачу. Обратите внимание, что правые части разностных уравнений системы (10.4) на каждом шаге записаны для левой границы шага. Такие разностные схемы называют явными, т. к. все значения Yi+1 и yi+1 находятся в левой части уравнений. Полученную явную разностную схему легко записать в матричной форме

Аz=В, (10.6)

где z — неизвестный вектор, получающийся объединением векторов Y и у. Решив систему (10.6), мы получим решение краевой задачи.

Примечание 2
Примечание 2

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

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

Задачи на собственные значения для ОДУ

10.3. Задачи на собственные значения для ОДУ

Задачи на собственные значения — это краевые задачи для системы ОДУ, в которой правые части зависят от одного или нескольких параметров X. Значения этих параметров неизвестны, а решение краевой задачи существует только при определенных ?k, которые называются собственными значениями (eigenvalues) задачи. Решения, соответствующие этим ?k, называют собственными функциями (eigenflmctions) задачи. Правильная постановка таких задач требует формулировки количества граничных условий, равного сумме числа уравнений и числа собственных значений. Физическими примерами задач на собственные значения являются, например, уравнение колебаний струны, уравнение Шредингера в квантовой механике, уравнения волн в резонаторах и многие другие.

С вычислительной точки зрения, задачи на собственные значения очень похожи на рассмотренные выше краевые задачи. В частности, для многих из них также применим метод стрельбы (см. разд. 10.1.2). Отличие заключается в пристрелке не только по недостающим левым граничным условиям, но еще и по искомым собственным значениям. В Mathcad для решения задач на собственные значения используются те же функции sbval и bvaifit. В их первый аргумент, т. е. вектор, присваивающий начальные значения недостающим начальным условиям, следует включить и начальное приближение для собственного значения.

Рассмотрим методику решения на конкретном примере определения собственных упругих колебаний струны. Профиль колебаний струны у(х) описывается линейным дифференциальным уравнением второго порядка

где р(х) и q(x) — жесткость и плотность, которые, вообще говоря, могут меняться вдоль струны. Если струна закреплена на обоих концах, то граничные условия задаются в виде у(0)=у(1)=0. Сформулированная задача является частным случаем задачи Штурма—Лиувилля. Поскольку решается система двух ОДУ, содержащая одно собственное значение x, то по идее задача требует задания трех (2+1) условий. Однако, как легко убедиться, уравнение колебаний струны — линейное и однородное, поэтому в любом случае решение у(х) будет определено с точностью до множителя. Это означает, что производную решения можно задать произвольно, например, у' (0)=1, что и будет третьим условием. Тогда краевую задачу можно решать как задачу Коши, а пристрелку вести только по одному параметру — собственному значению.

Процедура поиска первого собственного значения представлена в листинге 10.5.

Явная схема Эйлера

11.2.1. Явная схема Эйлера

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

Построение разностной схемы

Используем для решения уравнения теплопроводности шаблон, изображенный на Рисунок 11.6. Для дискретизации второй производной по пространственной координате необходимо использовать три последовательных узла, в то время как для разностной записи первой производной по времени достаточно двух узлов. Записывая на основании данного шаблона дискретное представление для (i, k) -го узла, получим разностное уравнение:

Приведем в разностной схеме (11.8) подобные слагаемые, перенеся в правую часть значения сеточной функции с индексом k (как часто говорят, с предыдущего слоя по времени), а в левую — с индексом k+i (т. е. со следующего временного слоя). Кроме этого, введем коэффициент с, который будет характеризовать отношение шагов разностной схемы по времени и пространству

Несколько забегая вперед, заметим, что значение

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

Классификация уравнений в частных производных

11.1.1. Классификация уравнений в частных производных

Постановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений), а также необходимого количества краевых условий (число и характер задания которых определяются спецификой уравнения). По своему названию уравнения должны содержать частные производные неизвестной функции и (или нескольких функций, если уравнений несколько) по различным аргументам, например, пространственной переменной х и времени t. Соответственно, для решения задачи требуется вычислить функцию нескольких переменных, например, u(x,t) в некоторой области определения аргументов 0<x<L и 0<t<T. Граничные условия определяются как заданные временные зависимости функции и, или производных этой функции, на границах расчетной области 0 и L, а начальные — как заданная и (х, 0).

Сами уравнения в частных производных (несколько условно) можно разделить на три основных типа:

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

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

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

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

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

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

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

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

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

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

Эллиптические уравнения

11.3.2. Эллиптические уравнения

Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач — двумерного уравнения Пуассона.

Это уравнение содержит вторые производные функции u(х,у) по двум пространственным переменным;

Уравнение Пуассона описывает, например, распределение электростатического поля u(х,у) в двумерной области с плотностью заряда f (x,y), или (см. разд. 11.1.2) стационарное распределение температуры u(х,у) на плоскости, в которой имеются источники (или поглотители) тепла с интенсивностью f (х,у) .

Примечание 1
Примечание 1

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

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

Корректная постановка краевой задачи для уравнения Пуассона требует задания граничных условий. В Mathcad решение ищется на плоской квадратной области, состоящей из (M+1)х(M+1) точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон упомянутого квадрата. Самый простой вариант — это нулевые граничные условия, т. е. постоянная температура по всему периметру расчетной области. В таком случае можно использовать встроенную функцию multigrid:

multigrid(F,ncycie) — матрица решения уравнения Пуассона размера (м+1)х(м+1) на квадратной области с нулевыми граничными условиями:

 F — матрица размера (M+1)х(M+1), задающая правую часть уравнения Пуассона;  ncycle — параметр численного алгоритма (количество циклов в пределах каждой итерации).

ВНИМАНИЕ!

Сторона квадрата расчетной области должна включать точно N=2n шагов, т. е. 2n+1 узлов, где n — целое число.

Параметр численного метода ncycle в большинстве случаев достаточно взять равным 2.

Неявная схема Эйлера

11.2.2. Неявная схема Эйлера

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

Построение неявной разностной схемы

Чтобы построить неявную разностную схему для уравнения диффузии, используем шаблон, изображенный на Рисунок 11.11, т. е. для дискретизации пространственной производной будем брать значения сеточной функции с верхнего (неизвестного) слоя по времени. Таким образом, разностное уравнение для (i,k)-ro узла будет отличаться от уравнения для явной схемы (11.8) только индексами по временной координате в правой части:

Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого 1-го узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции в узах шаблона показаны на Рисунок 11.11 в виде подписей, подобно тому, как это было сделано для явной схемы (см. Рисунок 11.6).

Пример уравнение диффузии тепла

11.1.2. Пример: уравнение диффузии тепла

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

Двумерное динамическое уравнение

Рассмотрим следующее параболическое уравнение в частных производных, зависящее от трех переменных — двух пространственных х и у, а также от времени t:

Примечание 1
Примечание 1

Выражение в скобках в правой части уравнения (сумму вторых пространственных производных функции и часто, ради краткости, обозначают при помощи оператора Лапласа: ?u).

Это уравнение называется двумерным уравнением теплопроводности или, по-другому, уравнением диффузии тепла. Оно описывает динамику распределения температуры u(x,y,t) на плоской поверхности (например, на металлической пластине) в зависимости от времени (Рисунок 11.1). Физический смысл коэффициента о, который, вообще говоря, может быть функцией как координат, так и самой температуры заключается в задания скорости перетекания тепла от более нагретых областей в менее нагретые. Функция ?(x,y,t,u) описывает приток тепла извне, т. е. источники тепла, которые также могут зависеть как от пространственных координат (что задает локализацию источников), так и от времени и от температуры и.

О возможности решения многомерных уравнений

11.2.3.О возможности решения многомерных уравнений

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

Можно ли при помощи Mathcad решать двумерные или трехмерные (пространственные) уравнения? С точки зрения программирования пользователем численных алгоритмов типа метода сеток принципиальных ограничений нет. Разумеется, если сначала аккуратно выписать разностную схему соответствующего многомерного дифференциального уравнения, то вполне возможно запрограммировать ее при помощи описанных нами средств. Самым главным противодействием будет существенное увеличение времени расчетов. Простая оценка необходимого количества операций показывает, что если ввод зависимости уравнения от второй пространственной координаты многократно увеличивает число разностных уравнений, которые должны решаться при реализации каждого шага по времени. К примеру, если используется пространственная сетка из 100 узлов по каждой координате, то вместо 102 разностных уравнений на каждом шаге придется решать уже 104 уравнений, т. е. объем вычислений сразу же возрастает в 100 раз. Вообще говоря, пакет Mathcad не является экономичной средой вычислений, и бороться с их сильно возрастающим объемом следует пользователю еще на этапе разработки алгоритма. Хорошим примером такой борьбы может служить применение специфических алгоритмов, типа метода прогонки (см. разд. "Алгоритм прогонки"этой главы).

Приведем некоторые дополнительные замечания, связанные с возможностью осуществить редукцию громоздких (в смысле организации вычислений) двумерных задач к более простым. Рассмотрим ради примера двумерное уравнение теплопроводности (11.1) без источника с нулевыми граничными условиями и некоторым начальным двумерным распределением температуры по расчетной поверхности:

Произведем дискретизацию данного уравнения по временной координате, заменяя первую производную ее разностным аналогом и несколько перегруппировывая слагаемые и множители:

Как вы видите, мы используем неявную разностную схему, заранее заботясь о том, чтобы разностная задача была более устойчивой. Здесь ui(x,y) — известная с предыдущего шага по времени функция двух пространственных переменных, а и1+1(х,у) — функция, подлежащая определению при реализации каждого шага по времени.

Посмотрим на полученную задачу с несколько другой стороны — а именно как на дифференциальное уравнение относительно неизвестной функции двух переменных ui+1(x,y). Подчеркнем, что такое уравнение получается для каждого шага по времени, т. е. для реализации всей разностной схемы требуется решить большое число таких уравнений.

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

Суммируя сказанное, можно констатировать, что если вы имеете запрограммированный алгоритм решения выписанного эллиптического уравнения, чуть более сложного, чем уравнение Пуассона, то его с легкостью можно использовать в качестве подпрограммы реализации разностной схемы двумерного уравнения теплопроводности. Забегая вперед, приходится отметить, что, к сожалению, встроенные функции Mathcad для решения уравнения Пуассона в данном случае не годятся в качестве такой подпрограммы, поскольку предполагают независимость источника от самой неизвестной функции, и могут справляться лишь с правой частью, зависящей только от пространственных координат (см. разд. 11.3.2).

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

Численное решение обратного уравнения

Рисунок 11.5. Численное решение обратного уравнения теплопроводности дает совершенно нефизичную картину динамики температуры (см. листинг 11.2 ниже с параметром D=-1)

Численное решение уравнения теплопроводности

Рисунок 11.10. Численное решение уравнения теплопроводности при помощи явной схемы Эйлера (см. листинг 11.1 ниже с временным шагом t=0.0015)

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

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

Примечание 4
Примечание 4

Как несложно убедиться, для t=0.0005 коэффициент Куранта C=0.4, для t=0.0010 он все еще меньше единицы: C=0.8, а для t=0. 0015 решение уже больше единицы: C=1.2, в связи с чем схема становится неустойчивой (см. Рисунок 11.10).

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

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

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

Дифференциальные уравнения в частных производных требуют нахождения функции не одной, как для ОДУ, а нескольких переменных, например, f (х,у) или f(x,t). Постановка задач (см. разд. 11.Т) включает в себя само уравнение (или систему уравнений), содержащее производные неизвестной функции по различным переменным (частные производные), а также определенное количество краевых условий на границах расчетной области.

Несмотря на то, что Mathcad обладает довольно ограниченными возможностями по отношению к уравнениям в частных производных, в нем имеется несколько встроенных функций (см. разд. 11.3). Решать уравнения в частных производных можно и путем непосредственного программирования пользовательских алгоритмов (см. разд. 11.2). Автор совершенно сознательно сначала рассматривает численные методы для решения уравнений в частных производных, а уже затем описывает предназначенные для этого встроенные функции, чтобы читатель ясно осознавал, каким образом Mathcad производит расчеты. "Слепое" использование встроенных функций для решения уравнений в частных производных не всегда бывает успешным, и ответственность за верный выбор их параметров часто ложится на пользователя, которому необходимо четко представлять основные принципы функционирования численных алгоритмов, примененных во встроенных функциях.

Физическая модель одномерного уравнения теплопроводности

Рисунок 11.3. Физическая модель одномерного уравнения теплопроводности

Физическая модель описываемая двумерным уравнением теплопроводности

Рисунок 11.1. Физическая модель, описываемая двумерным уравнением теплопроводности

Для того чтобы правильно поставить краевую задачу для двумерного уравнения теплопроводности, следует определить следующие дополнительные условия:

граничные условия, т. е. динамику функции u(x,y,t) и (или) ее производных на границах расчетной области;  начальное условие, т. е. функцию u(х, у, t).

Примечание 2
Примечание 2

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

Стационарное двумерное уравнение

Частный случай уравнения теплопроводности определяет стационарную, т. е. не зависящую от времени, задачу. Стационарное уравнение описывает физическую картину распределения температуры по пластине, не изменяющуюся с течением времени. Такая картина может возникнуть при условии, что стационарный источник тепла действует довольно продолжительное время, и переходные процессы, вызванные его включением, прекратились. Пример численного решения такого уравнения показан на Рисунок 11.2 в виде поверхности u(х,у).

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

Глава 11. Дифференциальные уравнения в частных производных

Дифференциальные уравнения в частных производных 11.1. О постановке задач 11.1.1. Классификация уравнений в частных производных 11.1.2. Пример: уравнение диффузии тепла 11.2. Разностные схемы 11.2.1. Явная схема Эйлера 11.2.2. Неявная схема Эйлера 11.2.3. О возможности решения многомерных уравнений 11.3. Встроенные функции для решения уравнений в частных производных 11.3.1. Параболические и гиперболические уравнения 11.3.2. Эллиптические уравнения

Содержание

График линий уровня решения уравнения Пуассона (продолжение листинга 11 7)

Рисунок 11.18. График линий уровня решения уравнения Пуассона (продолжение листинга 11.7)

Уравнение Пуассона с произвольными граничными условиями

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

relax(a,b,c,d,e,F,v,rjac) — матрица решения дифференциального уравнения в частных производных на квадратной области, полученного с помощью алгоритма релаксации для метода сеток:

 a,b,c,d,e — квадратные матрицы коэффициентов разностной схемы, аппроксимирующей дифференциальное уравнение;  F — квадратная матрица, задающая правую часть дифференциального уравнения;  v — квадратная матрица граничных условий и начального приближения к решению; rjac — параметр численного алгоритма (спектральный радиус итераций Якоби);

Параметр численного алгоритма характеризует скорость сходимости итераций. Он должен быть числом от о до 1. В матрице граничных условий v необходимо задать только граничные элементы, исходя из значения краевых условий по периметру расчетной области. Прочие (внутренние) элементы этой матрицы служат для задания начального приближения к решению.

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

ВНИМАНИЕ!

Все матрицы, задающие как коэффициенты разностной схемы а, b, с, d, e, граничные условия v, так и само решение F, должны иметь одинаковый размер (M+1)х(M+1), соответствующий размеру расчетной области. При этом целое число м обязательно должно быть степенью двойки: м=2п.

Решение уравнения Пуассона с тремя источниками разной интенсивности при помощи функции relax приведено в листинге 11.8.

График поверхности решения уравнения Пуассона (продолжение листинга 11 7)

Рисунок 11.17. График поверхности решения уравнения Пуассона (продолжение листинга 11.7)

Явная схема для линейного уравнения теплопроводности

Листинг 11.1. Явная схема для линейного уравнения теплопроводности

Начальное распределение температуры вдоль расчетной области и решение для двух моментов времени показано на Рисунок 11.7 сплошной, пунктирной и штриховой линиями соответственно. Физически такое поведение вполне естественно — с течением времени тепло из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает и размывается.

Неявная схема для линейного уравнения теплопроводности

Листинг 11.2. Неявная схема для линейного уравнения теплопроводности

Алгоритм прогонки

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

Примечание 2
Примечание 2

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

Основным вычислительным ядром программы, реализующей на Mathcad неявную разностную схему, было решение (на каждом временном слое) системы линейных алгебраических уравнений, задаваемых матрицей А. Заметим, что эта матрица, как говорят, имеет диагональное преобладание, а точнее, является трехдиагоналъной (Рисунок 11.13). Все ее элементы, кроме элементов на главной диагонали и двух соседних диагоналях, равны нулю. С точки зрения оптимизации быстродействия алгоритма применение встроенной функции isolve является весьма расточительным, поскольку основной объем арифметических операций, выполняемых компьютером (а он составляет, как нетрудно убедиться, величину порядка м2), сводится к непроизводительному перемножению нулей.

Алгоритм прогонки (продолжение листинга 11 2)

Листинг 11.3. Алгоритм прогонки (продолжение листинга 11.2)

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

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

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

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

Примечание 1
Примечание 1

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

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

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

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

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

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

Решение уравнения Пуассона с нулевыми граничными условиями

Листинг 11.7.Решение уравнения Пуассона с нулевыми граничными условиями

В первой строке листинга задается значение м=32, в двух следующих строках создается матрица правой части уравнения Пуассона, состоящая из всех нулевых элементов, за исключением одного, задающего расположение источника. В последней строке матрице с присваивается результат действия функции multigrid. Обратите внимание, первый ее аргумент сопровождается знаком "минус", что соответствует записи правой части уравнения Пуассона (11.11). Графики решения показаны на Рисунок 11.17 и 11.18 в виде трехмерной поверхности и линий уровня соответственно.

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

Листинг 11.7 содержит пример использования функции multigrid для расчета краевой задачи на области 33х33 точки и точечным источником тепла в месте, задаваемом координатами (15,20) внутри этой области.

Решение уравнения Пуассона с помощью функции relax

Листинг 11.8. Решение уравнения Пуассона с помощью функции relax

Первые три строки имеют тот же смысл, что и в предыдущем листинге. Только вместо одного источника тепла взято их другое распределение — один сильный источник, один более слабый и один сток тепла. В следующих шести строках задаются коэффициенты разностной схемы. Отложим их обсуждение до последнего раздела этой главы, ограничившись утверждением, что для решения уравнения Пуассона коэффициенты должны быть взяты именно такими, как показано в листинге 11.8. В предпоследней строке задана матрица нулевых граничных условий и нулевых начальных приближений, а в последней матрице с присваивается результат действия функции relax. График полученного решения в виде линий уровня показан на Рисунок 11.19.

Решение уравнения теплопроводности при помощи функции relax

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

Результат действия программы листинга 11.9 показан на Рисунок 11.21 в виде трехмерной поверхности. Если сравнить Рисунок 11.21 с Рисунок 11.4, полученным при расчетах по запрограммированной разностной схеме, то в графиках Рисунок 11.4 нетрудно узнать сечения этой поверхности плоскостями t=const. Еще раз подчеркнем, что использовать встроенную функцию можно только для тех уравнений, которые допускают построение разностной схемы типа "крест" (см. Рисунок 11.17) или составного фрагмента этой схемы.

Матрица системы линейных разностных уравнений для неявной схемы (листинг 11 2 для M=10)

Рисунок 11.13. Матрица системы линейных разностных уравнений для неявной схемы (листинг 11.2 для M=10)

Для отыскания решения линейных систем алгебраических уравнений имеется чрезвычайно эффективный алгоритм, называемый прогонкой, который позволяет уменьшить число арифметических операций на целый порядок, т. е. до значения порядка м. Это означает, что при использовании пространственных сеток с юоо узлами выигрыш во времени вычислений составит величину порядка ю3! Реализация данного алгоритма приведена в листинге 11.3, который является продолжением листинга 11.2, используя определенные в нем коэффициенты матрицы А, а также начальное условие.

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

О постановке задач

11.1.О постановке задач

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

Разностные схемы

11.2. Разностные схемы

Рассмотрим одномерное уравнение теплопроводности (11.3) и на его примере разберем наиболее часто использующийся для численного решения уравнений в частных производных метод сеток. Выпишем еще раз само уравнение

а также и начальное

u(x,0) = u0(x), (11.6)

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

u(0, t) = f0(t) , u(L, t) = f1(t), (11.7)

которые необходимы для правильной с математической точки зрения постановки задачи.

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

Таким образом, сначала следует покрыть расчетную область (x,t) сеткой и использовать затем узлы этой сетки для разностной аппроксимации уравнения. В результате, вместо поиска непрерывных зависимостей u(x,t) достаточно будет отыскать значения функции в узлах сетки (а ее поведение в промежутках между узлами может быть получено при помощи построения какой-либо интерполяции). По этой причине дискретное представление функции и часто называют сеточной функцией.

Поскольку уравнения в частных производных по определению зависят от производных неизвестных функций по нескольким переменным, то способов дискретизации этих уравнений может быть, как правило, несколько. Конфигурацию узлов, используемую для разностной записи уравнений в частных производных на сетке, называй шаблоном, а полученную систему разностных уравнений — разностной схемой. О принципах построения разностных схем и, в частности, о классах явных и неявных схем мы уже подробно говорили на примере краевых задач для ОДУ (см. разд. 10.4.1), поэтому, излишне не повторяясь, перейдем к рассмотрению типичных особенностей уравнений в частных производных, которые возникают при разработке и реализации разностных схем.

Решение линейного уравнения теплопроводности (продолжение листинга 11 1)

Рисунок 11.7. Решение линейного уравнения теплопроводности (продолжение листинга 11.1)

Примечание 2
Примечание 2

Несколько забегая вперед, заметим, что показанное на Рисунок 11.7 решение и соответствует коэффициенту Куранта С=0.4. Попробуйте осуществить расчет с увеличенным временным шагом, чтобы коэффициент с был больше 1, и посмотрите, что из этого получится (такой расчет и его объяснение приведены ниже в разд. "Устойчивость").

Нелинейное уравнение

Намного более интересные решения можно получить для нелинейного уравнения теплопроводности, например, с нелинейным источником тепла ф(u)=103(u-u3). Заметим, что в листинге 11.1 мы предусмотрительно определили коэффициент диффузии и источник тепла в виде пользовательских функций, зависящих от аргумента и, т. е. от температуры. Если бы мы собирались моделировать явную зависимость их от координат, то следовало бы ввести в пользовательскую функцию в качестве аргумента переменную х, как это сделано для источника тепла ф. Поэтому нет ничего проще замены определения этих функций с констант D(U)=1 и ф(х,u)=о на новые функции, которые станут описывать другие модели диффузии тепла. Начнем с того, что поменяем четвертую строку листинга 11.1 на ф(х,и)=ю3-(и-и3), не изменяя пока постоянного значения коэффициента диффузии.

Примечание 3
Примечание 3

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

Если осуществить расчеты с упомянутым источником (имеющим кубическую нелинейность), то получится очень интересное решение уравнения теплопроводности, имеющее профиль тепловых фронтов. С течением времени граница раздела высокой и низкой температуры распространяется в обе стороны от зоны первичного нагрева, оставаясь весьма четко выделенной (Рисунок 11.8).

Еще более неожиданные решения возможны при нелинейности также и коэффициента диффузии. Например, если взять квадратичный коэффициент диффузии D(x,u)=u2 (что с учетом его умножения на неизвестные функции создаст кубическую нелинейность уравнения), а также ф(х,u)=103-u3-5, то вы сможете наблюдать совсем иной режим горения среды. В отличие от рассмотренного эффекта распространения тепловых фронтов, горение оказывается локализованным в области первичного нагрева среды, причем температура в центре нагрева со временем возрастает до бесконечной величины (Рисунок 11.9). Такое решение описывает так называемый режим горения "с обострением".

Решение линейного

Рисунок 11.12. Решение линейного уравнения теплопроводности при помощи неявной схемы на первом слое по времени (листинг 11.2)

Решение одномерного уравнения теплопроводности (см листинг 11 1 ниже)

Рисунок 11.4. Решение одномерного уравнения теплопроводности (см. листинг 11.1 ниже)

Линейное и нелинейное уравнения

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

Линейная задача — если коэффициент диффузии о не зависит от температуры и и, кроме того, если источник тепла ф либо также не зависит от и, либо зависит от и линейно. В этом случае неизвестная функция u (x, t) и все ее производные входят в уравнение только в первой степени (линейно).  Нелинейная задача — если уравнение имеет нелинейную зависимость от u(x,t), т. е. или коэффициент диффузии зависит от и, и (или) источник тепла нелинейно зависит от и.

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

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

Примечание 3
Примечание 3

Различные линейные и нелинейные варианты рассматриваемого уравнения теплопроводности описывают различные модели физических сред, которые характеризуются определенными зависимостями D(u) и ф(и). В частности, для металлов в большинстве случаев можно считать, что D=const, в то время как для плазмы имеется специфическая зависимость коэффициента диффузии от температуры.

Обратное уравнение теплопроводности

Замечательными свойствами обладает так называемое обратное уравнение диффузии тепла, которое получается путем замены в исходном (прямом) уравнении переменной t на -t. Согласно постановке задачи, обратное уравнение теплопроводности описывает реконструкцию динамики профиля температуры остывающего стержня, если известно начальное условие в виде профиля температуры в некоторый момент времени после начала остывания. Таким образом, требуется определить, как происходило остывание стержня. Мы ограничимся самым простым линейным уравнением с D=const без источников тепла:

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

Если попробовать осуществить расчет обратного уравнения диффузии тепла по тем же самым алгоритмам, что и для обычных уравнений (для этого достаточно в листинге 11.1 или 11.2 заменить значение коэффициента диффузии на отрицательное число, например, D=-1), то мы получим заведомо нефизичное решение. Оно показано на Рисунок 11.5 в виде профилей распределения температуры для нескольких последовательных моментов времени. Как видно, решение выражается в появлении все более быстрых пространственных осцилляции профиля температуры для каждого нового момента времени. Очень существенно, что такое решение является не проявлением неустойчивости численного алгоритма (как для ситуации, рассмотренной в разд. "Устойчивость"этой главы), а определяется спецификой самой задачи.

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

Решение стационарного двумерного уравнения теплопроводности (см листинг 11 7 ниже)

Рисунок 11.2. Решение стационарного двумерного уравнения теплопроводности (см. листинг 11.7 ниже)

Как несложно сообразить, если искомая функция не зависит от времени, то частная производная по времени в левой части уравнения равна нулю, и само уравнение можно переписать (переобозначив ради упрощения ?<-?/D) следующим образом:

Полученное уравнение, согласно классификации предыдущего раздела, является эллиптическим. Его называют уравнением Пуассона, а для его решения в Matcad предусмотрены две встроенные функции. Если, к тому же, источники равны нулю, то уравнение (11.2), принимающее вид ?u=0, называют уравнением Лапласа.

Одномерное динамическое уравнение

Предположим, что мы рассматриваем задачу распределения тепла не по плоской поверхности, а по удлиненному телу типа металлического стержня (Рисунок 11.3). В этом случае зависимость от координаты у в общем уравнении теплопроводности пропадает, и получается одномерное уравнение:

Одномерное уравнение намного проще двумерного, поскольку объем вычислений для реализации алгоритма его численного решения не так велик. Типичное решение одномерного уравнения диффузии тепла с коэффициентом диффузии о=2, нулевым источником ф=о и начальным распределением температуры в форме нагретой центральной области стержня показано (в виде графика поверхности) на Рисунок 11.4.

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

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

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

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

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

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

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

Решение уравнения Пуассона с помощью функции relax (продолжение листинга 11 8)

Рисунок 11.19. Решение уравнения Пуассона с помощью функции relax (продолжение листинга 11.8)

Разностная схема для решения уравнения Пуассона

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

Начнем с пояснения выбора этих коэффициентов (см. листинг 11.8) для уравнения Пуассона. Согласно основным идеям метода сеток (см. разд. 11.2), для дискретизации обеих пространственных производных в уравнении (11.12) следует использовать по три соседних узла вдоль каждой из координат. Поэтому уравнение Пуассона (11.12) может быть записано в разностной форме при помощи шаблона типа "крест" (Рисунок 11.20). В этом случае, после приведения подобных слагаемых в разностных уравнениях коэффициенты разностной схемы будут такими, как показано возле узлов шаблона на этом рисунке (аналогичные коэффициенты для явной и неявных схем решения уравнения теплопроводности см. на Рисунок 11.6 и 11.11 соответственно).

Теперь, если вы сравните полученные числа с константами, которые присвоены элементам матриц-аргументов функции relax (см. листинг 11.8), то увидите, что они как раз и описывают вычисленные нами только что коэффициенты разностной схемы "крест". Таким образом, нетрудно сообразить, что с помощью встроенной функции relax можно решать и другие линейные дифференциальные уравнения в частных производных, которые можно аппроксимировать схемой типа "крест" или схемой, являющейся ее составной частью. Конечно, для того чтобы использовать эту встроенную функцию для другого уравнения, необходимо будет составить соответствующую разностную схему.

Решение уравнения теплопроводности

Рисунок 11.9. Решение уравнения теплопроводности с нелинейным источником и коэффициентом диффузии (режим локализации горения)

Читателю предлагается поэкспериментировать с этим и другими нелинейными вариантами уравнения теплопроводности. Существенно, что такие интересные результаты удается получить лишь численно, а в Mathcad только с применением элементов программирования.

Устойчивость

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

Слегка изменим соотношение шагов по времени и пространственной координате, произведя расчеты сначала с t=0.0005 (эти результаты показаны на Рисунок 11.7 выше), а также с t=0.0010 и t=0.0015. Результат применения разностной схемы Эйлера для t=0.0010ю примерно тот же, что и для меньшего значения т, приведенного на Рисунок 11.7. А вот следующее (казалось бы, незначительное) увеличение шага по времени приводит к катастрофе (Рисунок 11.10).

Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)

Рисунок 11.8. Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)

Решение уравнения теплопроводности с помощью функции relax (продолжение листинга 11 9)

Рисунок 11.21. Решение уравнения теплопроводности с помощью функции relax (продолжение листинга 11.9)

Решение волнового уравнения (продолжение листинга 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)

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

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

Шаблон аппроксимации явной схемы для уравнения теплопроводности

Рисунок 11.6. Шаблон аппроксимации явной схемы для уравнения теплопроводности

Множители для каждого из значений сеточной функции в узлах шаблона, соответствующие разностному уравнению (11.9), приведены рядом с каждой точкой шаблона на Рисунок 11.6. Фактически геометрия шаблона и эти множители задают построенную нами разностную схему.

Несложно убедиться в том, что для получения замкнутой системы разностных алгебраических уравнений систему (11.9) необходимо дополнить дискретным представлением начального и граничных условий (11.6) и (11.7). Тогда число неизвестных будет в точности равно числу уравнений, и процесс формирования разностной схемы будет окончательно завершен.

Примечание 1
Примечание 1

Важно подчеркнуть, что возможная нелинейность полученной системы алгебраических уравнений определяется зависимостями от температуры функций D(u) и ф (u), т. е. как коэффициент диффузии, так и источник тепла могут быть функциями сеточной функции ui,k.

Если присмотреться к разностным уравнениям (11.9) повнимательнее, то можно сразу предложить несложный алгоритм реализации этой разностной схемы. Действительно, каждое неизвестное значение сеточной функции со следующего временного слоя, т. е. левая часть соотношения (11.9) явно выражается через три ее значения с предыдущего слоя (правая часть), которые уже известны. Таким образом, в случае уравнения теплопроводности нам очень повезло. Для расчета первого слоя по времени следует попросту подставить в (11.9) начальное условие (известные значения и с нулевого слоя в узлах сетки), для расчета второго слоя достаточно использовать вычисленный таким образом набор и с первого слоя и т. д. Из-за того что разностная схема сводится к такой явной подстановке, ее и называют явной, а благодаря пересчету значений с текущего слоя через ранее вычисленные слои — схемой бегущего счета.

Линейное уравнение

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

Решение системы разностных уравнений (11.9) для модели без источников тепла, т.е. ф(x,T,t)=0 и постоянного коэффициента диффузии D=const приведено в листинге 11.1. В его первых трех строках заданы шаги по временной и пространственной переменным t и А, а также коэффициент диффузии о, равный единице. В следующих двух строках заданы начальные (нагретый центр области) и граничные (постоянная температура на краях) условия соответственно. Затем приводится возможное программное решение разностной схемы, причем функция пользователя v(t) задает вектор распределения искомой температуры в каждый момент времени (иными словами, на каждом слое), задаваемый целым числом t.

Шаблон аппроксимации уравнения Пуассона "крест"

Рисунок 11.20. Шаблон аппроксимации уравнения Пуассона "крест"

Решение уравнения диффузии тепла при помощи функции relax

Приведем пример применения встроенной функции relax для решения другого уравнения в частных производных (т. е. не уравнения Пуассона, для которого она изначально предназначена). Вычислим при помощи этой функции решение уже хорошо нам знакомого однородного линейного уравнения теплопроводности (см. разд. 11.1.2). Будем использовать явную разностную схему, шаблон которой изображен на Рисунок 11.6. Для того чтобы "приспособить" для явной схемы функцию relax, требуется только задать ее аргументы в соответствии с коэффициентами, показанными на шаблоне (см. Рисунок 11.6). Программа, реализующая таким способом явную схему, представлена на листинге 11.9. Число Куранта в этом листинге обозначено переменной с, как и положено явной разностной схеме, она выдает устойчивое решение только для C<1.

Шаблон неявной схемы для уравнения теплопроводности

Рисунок 11.11. Шаблон неявной схемы для уравнения теплопроводности

Очень важно, что если само уравнение теплопроводности линейно, то с в левой части разностного уравнения является константой, а ф в его правой части может зависеть только от первой степени и. Поэтому система уравнений (11.10) для всех пространственных узлов 1=1.. .м-l является линейной системой, что существенно упрощает ее решение (поскольку известно, что для линейных систем с ненулевым определителем решение существует и является единственным). Напомним, что для получения замкнутой системы линейных уравнений необходимо дополнить данный набор разностных уравнений граничными условиями, т. е. известными значениями сеточной функции для i=0 и i=M.

Примечание 1
Примечание 1

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

Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования Mathcad и встроенной функции решения системы линейных уравнений isolve. Один из возможных способов решения предложен в листинге 11.2. Большая часть этого листинга является вводом параметров задачи (шагов, начальных и граничных условий), и только в последней его строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения системы линейных уравнений isolve). В нескольких предыдущих строках листинга (после расчета коэффициента Куранта сои) формируется матрица системы уравнений А, которая записывается в подходящем для Mathcad виде, как это сделано в листинге 11.2. Как несложно убедиться, столбец правых частей разностных уравнений выражается вычисленными значениями сеточной функции с предыдущего слоя.

Результаты расчетов по неявной схеме показаны на Рисунок 11.12 и, как видно, они дают примерно те же результаты, что и в случае применения явной схемы (см. Рисунок 11.7). Обратите внимание, что решение устойчиво при любых значениях коэффициента Куранта (в том числе и больших 1), поскольку, как следует из соответствующих положений теории численных методов, неявная схема является безусловно-устойчивой.

Встроенные функции для решения уравнений в частных производных

11.3. Встроенные функции для решения уравнений в частных производных

Как видно из предыдущего раздела, с уравнениями в частных производных вполне можно справиться и не прибегая к специфическим средствам Mathcad. Между тем в версиях Mathcad 11 (и выше) имеется несколько встроенных функций, при помощи которых можно автоматизировать процесс решения дифференциальных уравнений в частных производных. Рассмотрим в данном разделе основные аспекты их применения, отмечая не только инструкции по их применению, описанные разработчиками Mathcad, но также и некоторые "скрытые" возможности этих функций.

---

---