- •1. ЭЛЕМЕНТАРНАЯ ТЕОРИЯ ПОГРЕШНОСТЕЙ; ВЫЧИСЛИТЕЛЬНЫЕ ЗАДАЧИ, МЕТОДЫ И АЛГОРИТМЫ
- •1.1. Источники и классификация погрешностей результата численного эксперимента
- •1.2. Погрешности чисел
- •1.3. Погрешности арифметических операций
- •1.4. Погрешности функций
- •1.5. Особенности машинной арифметики
- •1.6. Лабораторная работа № 1. Определение абсолютной и относительной погрешностей приближенных чисел. Оценка погрешностей результата
- •1.7. Корректность вычислительной задачи
- •1.8. Обусловленность вычислительной задачи
- •1.9. Вычислительные методы, их классификация
- •2. ПРИБЛИЖЕНИЕ ФУНКЦИЙ
- •2.1. Задача приближения функций
- •2.2. Интерполяция обобщенными многочленами
- •2.3. Полиномиальная интерполяция. Многочлен Лагранжа
- •2.4. Погрешность интерполяции
- •2.5. Конечные разности и их свойства
- •Доказательство
- •2.6. Разделенные разности и их свойства
- •2.9. Лабораторная работа № 2. Интерполирование и экстраполирование данных. Интерполяционный многочлен Лагранжа
- •2.10. Интерполяционный многочлен Ньютона с конечными разностями
- •2.11. Лабораторная работа № 3. Интерполирование и экстраполирование данных. Интерполяционный многочлен Ньютона
- •2.12. Интерполяционные формулы Гаусса, Стирлинга и Бесселя
- •3. МЕТОД НАИМЕНЬШИХ КВАДРАТОВ И СПЕЦИАЛЬНЫЕ ИНТЕРПОЛЯЦИОННЫЕ МНОГОЧЛЕНЫ
- •3.1. Постановка задачи и вывод формул метода наименьших квадратов
- •3.3. Глобальная полиномиальная интерполяция
- •3.4. Чувствительность интерполяционного многочлена к погрешностям входных данных
- •3.5. Многочлены Чебышева
- •3.6. Решение задачи минимизации оценки погрешности
- •3.8. Лабораторная работа №5. Экономизация степенных рядов
- •3.9. Локальная интерполяция
- •3.10. Сплайны, их свойства и построение
- •3.11. Погрешность приближения кубическими сплайнами
- •3.13. Тригонометрическая интерполяция. Дискретное преобразование Фурье и его реализация на ЭВМ
- •3.14. Матричная форма записи дискретного преобразования Фурье (ДПФ)
- •3.15. Алгоритм реализации ДПФ
- •3.16. Пример реализации алгоритма ДПФ при
- •3.17. Лабораторная работа № 7. Дискретное преобразование Фурье
- •4. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
- •4.1. Простейшие формулы численного дифференцирования для первой производной
- •4.2. Формулы численного дифференцирования для второй производной
- •4.3. Формулы численного дифференцирования, основанные на интерполяции алгебраическими многочленами
- •4.4. Обусловленность формул численного дифференцирования
- •4.5. Простейшие квадратурные методы численного интегрирования
- •4.6. Оценка погрешностей простейших квадратурных формул
- •4.7. Квадратурные формулы интерполяционного типа
- •4.8. Квадратурные формулы Гаусса
- •4.9. Лабораторная работа № 8. Численное дифференцирование и численное интегрирование функций
- •5. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ И ПРОБЛЕМЫ СОБСТВЕННЫХ ЗНАЧЕНИЙ
- •5.1. Нормы векторов и матриц и их свойства
- •5.2. Обусловленность задачи решения системы линейных алгебраических уравнений
- •5.3. Метод Гаусса (схема единственного деления)
- •5.4. Метод прогонки
- •5.5. Метод простых итераций
- •5.6. Сходимость метода простых итераций
- •5.10. Постановка задачи нахождения собственных чисел
- •5.11. Подобные матрицы
- •5.12. Локализация собственных значений
- •5.13. Степенной метод
- •5.14. Вычисление собственных векторов методом обратных итераций
- •6. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ И СИСТЕМ
- •6.1. Решение нелинейных уравнений
- •6.2. Метод Ньютона для уравнений
- •6.3. Сходимость метода Ньютона и трудности его применения
- •6.4. Метод Ньютона решения систем нелинейных уравнений
- •6.6. Модификации метода Ньютона
- •6.7. Лабораторная работа № 11. Решение систем нелинейных уравнений методом Ньютона
- •7. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
- •7.1. Задача Коши для дифференциального уравнения первого порядка
- •7.2. Численные методы решения задачи Коши. Основные понятия и определения
- •7.3. Решение с помощью рядов Тейлора
- •7.5. Анализ ошибок, возникающих при использовании методов Рунге - Кутты
- •7.6. Методы прогноза и коррекции
- •7.7. Сравнение методов
- •7.8. Лабораторная работа № 12. Методы интегрирования обыкновенных дифференциальных уравнений
- •7.9. Решение задачи Коши для систем обыкновенных дифференциальных уравнений
- •7.11. Лабораторная работа № 13. Численное интегрирование систем дифференциальных уравнений первого порядка
- •8. ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ (УРАВНЕНИЯ МАТЕМАТИЧЕСКОЙ ФИЗИКИ)
- •8.1. Классификация уравнений математической физики
- •8.2. Простейшие задачи, приводящие к дифференциальным уравнениям в частных производных
- •8.4. Уравнения параболического типа. Явные и неявные схемы
- •Доказательство
- •8.5. Уравнения гиперболического типа
- •8.6. Уравнения эллиптического типа
- •8.7. Свойства разностных схем для дифференциальных уравнений: способность аппроксимировать исходную дифференциальную задачу, устойчивость и сходимость
- •8.8. Некоторые обобщения
- •8.9. Лабораторная работа № 14. Решение задачи Дирихле для уравнения Лапласа методом сеток
- •8.10. Лабораторная работа № 15. Решение однородного уравнения колебаний струны методом сеток по неявной схеме.
7.5. Анализ ошибок, возникающих при использовании методов Рунге - Кутты
Один из серьезных недостатков методов Рунге - Кутты состоит в отсутствии простых способов оценки ошибки интегрирования. Все же без некоторой оценки ошибки трудно правильно выбрать величину шага интегрирования h. Применим так называемый принцип Рунге.
Пусть yn(T ) есть точное решение дифференциального уравнения y/ = f (x, y) при x = x0 + nh. Тогда для метода Рунге - Кутты, описываемого формулами (7.4.13), справедлива следующая оценка погрешности:
|
|
|
|
h |
|
|
|
|
h |
|
|
|
|
|
(h) |
|
|
|
2 |
|
|
|||||
|
|
|
|
|
||||
|
|
|
||||||
ε 2 |
|
= |
yn |
|
− yn |
, |
||
|
2 p −1 |
|||||||
|
|
|
|
|
где yn(h ) - приближение к точному решению yn(T ), вычисленное с шагом
|
h |
h |
|
h |
|
||
|
|
|
|
|
|
|
|
|
|
|
|||||
приближение с шагом |
, ε |
2 |
|
= yn(T ) − yn |
2 |
. |
|
2 |
|
|
|
||||
|
|
|
|
|
|
|
(7.5.1)
h
h , yn 2 - такое же
Так как для метода, описываемого формулами (7.4.13), p = 4, то
h |
|
|
1 |
|
h |
|
|
|
|||
|
|
|
|
|
|
|
|
||||
|
|
|
|
||||||||
ε |
2 |
|
= |
|
|
yn |
2 |
|
− yn(h) . |
(7.5.2) |
|
|
|
15 |
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
Формула (7.5.1) выведена в предположении, что на каждом шаге интегрирования допускается погрешность, приблизительно пропорциональная h p+1 , то есть yn(T ) = yn(h ) + Kh5 , что справедливо для достаточно гладких функций. Таким образом, ошибка интегрирования в предположении, что y(5)(x) практически постоянна, равна:
h |
|
|
|
1 |
|
h |
|
|
|||
|
|
|
|
|
|
|
|
||||
|
|
|
|
||||||||
ε |
2 |
|
= Kh5 |
= |
|
|
y |
2 |
|
− y(h) . |
|
|
|
15 |
|
|
|||||||
|
|
|
|
|
|
n |
|
|
n |
||
|
|
|
|
|
|
|
|
|
|
|
|
Это довольно точная оценка, однако для ее использования необходимо вычислять решение дважды.
Предложено несколько полуэмпирических критериев смены шага и выбора оптимального шага интегрирования при условии достижения заданной точности. Например, использу-
kn(2) − kn(3)
ется такое оценочное правило: если ( ) ( ) достаточно велико (обычно больше несколь- kn1 − kn2
ких сотых), то шаг интегрирования необходимо уменьшить.
Существуют более точные методы оценки погрешности интегрирования, основанные на использовании для контроля точности двух различных методов Рунге - Кутты. Один из самых эффективных - метод Рунге - Кутты - Фельберга. В этом методе для оценки погрешности метода пятого порядка используются формулы метода четвертого порядка точности, причем на одном шаге интегрирования требуется всего лишь шесть вычислений значений правой части f (x, y).
7.6. Методы прогноза и коррекции
Отличительной чертой методов Рунге - Кутты является то, что при вычислении следующей точки (xn+1 , yn+1 ) используется информация только о предыдущей (xn , yn ). Зато при-
ходится вычислять значение правой части в нескольких промежуточных точках. Это нерационально, поскольку при налаженном процессе интегрирования имеется дополнительная информация - предыдущие точки решения.
169
Рассмотрим общую идею группы методов, известных в литературе под названием «методов прогноза и коррекции». Как ясно из названия, вначале «предсказывается» значение yn+1 , а затем используется тот или иной метод для «корректировки» этого значения. Эту же
формулу можно использовать сколько угодно раз для повторной корректировки уже скорректированного значения yn+1 .
Для демонстрации основных идей метода можно использовать для прогноза формулу
любого метода численного интегрирования. Воспользуемся формулой второго порядка |
|
yn(0+)1 = yn−1 + 2hf (xn , yn ), |
(7.6.1) |
где индекс (0) означает исходное приближение к yn+1 . Непосредственно из написанной формулы следует, что с ее помощью нельзя вычислить, например, y1 , ибо для этого потребовалась бы точка, расположенная перед начальной точкой y0 . Чтобы начать решение по методу
прогноза и коррекции, часто используется метод Рунге - Кутты.
Геометрически предсказание сводится к тому, что находится угол наклона касательной в точке (xn , yn ) (см. левый рисунок). После этого через точку (xn−1 , yn−1 ) проводится
|
|
|
L |
|
L |
L |
L |
|
|
|
1 |
|
2 |
1 |
L3 |
|
|
|
|
|
|
|
|
|
|
|
(xn+1 , yn(0+)1 ) |
|
|
(xn+1 , yn(1+)1 ) |
|
|
L |
|
|
|
|
||
|
y = f (x) |
|
|
|
y = f (x) |
|
|
|
h |
h |
|
h |
h |
|
|
xn−1 |
xn |
|
|
xn+1 |
xn−1 |
xn |
xn+1 |
прямая L , параллельная |
L |
. Предсказанное значение y(0) |
будет расположено там, где пря- |
||||
|
|
|
1 |
|
n+1 |
|
|
мая L |
пересечется с абсциссой x = xn+1. |
Скорректируем теперь предсказанное значение. |
Вычислим наклон касательной в точке (xn+1 , yn(0+)1 ). Эта касательная изображена на правом ри-
сунке и обозначена L2 . Усредним тангенсы углов наклона прямых L1 и L2 |
и получим пря- |
||||||
мую L . Наконец, проведем через точку (x |
n |
, y |
n |
) прямую L , параллельную L , и точка пере- |
|||
|
|
|
|
3 |
|
||
сечения этой прямой с абсциссой x = xn+1 даст новое приближение (xn+1 , yn(1+)1 ). |
|||||||
Это значение вычисляется по формуле |
|
|
|
||||
yn(1+)1 = yn + |
h |
[f (xn , yn )+ f (xn+1 , yn(0+)1 )]. |
(7.6.2) |
||||
|
|||||||
2 |
|
|
|
|
|
|
В общем случае i -е приближение, очевидно, будет находиться таким образом:
yn(i+)1 |
= yn + h [f (xn , yn )+ f (xn+1 , yn(i+−11))]. |
(7.6.3) |
|||
|
2 |
|
|
|
|
Итерационный процесс можно прекратить, когда |
|
||||
|
|
yn(i++11) − yn(i+)1 |
|
< ε. |
(7.6.4) |
|
|
|
К настоящему времени разработано много методов прогноза и коррекции разных порядков, обеспечивающих очень большую точность решения. В некоторых из них по заданной точностипредусмотренавтоматическийвыбор шага и автоматический выбор порядка метода.
170
7.7. Сравнение методов
а). Методы Рунге - Кутты.
1.В этих методах используется информация только об очередной точке решения, поэтому с их помощью можно начинать решение.
2.По той же самой причине приходится многократно вычислять функцию f (x, y).
3.Методы этой группы позволяют очень легко менять величину шага h.
4.При использовании этих методов трудно получить оценку ошибки интегрирования. б). Методы прогноза и коррекции.
1. Так как в этих методах используется информация о ранее вычисленных точках решения, то с их помощью нельзя начать решение.
2. Из-за использования информации о ранее вычисленных точках методы этой группы более экономичны по затратам машинного времени.
3. При любом изменении величины шага h приходится временно возвращаться к методам Рунге - Кутты.
4. В качестве побочного продукта вычислений получается хорошая оценка ошибки интегрирования.
Как всегда, наиболее целесообразным является использование при решении практи-
ческих задач комбинации этих двух методов. |
|
Пример. Методом Рунге - Кутты с шагом h = 0.1 найти на отрезке [0, 0.3] решение |
|
следующего дифференциального уравнения: y/ = cos bx |
с начальным условием |
a + y2 |
|
y(0)= 0, a =1.4, b = 2.6.
|
Обозначим через yi |
приближенное значение решения в точке xi . Формулы метода: |
||||||||||||||||||||||||||||||||||||
|
|
|
|
|
y |
i+1 |
= y |
i |
+ hk |
, |
k |
i |
= |
1 (k (1) + |
2k (2) + |
2k (3) + k (4)), |
||||||||||||||||||||||
|
|
|
|
|
|
|
|
i |
|
|
|
|
|
|
|
6 |
i |
|
|
|
|
|
|
i |
|
i |
i |
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
k (1) |
|
|
|
|
|
|
), |
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
= f (x |
, |
y |
|
|
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
i |
|
|
|
|
i |
|
|
i |
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
k (2) |
= f |
x |
i |
+ h |
, y |
i |
+ h k (1) , |
(7.7.1) |
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
i |
|
|
|
|
|
|
|
2 |
|
|
|
|
|
2 |
i |
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
h |
|
|
|
|
|
|
h |
(2) |
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
(3) |
= f |
|
|
|
, yi |
+ |
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
ki |
|
xi + |
2 |
2 |
ki |
, |
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
(4) |
|
|
|
|
|
|
|
|
|
|
|
|
(3) |
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
ki |
|
= f (xi + h, yi + hki |
). |
|
|
|
|
|
||||||||||||||||||
|
Все вычисления удобно располагать по схеме, указанной в таблице. |
|||||||||||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
i |
|
|
|
x |
|
|
|
|
|
|
|
|
y |
|
|
|
|
|
|
|
|
|
|
|
|
|
k |
|
|
y |
|
|
|||||
|
0 |
|
|
|
x0 |
|
|
|
|
|
|
|
|
y0 |
|
|
|
|
|
|
|
|
|
|
k0(1) |
|
|
k0(1) |
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
x0 |
+ |
h |
|
|
|
|
|
y0 |
+ |
|
h |
k0(1) |
|
|
|
|
|
|
|
|
|
k0(2) |
|
|
2k0(2) |
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||||||
|
|
|
|
|
2 |
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
x0 |
+ |
h |
|
|
|
|
|
y0 |
+ |
h |
k0(2) |
|
|
|
|
|
|
|
|
|
k0(3) |
|
|
2k0(3) |
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
|
|
|
|
|
2 |
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||||
|
|
|
|
x0 + h |
|
|
|
|
y0 + hk0(3) |
|
|
|
|
|
|
|
|
|
k0(4) |
|
|
k0(4) |
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
y0 |
1 |
x1 |
y1 |
Она заполняется следующим образом.
1. Записываем в первой строке x0 , y0 , вычисляем f (x0 , y0 ) и заносим в таблицу в качестве k0(1).
171
|
2. |
Для |
второй |
строки |
вычисляем |
x |
0 |
|
+ |
h |
|
и |
y |
|
+ |
h |
k (1) , |
затем |
находим |
|||||||||||
|
|
h |
|
h |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
0 |
2 |
0 |
|
|
||||||
|
|
|
(1) |
|
|
|
|
|
(2) |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
f x0 |
+ |
|
, y0 + |
|
k0 |
|
и записываем в качестве k0 |
|
. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
2 |
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
3. |
Находим |
y0 |
+ h k0(2) , вычисляем |
f x0 |
+ |
h |
, y0 + |
h |
k0(2) |
|
и записываем в таблицу на |
||||||||||||||||||
|
2 |
2 |
||||||||||||||||||||||||||||
место k0(3) . |
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
+ hk0(3) , |
|
|
||||||
|
4. |
Записываем |
|
в четвертой строке |
|
|
x0 + h |
|
и |
y0 |
затем |
находим |
||||||||||||||||||
k0(4 ) = f (x0 + h, y0 + hk0(3)). |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
5. |
Суммируем числа, стоящие в столбце y , делим на шесть и заносим в таблицу в |
||||||||||||||||||||||||||||
качестве |
y0 . |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
6. |
Вычисляем x1 |
= x0 + h , |
y1 = y0 + |
y0 . Затем все вычисления необходимо повторить |
|||||||||||||||||||||||||
с пункта 1, принимая за начальную точку (x1, y1 ). |
Для контроля правильности выбора ша- |
|||||||||||||||||||||||||||||
га h рекомендуется вычислять дробь θ = |
|
ki(2) − ki(3) |
|
|
|
. Величина θ не должна превышать двух |
||||||||||||||||||||||||
|
|
|||||||||||||||||||||||||||||
|
|
|
|
|
||||||||||||||||||||||||||
|
ki(1) − ki(2) |
|
|
|||||||||||||||||||||||||||
|
|
|
- трех сотых. В противном случае шаг h следует уменьшить.
i |
|
x |
|
|
y |
|
|
k |
|
|
|
y |
|
|
θ |
|
0 |
|
0.00 |
|
|
0.000000 |
|
|
0.714286 |
|
|
|
0.14286 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.05 |
|
0.035714 |
|
0.707614 |
|
1.415228 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.05 |
|
0.035381 |
|
0.707626 |
|
1.415252 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.10 |
|
0.070763 |
|
0.687818 |
|
0.687818 |
0.0018 |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.705431 |
|
|
|
|
||
1 |
|
0.10 |
|
|
0.705431 |
|
|
0.509261 |
|
|
|
0.509262 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.15 |
|
0.730894 |
|
0.478185 |
|
|
|
0.956370 |
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.15 |
|
0.729340 |
|
0.478747 |
|
0.957494 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.20 |
|
0.753306 |
|
0.441084 |
|
0.441084 |
0.0018 |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.477368 |
|
|
|
|
||
2 |
|
0.20 |
|
|
1.182799 |
|
|
0.310045 |
|
|
|
0.310045 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.25 |
|
1.198301 |
|
0.280714 |
|
0.561428 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.25 |
|
1.196835 |
|
0.281062 |
|
0.562124 |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.30 |
|
1.210905 |
|
0.248026 |
|
0.248026 |
0.0012 |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.280270 |
|
|
|
|
3 |
|
0.30 |
|
|
1.463069 |
|
|
|
|
|
|
|
|
|
|
|
7.8. Лабораторная работа № 12. Методы интегрирования обыкновенных дифференциальных уравнений
Существует большое число методов приближенного решения дифференциальных уравнений, основанных на самых различных идеях. Численные методы дают приближенное решение y(x) в виде таблицы значений y1 , y2 ,..., yn в точках x1 , x2 ,..., xn .
Простейшим методом численного интегрирования дифференциального уравнения первого порядка y/ = f (x, y), удовлетворяющего начальному условию y(x0 )= y0 , является метод Эйлера. В нем величины yi вычисляются по формуле
xi = x0 + ih, i = 0,1,2,...
yi+1 = yi + hf (xi , yi ).
172
Метод Эйлера относится к группе одношаговых, в которых для расчета точки (xi+1 , yi+1 ) требуется информация только о последней вычисленной точке (xi , yi ). Геометри-
ческая интерпретация метода изложена в подразд. 7.4.
В среде Mathcad имеется тринадцать встроенных функций решения дифференциальных уравнений и систем ( задача Коши, краевая задача, уравнения в частных производных). Самая употребительная из них - rkfixed, в которую заложен метод Рунге – Кутты четвертого порядка с постоянным шагом. Подпрограммы для метода Эйлера нет из-за его низкой точности. Формулы метода Эйлера настолько просты, что вычисления по ним можно организовать
с помощью дискретной переменной. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
Рассмотрим пример. |
Решим |
|
задачу |
|
|
Коши |
для дифференциального |
уравнения |
|||||||||||||||||||||||||
y/ = cos(x − y)+ |
1.25y |
|
|
при y(0)= 0. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
1.5 + x |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
Введем начало программы |
|
|
|
|
|
|
|
y |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
ORIGIN := 1 |
f (x, y):= cos(x − y)+1.25 |
|
|
a := 0 |
b := 1 |
|
h := 0.1 |
y0 := 0 |
|||||||||||||||||||||||||
|
|
1.5 + x |
|||||||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
i := 1...11 |
xi := a |
+ h (i −1) |
|
|
|
|
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
i := 2...11 |
|
yi |
:= yi−1 + h f (xi−1 , yi−1 ) |
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
xT |
|
|
0.0 |
|
0.1 |
|
|
|
0.2 |
|
|
0.3 |
|
0.4 |
|
0.5 |
|
|
0.6 |
|
|
0.7 |
|
0.8 |
|
|
0.9 |
|
|
1.0 |
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||||||
|
yT |
|
|
0.000 |
|
0.100 |
|
|
|
0.208 |
|
|
0.323 |
|
0.445 |
|
0.575 |
|
0.710 |
|
|
0.852 |
|
0.999 |
|
|
1.152 |
|
|
1.308 |
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Аналогичного результата можно достигнуть, введя подпрограмму, реализующую формулы (7.4.1) метода Эйлера:
z1:= Eiler(y0,0,1,0.1, f )
173
Минимальная переделка подпрограммы позволяет запрограммировать исправленный и модифицированный метод Эйлера по формулам (7.4.4) и (7.4.6).
z2 := Ieiler(y0,0,1,0.1, f ) z3 := Meiler(y0,0,1,0.1, f )
Обратимся теперь к средствам пакета Mathcad. Для решения обыкновенных «неособенных» дифференциальных уравнений здесь используются две функции rkfixed и Rkadapt,
реализующие метод Рунге – Кутты четвертого порядка с постоянным и переменным шагом. |
|
Набор параметров у этих подпрограмм одинаков: |
rkfixed(y, a,b, n, f ), где y = y(x0 ), |
a и b - начало и конец интервала интегрирования, n - |
число точек и, следовательно, шаг |
h = b −n a , f (x, y)- правая часть дифференциального уравнения. Несмотря на то что при ре-
шении дифференциального уравнения функция Rkadapt использует переменный шаг, она
тем не менее представляет ответ для n точек, находящихся на одинаковом расстоянии друг |
||||||||||||
от друга, равном h = |
b − a |
. |
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|||
|
n |
|
|
|
|
|
|
|
|
|
||
Вводим следующую часть программы: |
|
yy1 |
|
|||||||||
yy1 := 0 fun(x, yy) |
:= cos(x − yy1 )+1.25 |
|
||||||||||
1.5 |
+ x |
|||||||||||
z4 := rkfixed(yy,0,1,10,fun) |
|
|
|
|||||||||
z5 := rkfixed(yy,0,1,20,fun) |
||||||||||||
z4 := Rkadapt(yy,0,1,10,fun) |
|
|
|
|
|
|
||||||
i :=1...10 Er1 := |
|
(z4 2 ) −(z5 2 ) |
2 i |
|
|
|
|
|
||||
|
|
|
|
|
|
|||||||
|
i |
|
|
i |
|
|
|
|
|
|
||
|
|
|
|
|
|
174 |
|
|
errRKF := |
max(Er1) |
errRKF = 5.291 10−3 |
|
15 |
|
||
|
|
|
Обращение к rkfixed с h1 = 0.1 и h2 = 0.05 сделано для оценки погрешности интег-
рирования по правилу Рунге (7.5.2). errRKF - оценка погрешности.
На данном учебном примере, конечно, невозможно оценить выгоды использования функции Rkadapt вместо rkfixed. В более сложных случаях Rkadapt решает уравнение более точно и быстро. Даже в этом примере точность решения по Rkadapt выше, чем по rkfixed. Это видно из следующей таблицы:
x |
|
0.1 |
|
0.2 |
|
0.3 |
|
0.4 |
|
0.5 |
z4 |
|
0.1040989 |
|
0.2161356 |
|
0.3357322 |
|
0.4625076 |
|
0.5960572 |
|
|
|
|
|
|
|
|
|
|
|
z5 |
|
0.1040990 |
0.2161359 |
0.3357326 |
0.4625081 |
0.5960578 |
||||
|
|
|
|
|
|
|
|
|
|
|
z44 |
0.1040990 |
0.2161359 |
0.3357326 |
0.4625081 |
0.5960571 |
|||||
x |
|
0.6 |
|
0.7 |
|
0.8 |
|
0.9 |
|
1.0 |
z4 |
|
0.7359363 |
|
0.8816484 |
|
1.0326377 |
|
1.1882891 |
|
1.3479326 |
|
|
|
|
|
|
|
|
|
|
|
z5 |
|
0.7359370 |
0.8816491 |
1.0326386 |
1.1882900 |
1.3479335 |
||||
|
|
|
|
|
|
|
|
|
|
|
z44 |
0.7359371 |
0.8816492 |
1.0326386 |
1.1882901 |
1.3479336 |
В заключение приведем подпрограмму RGK, реализующую формулы (7.4.13) метода Рунге – Кутты четвертого порядка. Далее, как и в предыдущем случае, произведена оценка погрешности по правилу Рунге и даны графики полученных решений. Программа метода Рунге – Кутты четвертого порядка (RGK) приведена после графиков проинтегрированных функций.
z6 := RGK (y0,0,1,0.1, f ) z7 := RGK (y0,0,1,0.05, f ) i := 1...10 Er2i := (z6 2 )i − (z7 2 )2 i
errRGK := max(Er2)
15
Задание № 1. Решить заданное дифференциальное уравнение первого порядка методом Эйлера и Рунге – Кутты четвертого порядка на отрезке x [0,1] с шагом h = 0.1 и оце-
нить погрешность интегрирования по правилу Рунге. |
|||||||||||||||||
1. |
y/ =1 + 0.2 y sin x − y2 , y |
(0) |
= 0.1. |
||||||||||||||
2. |
y/ = cos(x + y)+ 0.5(x − y), |
y(0)= 0. |
|||||||||||||||
3. |
y/ = |
cos x |
|
− 0.5y2 , |
|
y(0)= 0.2. |
|||||||||||
|
|
|
|||||||||||||||
|
|
|
x +1 |
|
|
|
|
|
|
||||||||
4. |
y/ = (1 − y2 )cos x + 0.6 y, |
y(0)= 0. |
|||||||||||||||
5. |
y/ =1+0.4 ysin x −1.5y2 , |
y(0)=1. |
|||||||||||||||
6. |
y/ = |
cos y |
+ 0.3y2 , |
|
y(0)= 0. |
|
|||||||||||
|
|
|
|
||||||||||||||
|
|
|
x + 2 |
|
|
|
|
|
|
||||||||
7. |
y/ = cos(1.5x + y)+ |
(x − y), |
y(0)= 0.5. |
||||||||||||||
8. |
y/ = 1 − sin(x + y)+ |
|
0.5y |
|
, y(0)= 0.3. |
||||||||||||
|
x + 2 |
||||||||||||||||
|
|
|
|
cos y |
|
|
|
|
|
|
|||||||
9. |
y/ = |
|
+ 0.1y2 |
, y(0)= 0. |
|||||||||||||
|
|
|
|
||||||||||||||
|
1.5 + x |
|
|
|
|
|
|
||||||||||
10. |
y/ = 0.6 sin x −1.25y2 +1, y(0)=1. |
||||||||||||||||
11. |
y/ = cos(2x + y)+1.5(x − y), y(0)= 0.1. |
||||||||||||||||
12. |
y/ =1 − |
0.1y |
− sin(2x + y), |
y(0)= 0.5. |
|||||||||||||
x + 2 |
|||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
||||||
13. |
y/ = |
|
cos y |
|
− 0.1y2 , y(0) |
= 0.2. |
|||||||||||
1.25 + x |
|||||||||||||||||
|
|
|
|
|
|
|
|
y(0)= 0.5. |
|||||||||
14. |
y/ =1 + 0.8y sin x − 2 y2 , |
||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
176 |