Задачи Коши для обыкновенных дифференциальных уравнений
Дифференциальные уравнения – это уравнения, в которых неизвестными являются не переменные (т.е. числа), а функции одной или нескольких переменных. Эти уравнения (или системы) включают соотношения между искомыми функциями и их производными. Если в уравнения входят производные только по одной переменной, то они называются обыкновенными дифференциальными. В противном случае говорят об уравнениях в частных производных. Таким образом, решить (иногда говорят проинтегрировать) дифференциальное уравнение – значит, определить неизвестную функцию на определенном интервале изменения ее переменных.
Обыкновенные дифференциальные уравнения (ОДУ) с неизвестной функцией y(t), в которое входят производные этой функции вплоть до y(N)(t), называется ОДУ N-го порядка. В частности, уравнение первого порядка может по определению содержать помимо самой искомой функции y(t) только ее первую производную у'(t), второго порядка – у'(t) и у»(t) и т.д. В подавляющем большинстве случаев дифференциальное уравнение можно записать в стандартной форме (форме Коши): у'(t) = f(y(t), t). Уравнение второго порядка может содержать, помимо самой функции, ее первую и вторую производные и т.д.
Имеются два типа задач, которые возможно решать с помощью MathCAD:
- задачи Коши – для которых заданы значения функций в начальной точке интервала интегрирования уравнения;
- краевые задачи для которых заданы определенные соотношения сразу на обоих границах интервала.
MathCAD умеет решать только такие системы дифференциальных уравнений, которые могут быть представлены в стандартной форме (форме Коши). Для неизвестных функций y0(t) ,y1(t),…..yN-1(t) система ОДУ должна быть записана в форме:
у0‘(t) = f0(y0(t) ,y1(t),…..yN-1(t), t);
у1‘(t) = f0(y0(t) ,y1(t),…..yN-1(t), t);
…………….. (1)
уN-1‘(t) = f0(y0(t) ,y1(t),…..yN-1(t), t),
что эквивалентно следующему векторному представлению:
Y'(t) = F(Y(t), t), (2)
где Y и Y’ – соответствующие неизвестные векторные функции переменой t размера Nx1, F – векторная функция того же размера и количества переменных (N+1) (N компонент вектора и, возможно, t). Именно векторное представление (2) используется для ввода системы ОДУ в среде MathCAD.
Для того чтобы определить задачу Коши для системы из N ОДУ первого порядка, следует определить еще ровно N начальных условий, задающих значение каждой из функций yi(t0) в начальной точке интегрирования системы t0. В векторной форме они могут быть записаны в виде:
Y(t0) = B, (3)
где В – вектор начальных условий размера Nx1, составленный из yi(t0).
Стандартные процедуры MathCAD применимы для систем ОДУ первого порядка, записанных в форме (2)–(3).
Дифференциальное уравнение N-го порядка
Для решения ОДУ порядка N≥1 в MathCAD предусмотрены две возможности:
- вычислительный блок Given/Odesolve – решение имеет вид функции от t;
- встроенные функции решения систем ОДУ, причем уравнения высших порядков необходимо предварительно свести к эквивалентной системе уравнений первого порядка – в этом случае решение имеет формат вектора.
Первый путь предпочтительнее из соображений наглядности представления задачи и результатов, второй дает пользователю больше рычагов воздействия на параметры численного метода.
Вычислительный блок для решения ОДУ, реализующий численный метод Рунге–Кутта, состоит из трех частей:
Given – ключевое слово;
ОДУ и начальные условия в формате y(t0) = b, записанные с помощью логических операторов, которые должны набираться на панели инструментов Boolean (Булевы операторы);
Odesolve(t, t1) – встроенная функция для решения ОДУ относительно переменной t на интервале (t0, t1), причем t0<t1.
На запись ОДУ в пределах вычислительного блока накладывается несколько
ограничений. Во-первых, ОДУ должно быть линейно относительно старшей производной, т.е. фактически должно быть поставлено в стандартной форме. Во-вторых, начальные условия должны иметь форму
у(t0) = b или у(N)(t0) = b.
Рис. 9.1. Пример решения задачи Коши
Предпочтительнее задание функции Odesolve(t, t1, step) с тремя параметрами, где step – внутренний параметр численного метода, определяющий количество шагов, в которых метод Рунге–Кутта будет рассчитывать решение дифференциального уравнения. Чем больше step, тем с лучшей точностью будет получен результат, но тем больше времени будет затрачено на его поиск. Помните, что подбором этого параметра можно заметно ускорить расчеты без существенного ухудшения их точности.
Пример решения задачи Коши для ОДУ второго порядка – системы генератора с затуханием (с параметром квадратичной нелинейности γ), приведен на рис. 9,1.
Символ производной допускается вводить как средствами панели Calculus (Вычисления), так и в виде штриха, набрав его с помощью сочетания клавиш <Ctrl>+<F7>.
Результатом применения блока Given/Odesolve является функция y(t), определенная на промежутке (t0, t1). Следует воспользоваться обычными средствами MathCAD, чтобы построить ее график или получить значение функции в какой-либо точке указанного интервала.
Пользователь имеет возможность выбирать между двумя модификации численного метода Рунге–Кутта. Для смены метода необходимо нажат правой кнопки мыши на области функции Odesolve вызвать контекстное меню и выбрать в нем один из трех пунктов: Fixed (С фиксированным шагом), Adaptive (Адаптивный) или Stiff (Для жестких ОДУ).
Система N дифференциальных уравнений
При помощи MathCAD можно решать системы N>1 ОДУ первого порядка, если они записаны в стандартной форме (Коши) в виде векторного соотношении Y'(t)=F,(Y(t), t).
В MathCAD имеется несколько встроенных функций, которые позволяют решать задачу Коши различными численными методами. Для "хороших" нежестких систем ОДУ применяются следующие функции:
rkfixed(y0, t0, t1, M, D) – метод Рунге–Кутта с фиксированным шагом;
Rkadapt(y0, t0, t1, M, D) – метод Рунге–Кутта с переменным шагом;
Bulstoer(y0, t0, t1, M, D) –метод Булирша–Штера:
у0 – вектор начальных значений в точке to размера Nx1;
t0 – начальная точка расчета;
t1 – конечная точка расчета;
М – число шагов, на которых численный метод находит решение;
D – векторная функция размера Nx1 двух аргументов – скалярного t и векторного у. При этом у – искомая векторная функция аргумента t того же размера Nx1.
Необходимо соблюдать регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций MathCAD, например, Find ≡ find.
Каждая из приведенных функций выдает решение в виде матрицы размера (М+1)x(N+1). В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных N столбцах – значения искомых функций y0(t), y1(t), …, yN-1(t), рассчитанные для этих значений аргумента. Поскольку всего точек
(помимо начальной) М, то строк в матрице решения будет всего М+1.
На рис. 56 приведен пример решения системы ОДУ генератора с затуханием при помощи функции rkfixed, результат расчетов представлен как содержимое матрицы и в виде графика. Точки, в которых получено решение, отмечены на графике кружками. Чтобы использовать другой численным алгоритм, достаточно поменять имя функции rkfixed в последней строке программы на другое (на практике более эффективны функции Rkadapt и Bulstoer).
Первая строка листинга представляет задание параметров модели, вторая – начального условия задачи Коши, а в третьей строке листинга определено число шагов, на которых рассчитывается решение. Самая важная – это предпоследняя строка листинга, в которой, собственно, определяется система ОДУ. Последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (0,40).
Функция D, входящая в число параметров встроенных функций для решения ОДУ, должна быть функцией обязательно двух аргументов. Второй ее аргумент должен быть вектором того же размера, что и сама функция D. Точно такой же размер должен быть и у вектора начальных значений у0. Векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки Subscript (Нижний индекс) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>.
Рис. 9.2. Пример решения
Размер полученной матрицы, представляющей решение, равен (M+1)x(N+1), т.е. 51х3. Просмотреть все компоненты матрицы u, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки.
Рис. 9.3. Уравнение
Решение одного уравнения (N=1)
Вычислительный блок Given/Odesolve выигрывает в простоте и в наглядности, однако иногда предпочтительнее решать ОДУ первого порядка с помощью встроенных функций rkfixed, Rkadapt или Bulstoer, например, при следующих обстоятельствах:
- одно ОДУ решается в контексте решения более сложных задач, в которые входят системы дифференциальных уравнений – в этом случае может потребоваться единый стиль программирования;
- ответ предпочтительнее получить в виде вектора, а не функции и т.п.
При решении одного ОДУ и само уравнение, и начальное условие можно задавать в скалярной форме (рис. 57). Результат выдается в виде матрицы размерности Мх2, которая состоит из двух столбцов: в одном находятся значения аргумента t (от t0 до t1 включительно), в другом – соответствующие значения искомой функции у(t).
Решение систем ОДУ в заданной точке
При решении дифференциальных уравнений требуется определить значения искомых функций не на всем интервале (t0, t1), а только в одной его последней точке. Известно, что для широкого класса ОДУ одна и та же система при разных начальных условиях при t → ∞ приходит в одну и ту же точку. Поэтому часто нужно определить именно эту точку.
Такая задача требует меньше ресурсов компьютера, чем решение системы ОДУ на всем интервале, поэтому в MathCAD имеются модификации встроенных функций Rkadapt и Bulstoer. Они имеют несколько другой набор параметров и работают значительно быстрее своих аналогов:
rkadapt(у0, t0, t1, асс, D, k, s) – метод Рунге–Кутта с переменным шагом;
bulstoer(у0, t0, t1, асс, D, k, s) – метод Булирша–Штера:
у0 – вектор начальных значений в точке to;
t0, t1 – начальная и конечная точки расчета;
асc – погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение; рекомендуется выбирать значения погрешности в районе 0.001);
d – векторная функция, задающая систему ОДУ;
k – максимальное число шагов, на которых численный метод будет находить решение;
s – минимально допустимая величина шага.
Вместо числа шагов на интервале интегрирования ОДУ в этих функциях необходимо задать точность расчета численным методом значения функций в последней точке. В этом смысле параметр асс похож на константу TOL, которая влияет на большинство встроенных численных алгоритмов MathCAD. Количество шагов и их расположение определяются численным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственно повлиять на разбиение интервала на шаги. Параметр k служит для того, чтобы шагов не было чрезмерно много, причем нельзя сделать k >1000. Параметр s – для того, чтобы ни одни шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма. Эти параметры следует задавать явно, исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая.
Рис. 9.4. Тестовый расчет
Пример использования функции bulstoer для той же модели линейного генератора приведен на рис. 58. В его первых двух строках, как обычно, определяется система уравнений и начальные условия; в следующей строке матрице и присваивается решение, полученное с помощью bulstoer. Структура этой матрицы в точности такая же, как и в случае решения системы ОДУ посредством уже рассмотренных нами ранее встроенных функций, однако в данном случае нам интересна только последняя точка интервала. Поскольку сделанное численным методом количество шагов, размер матрицы u, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строке листинга, присваивающей это число
переменной М, в этой же строке оно выведено на экран.
В предпоследней программе осуществлен вывод решения системы ОДУ на конце интервала, т. е. в точке t=50 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки). Чтобы попробовать альтернативный численный метод, достаточно заменить имя функции bulstoer на rkadapt.
Функции bulstoer и rkadapt (те, что пишутся со строчной буквы) не предназначены для нахождения решения в промежуточных точках интервала, хотя они и выдают их в матрице-результате.
О численных методах
Все численные методы решения ОДУ основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации получаются алгоритмы различной точности и быстродействия. В MathCAD использован наиболее популярный алгоритм Рунге–Кутта четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, тем более что сделать это очень просто благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.
Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо либо существуют участки медленных и быстрых его изменений. Метод Рунге–Кутта с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений – частыми. В результате для достижения одинаковой точности требуется меньшее число шагов, чем для rkfixed. Метод Булирша–Штера (Bulstoer) часто оказывается более эффективным для поиска гладких решений.
Имея в виду сделанные замечания, приведем краткую сводку алгоритмов решения задач Коши для ОДУ, отмечая, какие из встроенных функций следует использовать в конкретных случаях:
для решения единственного уравнения (любого порядка) используйте
- вычислительный блок Given/Odesolve;
- для стандартных нежестких систем используйте алгоритм Булирша– Штера (Bulstoer);
- для систем с участками быстро и медленно меняющихся решений используйте адаптивный алгоритм Рунге–Кутта (Rkadapt);
- в учебных целях и для решения несложных задач можно использовать алгоритм Рунге–Кутта с фиксированным шагом (Rkfixed);
- для получения решения в одной конечной точке интервала используйте (в зависимости от перечисленных классов задач) одну из встроенных функций с именем, начинающимся со строчной буквы.
Тем не менее, в MathCAD можно реализовать любой численный метод. Обычно решение заключается в нахождении ряда значений хi и уi искомой зависимости у(х) при i, изменяющемся от 0 до N при шаге изменения х, равном h. Рассмотрим три наиболее распространенных способа решения дифференциального уравнения, при которых h = const. Значения х0 и у0 должны быть известны как начальные условия, без них невозможно единственное решение.
Аналитическое решение обыкновенного дифференциального уравнения имеет вид: у(х)=А·ехр(х2/2). Погрешность метода определяется выражением: hn+1, где n – порядок метода. Простой метод Эйлера имеет первый порядок, при этом погрешность метода составляет h2. Будучи самым простым из известных, этот метод имеет и наибольшую погрешность. При малых h и достаточно гладких решениях этот метод вполне приемлем для решения многих практических задач.
Простой метод Эйлера реализуется применением на каждом шаге вычислений следующих итерационных выражений:
Для уменьшения погрешности решения следует применять методы более высокого порядка. Модифицированный метод Эйлера является методом второго порядка, и его погрешность пропорциональна h3.
Данный метод реализуется следующими итерационными выражениями на каждом шаге вычислений:
При высоких требованиях к точности решения можно воспользоваться методом Рунге-Кутта четвертого порядка. При нем погрешность пропорциональна h5. Этот метод реализуется с помощью следующих используемых на каждом шаге вычислений выражений:
Решение дифференциального уравнения второго порядка рассмотрим методом Рунге-Кутта четвертого порядка. Итерационные уравнения имеют вид:
Решить систему из двух дифференциальных уравнений простым методом Эйлера можно по итерационным уравнениям в векторной форме:
Модифицированный метод Эйлера имеет следующие уравнения:
Результаты решения необходимо приводить в двух видах: в табличном виде или в виде графика. Для методов Эйлера графическое решение приводят в сравнении с аналитическим решением.