Приближенные методы решения дифференциальных уравнений

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

На заметку
К сожалению, далеко не в каждом уравнении такой малый параметр можно выделить.

В качестве примера рассмотрим следующую задачу.

Задача 5.7

Найти приближенное решение в виде многочлена второго порядка по малому параметру для задачи Коши.

Задаем исходное уравнение.

Поскольку это уравнение первого порядка, начальное условие только одно.

Строго говоря, данное уравнение вычислительным ядром Maple решается точно. Ниже приведена соответствующая команда, однако, без указания результата. Причина проста -— результат этот весьма нетривиален. Желающие могут убедиться в этом самостоятельно. Последнее, кстати, является свидетельством того, что точный результат искать не всегда полезно; зачастую достаточно ограничиться приближенным решением — оно может оказаться вполне приемлемым по точности и в то же время простым и наглядным.

Решение ищем в виде ряда по малому параметру — в данном случае это е. Поэтому у(х) представляем в следующем виде.

Внимание!
Представленная выше команда, с помощью которой функции у(х) присваивается значение, на самом деле функцию не определяет. Если в командной строке ввести команду у(х), в области вывода появится уО(х) + еу\(х) + егу2(х). Но если команду заменить, скажем, на y(t), результат будет y(t). Другими словами, присваивание выполняется "на уровне названий", а у(х) в данном случае является не результатом действия оператора у() на переменную х, а названием переменной.

Уравнение, таким образом, будет иметь следующий вид.

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

Теперь можно собрать слагаемые при соответствующих степенях малого параметра (epsilon).

Неизвестные функции yO(x), yl(x) и у2(х) выбираются так, чтобы коэффициенты при степенях малого параметра (не всех, а только до второго включительно, до которого раскладывается в ряд искомая функция) равнялись нулю.

Отсюда получаем первое уравнение (для определения уО(х)), приравняв правую часть уравнения Eq (это коэффициент при нулевой степени параметра epsilon) к нулю.

Приравнивая к нулю коэффициент при первой степени малого параметра, получаем еще одно уравнение.

В этом уравнении, помимо yl(x), присутствует и уО(х). Поэтому, прежде чем решать уравнение Eq_l, необходимо сначала решить уравнение Eq 0, найти тем самым уО(х), после чего можно будет решать Eq_l относительно yl(x).

На заметку
Коэффициент при первой степени параметра epsilon определяется процедурой coef f (). Первым ее аргументом указана левая часть уравнения Eq (lhs(Eq)) — именно в левой части содержатся слагаемые с epsilon. Второй параметр (epsilon) является указанием на то, что нужно выделить коэффициент при первой степени epsilon. Приравняв этот коэффициент к нулю, получаем нужное уравнение.

Практически точно так же получаем последнее дифференциальное уравнение — приравниваем к нулю коэффициент при второй степени epsilon. Отличие от предыдущего случая состоит в том, что вторым параметром в процедуре coeff() указано epsilonA2.

Решать уравнения Eq_0, Eq_l и Eq 2 следует строго в той последовательности, в какой они вводились. Причина очевидна — каждое последующее уравнение содержит, помимо неизвестной функции, еще и функции, которые определяются предыдущими уравнениями.

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

Итак, получаем следующее.

Ниже, согласно полученному решению, выполняется присваивание.

Внимание!
В этом месте функция уО(х) "технически" становится переменной!

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

Однако если в приведенном выше вызове изменить аргумент (yO(t)), нужного значения (1/t) не получим.

Решаем теперь уравнение Eg 1 и находим yl(x).

В данном случае, перед присваиванием, раскладываем выражение для yl(x) на сумму дробей (процедура expand(%)).

Наконец, находим у2(х). Последовательность действий такая же.

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

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

Задача 5.8

Найти приближенное решение в виде многочлена второго порядка по малому параметру для задачи Коши 1 + х + еу, у(0) = sin(s).

Как и прежде, в первую очередь задаем само уравнение.

Начальные условия в этом случае определяются через малый параметр.

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

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

Слагаемые добавлять будем с помощью оператора цикла.

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

Исходное уравнение при этом будет иметь такой вид.

Преобразуем это уравнение, перенеся все слагаемые в левую часть.

