Лабораторная работа 9

 

Численное решение
обыкновенных дифференциальных уравнений

(задача Коши)

 

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

 

Количество аудиторных часов, отводимых на работу: 4

 

 

Перед выполнением работы изучите материал первого и  второго параграфа шестой главы данного пособия и ответьте на вопросы для самоконтроля (1-6). Расчетные программы задания 3 могут быть составлены на любом алгоритмическом языке (Бейсик, Паскаль, Си, Фортран).

 

Решение обыкновенных дифференциальных уравнений в Mathcad

(краткие теоретические сведения)

 

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

Использование встроенных функций и самостоятельную реализацию алгоритмов численных методов рассмотрим на примере решения простейшего уравнения:

,                                               (1)

на отрезке  с начальным условием .

Точное решение этого уравнения известно:­ .

Решение уравнения (1) методом Эйлера можно получить непосредственно по формуле (5.8). Пример решения приведен ниже.

 

 

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

 

Перейдем к рассмотрению некоторых встроенных в Mathcad функций для решения обыкновенных дифференциальных уравнений.

Функция Odesolve

Функция Odesolve([vector], x, b, [nstep]) – возвращает решение обыкновенного дифференциального уравнения с заданными начальными (задача Коши) или граничными (краевая задача) условиями. Функция Odesolve используется совместно с ключевым словом Given, организующим «вычислительный блок».

Параметры:

·        vector – вектор искомых функций (задается только в случае решения систем дифференциальных уравнений);

·        x – независимая переменная;

·        b – конечная точка решения (начальная точка указывается внутри вычислительного блока);

·        nstep – необязательный параметр, задающий количество шагов интегрирования (по умолчанию используется значение 1000).

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

Пример решения уравнения (1) с помощью функции Odesolve:

Ø Замечание. При записи дифференциального уравнения и начальных условий внутри блока вычислений (Given …) необходимо использовать знак логического равенства (Ctrl+=).<

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

Численный метод, по которому работает функция Odesolve, можно выбирать из контекстного меню. В зависимости от версии пакета Mathcad это может быть один из методов Адамса или Рунге-Кутта 4-го порядка с фиксированным и переменным шагом, а также другие методы.

Функция rkfixed

Функция rkfixed(ic,a,b,nstep,D) – возвращает решение обыкновенного дифференциального уравнения (или системы уравнений) первого порядка  с заданными начальными условиями. Функция использует метод Рунге-Кутта 4-го порядка точности с фиксированным шагом.

Параметры:

·        ic –вектор начальных условий;

·        а – начальное значение независимой переменной;

·        b – конечное  значение независимой переменной;

·        nstep –параметр, задающий количество шагов интегрирования;

·        Dвектор правых частей дифференциальных уравнений.

Пример решения уравнения (1) с помощью функции rkfixed:

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

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

Функция Rkadapt

Функция Rkadapt (ic,a,b,nstep,D) – возвращает решение обыкновенного дифференциального уравнения (или системы уравнений) первого порядка  с заданными начальными условиями. Функция использует адаптивный метод Рунге-Кутта 4-го порядка точности с изменяющимся шагом и по своему использованию полностью аналогична функции rkfixed.

 

 

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

 

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

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

,                                         (2)

с начальными условиями ,  на интервале .

 

 

 

 

 

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

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

                                                 (3)

Начальные условия для системы (3) можно переписать в виде:

, .

Проиллюстрируем решение системы дифференциальных уравнений с помощью функции Odesolve на примере системы (3).

 

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

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

Поясним смысл основных действий. После определения параметров задачи, задается вектор начальных условий IC. Элементы этого вектора – значения искомых функций в начальной точке . Далее определяется функция-вектор D(t,Y), содержащая правые части дифференциальных уравнений. Элементы вектора Y – это искомые функции. В данном случае  обозначает функцию y, а  – функцию v. После этого вызывается функция rkfixed. При вызове функции указывается вектор начальных условий – IC, начальная – 0 и конечная – 1 точки отрезка, на котором ищется решение, количество шагов интегрирования – 100 и вектор правых частей дифференциальных уравнений – D. Результат записывается в матрицу S. Первый (нулевой) столбец этой матрицы содержит значения независимой переменной t, второй и третий столбцы – значения искомых функций  и , соответственно.

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

 

Практические задания

 

Задание 1. Решить обыкновенное дифференциальное уравнение первого порядка методом Эйлера, модифицированным методом Эйлера и методом Рунге-Кутта 4-го порядка с постоянным шагом. Оценить погрешность разных методов и исследовать влияние шага интегрирование на величину погрешности. Уравнение и начальные условия возьмите из таблицы 1.

 

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

Для оценки погрешности проделайте следующее (для каждого метода):

– найдите общее решение дифференциального уравнения аналитическим методом (интегрированием);

– найдите частное решение, соответствующее заданным начальным условиям;

– определите в документе Mathcad функцию, соответствующую аналитическому (точному) решению, например yt(x) := …;

– вычислите погрешность приближенного решения, например так: erri := yt(xi) – yi, где yi – приближенное решение, а xi – узлы, в которых найдено приближенное решение уравнения;

– постройте графики erri как функции от xi для разных значений шагов интегрирования (например, h = 0.1, 0.05 и 0.01).

 

 

Задание 2. Решить обыкновенное дифференциальное уравнение второго порядка методом Эйлера и методом Рунге-Кутта 4-го порядка с постоянным шагом. Уравнение и начальные условия возьмите из таблицы 2.

 

Пояснения к выполнению задания. Решите заданное уравнение с помощью функции Odesolve. Преобразуйте заданное уравнение к системе двух дифференциальных уравнений первого порядка. Решение методом Эйлера выполните аналогично разобранному в теоретической части примеру. Решение методом Рунге-Кутта выполните с помощью встроенной в Mathcad функции rkfixed. Получите решения методами Эйлера и Рунге-Кутта при различных значениях шага интегрирования (например, h = 0.1, 0.05 и 0.01). Сравните решения при одинаковых значениях шага.

 

Задание 3. Составьте программы для решения дифференциального уравнения второго порядка методами Эйлера и Рунге-Кутта. Приближенные значения искомой функции запишите в текстовый файл. Графически изобразите решение с помощью программы Origin. Уравнение и начальные условия возьмите из таблицы 2.

 

Таблица 1.

Вариант

Уравнение

Отрезок

Начальные условия

1

 

 

 

2

 

 

 

3

 

 

 

4

 

 

 

5

 

 

 

6

 

 

 

7

 

 

 

8

 

 

 

9

 

 

 

10

 

 

 

11

 

 

 

12

 

 

 

13

 

 

 

14

 

 

 

15

 

 

 

16

 

 

 

17

 

 

 

18

 

 

 

19

 

 

 

20

 

 

 

 

 

Таблица 2.

Вариант

Уравнение

Отрезок

Начальные условия

1

2

3

4

5

6

7

 

 

 

8

 

 

 

9

 

 

 

10

 

 

 

11

 

 

 

12

 

 

 

13

 

 

 

14

 

 

 

15

 

 

 

16

 

 

 

17

 

 

 

18

 

 

 

19

 

 

 

20