Далее левую часть уравнения Eq раскладываем в степенной ряд по параметру epsilon в окрестности нуля (команда series(lhs(Eq),epsilon=0,N+l)). Поскольку нас интересуют степени epsilon вплоть до N, в процедуре series () третьим параметром указано N+1 (на единицу больше) — этот параметр определяет порядок остатка, а сам ряд имеет порядок на единицу меньше. Разложение преобразуем в полиномиальный вид (процедура convert()). Чтобы уравнение осталось уравнением, полученное выражение приравниваем к нулю.

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

Названия этих переменных являются объединением имени "Eq" и индекса i, который изменяется в диапазоне от 0 до N.

Коэффициент при i-й степени epsilon определяется процедурой coeff(). Степень epsilon указывается третьим параметром процедуры.

На следующем этапе выполняем разложение в ряд по малому параметру начальных условий. В данном случае начальное условие одно, и записано оно в переменную inCon. Эта переменная — равенство. Поэтому переносим все слагаемые в одну сторону — это делается командой (lhs(InCon)-rhs(InCon)) — и полученное выражение раскладывается в ряд по epsilon в окрестности нуля.

Теперь выражение трансформируется в полином.

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

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

Теперь выполним замену.

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

Решаем первое уравнение (приняв во внимание начальные условия).

Описываем процедуру spDiffSol(), которая будет иметь пять параметров: решаемое дифференциальное уравнение Eq, начальные условия InCon, функция fnc (с указанием аргумента!), относительно которой решается уравнение, малый параметр epsilon и целое неотрицательное число п, определяющее порядок разложения по малому параметру при нахождении решения.

После каждого из перечисленных параметров стоит двойное двоеточие (::), после которого указывается тип параметра. Так, equation (этот тип имеют первый и второй параметры) определяет уравнение. Третий параметр имеет тип функция (function). Малый параметр имеет символьный тип (symbol) и, наконец, последний, пятый, параметр — целое неотрицательное число (тип nonnegint). Если в дальнейшем при вызове процедуры хотя бы один из ее параметров будет иметь иной тип, появится сообщение об ошибке. Если напрямую тип не указывать, будет предпринята попытка выполнить процедуру и, скорее всего, ошибка возникнет непосредственно в процессе выполнения — такие ошибки отслеживать сложнее.

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

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

Первой определяется переменная х, которой в качестве значения присваивается операнд указанной аргументом функции fnc (op(fnc)), т.е. ее аргумент. Переменной у присваивается нулевой операнд функции fnc (команда op(0,fnc)); для функциональных зависимостей таким операндом является имя функции. Эти переменные, следовательно, определены так, чтобы тождественно выполнялось равенство y(x)=fnc. Кроме того, для удобства малый параметр запишем в переменную.

На следующем этапе в теле процедуры определяется другая процедура eps() — локальная, т.е. вне процедуры spDiffSol данная процедура будет недоступна! — при действии которой на аргумент t (т.е. процедурой определяется оператор) этот аргумент умножается на первый параметр процедуры г, возведенный в степень, определяемую вторым параметром процедуры i. После определения процедуры eps() формируется сумма F. Сумма эта операторная: слагаемыми суммы являются операторы. Каждый такой оператор есть композиция (оператор композиции) двух операторов: a[i] и eps(z,i), где i — индексная переменная, изменяющаяся в диапазоне от 0 до п. Действие оператора F на аргумент будет, согласно его описанию, следующим: на аргумент сначала действует оператор a[i], а результат этого действия умножается на z(в этом состоит действие оператора eps(z,i)). Общий результат формируется как сумма от действия каждого из слагаемых.

На заметку
В выражении для F использован оператор композиции §. Если G и f — функциональные операторы (т.е. G(x) и f (x) являются функциями), то оператор L, определенный как L:=G§f, действует следующим образом: L(x)=G(f (х)).

Описанные выше функциональные операторы a[i] будут формировать коэффициенты разложения искомой функции в ряд по малому параметру, и именно их нужно определить в процессе решения.

Далее вводится переменная eq, значение которой формируется следующим образом. Сначала в исходном решаемом уравнении Eq все слагаемые переносятся в левую часть (команда lhs(Eq)-rhs(Eq)), а в полученном таким образом выражении функция (fnc) заменяется с помощью процедуры subs() разложением в ряд. Разложение в ряд получается в результате действия оператора F на аргумент функции fnc (переменная х, а результат равен F(x)). После этого выражение eq раскладывается в ряд (процедура series()) по малому параметру в окрестности нуля (это нужно сделать на тот случай, если в выражении присутствуют отличные от степенных функции малого параметра). Остаток ряда должен иметь порядок п+1. Далее выполняется преобразование в полиномиальный вид (процедура convert ()). Из этого выражения будут сформированы уравнения для определения коэффициентов a[i].

Правая часть уравнения inCon переносится влево, и это выражение присваивается переменной S. После этого в выражении S оператор у() заменяется оператором F(). Для того чтобы оператор F возымел действие, вызывается процедура value () и результат ее выполнения присваивается переменной S. Таким образом, переменная S является тем выражением, которое определяет начальное условие и в котором действие оператора у() заменено действием операторного ряда F().

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

Выражение с начальными условиями далее раскладывается в ряд по малому параметру и трансформируется в полином.

После выполнения описанных выше преобразований с учетом начального условия, приступаем непосредственно к решению задачи. Для этого задействуем оператор цикла. В теле этого оператора первой командой формируется уравнение для определения очередного коэффициента a[i](x). С этой целью приравниваем к нулю полиномиальные коэффициенты разложения выражения eq по малому параметру. Полученные таким образом уравнения присваи-I ваем в качестве значения переменным, названия которых формируются путем | объединения базового имени eqn и индексной переменной i, принимающей значения в диапазоне от 0 до п. Точно так же поступаем и с начальными условиями, только в этом случае к нулю приравниваются коэффициенты в выражении S, а соответствующие уравнения записываются в переменные КО, IC1, .., 1Сп. Названия переменных создаются, как и ранее, посредством оператора конкатенации (| |). Наконец, уравнение (eqn| |i) и соответствующее ему граничное условие (1С 11 i) объединяются во множество, которое присваивается в качестве значения переменной Eq_and_ic| |i.

Далее следует процедура dsolve(), с помощью которой уравнение eqn| |i с граничным условием IC||i (они являются элементами множества Eq_and_IC11i) решается относительно функции a[i)(x). Поскольку результат решения представляется в виде равенства, после нахождения решения следует, согласно равенству-решению, определить функцию a[i](x) (если точнее, то оператор a[ij). Результат действия оператора a[i] на аргумент х определяется правой частью равенства, полученного в результате решения уравнения. Ссылка на этот результат выполнена посредством переменной среды %, а правая часть равенства, соответственно, возвращается командой rhs(%). Чтобы из этого выражения выделить оператор, используем процедуру unapply(). На этом команды оператора цикла заканчиваются.

По завершении описанного выше оператора цикла все уравнения решены, и в переменную Res записывается результат вычисления F(x).

Внимание!
Если не использовать процедуру eval(), а просто указать Res:=F(x), то результат выполнения процедуры будет представляться через коэффициенты a[i], без непосредственного вычисления результата их действия на аргумент х.

Последней командой процедуры является Res, в результате выполнения которой на экран выводится приближенное решение уравнения.

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

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

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

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

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

Задача 5.9

Найти решение в виде ряда для уравнения ху"(х)+у'(х)+ху(х) = чальными условиями у(0) = 1(0) = 0. с начальными условиями y(0)=1, y(0)=0.

Это уравнение имеет решение — функцию Бесселя j(x). Ниже попытаемся найти приближенное решение и сравнить его с точным. Как обычно, задаем уравнение.

Начальные условия для этого уравнения имеют следующий вид.

Если решать данную задачу с помощью процедуры dsolvef), получим следующий результат.

Чтобы получить приближенное решение в виде ряда, используем все ту же процедуру dsolve() практически с теми же параметрами; в конце добавлена опция series. Решение будет представлено рядом.

Количество слагаемых в ряде определяется значением переменной среды Order; значение этой переменной задает порядок остатка. Например, если нужно, чтобы последнее слагаемое было степени 9 по х, переменной Order присваиваем значение на единицу больше (т.е. 10).

В этом случае приближенное решение будет следующим.

То, что последнее слагаемое имеет степень 8, объясняется просто: в разложении присутствуют только четные степени х, т.е. коэффициенты при нечетных степенях — нули.

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

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

Ниже показаны кривые для точного и приближенного решений.

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

Содержание раздела

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