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

Дифференциальным уравнением первого порядка называется уравнение вида F(x,y,у’)=0 или у’=f(x,y). Функция y(x), при подстановке которой уравнение обращается в тождество, называется решением дифференциального уравнения.

Рассмотрим несколько численных методов решения дифференциальных уравнений первого порядка. Описание численных методов приводится для уравнения в виде у’=f(x,y).

Рассмотрим два варианта вывода расчетных формул

Видео:18+ Математика без Ху!ни. Дифференциальные уравнения.Скачать

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

Введение:

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

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

Возникает необходимость использовать численные методы, наиболее известным из которых является метод Рунге — Кутты [1]. Что касается Python, то в публикациях по численным методам, например [2,3], данных по применение Рунге — Кутты крайне мало, а по его модификации — методу Рунге-Кутта-Фельберга вообще нет.

В настоящее время, благодаря простому интерфейсу, наибольшее распространение в Python имеет функцию odeint из модуля scipy.integrate. Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие.

Целью настоящей публикации является сравнительный анализ перечисленных средств численного решения систем дифференциальных уравнений с модифицированным автором под Python методом Рунге-Кутта-Фельберга. В публикации так же приведены решения по краевым задачам для систем дифференциальных уравнений (СДУ).

Краткие теоретические и фактические данные по рассматриваемым методам и программным средствам для численного решения СДУ

Для одного дифференциального уравнения n – го порядка, задача Коши состоит в нахождении функции, удовлетворяющей равенству:

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

Перед решением эта задача должна быть переписана в виде следующей СДУ

(1)

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

Модуль имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (задача Коши). Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат odeint(func, y0, t[,args=(), . ]) Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2. yn], а второй – имя независимой переменной.

Функция func должна возвращать список из n значений функций при заданном значении независимого аргумента t. Фактически функция func(y,t) реализует вычисление правых частей системы (1).

Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений при t=t0.

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

Функция odeint() возвращает массив размера len(t) x len(y0). Функция odeint() имеет много опций, управляющих ее работой. Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений ei для каждого значения yi по формуле

Они могут быть векторами или скалярами. По умолчанию

Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Она создает объект ОДУ (тип scipy.integrate._ode.ode). Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы. Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей.

Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций. Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши.

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

При численном решении задачи Коши

(2)

(3)

по известному решению в точке t =0 необходимо найти из уравнения (3) решение при других t. При численном решении задачи (2),(3) будем использовать равномерную, для простоты, сетку по переменной t с шагом т > 0.

Приближенное решение задачи (2), (3) в точке обозначим . Метод сходится в точке если при . Метод имеет р-й порядок точности, если , р > 0 при . Простейшая разностная схема для приближенного решения задачи (2),(3) есть

(4)

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

Явные схемы второго и более высокого порядка аппроксимации удобно строить, ориентируясь на метод предиктор-корректор. На этапе предиктора (предсказания) используется явная схема

(5)

а на этапе корректора (уточнения) — схема

В одношаговых методах Рунге—Кутта идеи предиктора-корректора реализуются наиболее полно. Этот метод записывается в общем виде:

(6),

Формула (6) основана на s вычислениях функции f и называется s-стадийной. Если при имеем явный метод Рунге—Кутта. Если при j>1 и то определяется неявно из уравнения:

(7)

О таком методе Рунге—Кутта говорят как о диагонально-неявном. Параметры определяют вариант метода Рунге—Кутта. Используется следующее представление метода (таблица Бутчера)

Одним из наиболее распространенных является явный метод Рунге—Кутта четвертого порядка

(8)

Метод Рунге—Кутта— Фельберга

Привожу значение расчётных коэффициентов метода

(9)

С учётом(9) общее решение имеет вид:

(10)

Это решение обеспечивает пятый порядок точности, остаётся его адаптировать к Python.

Вычислительный эксперимент по определению абсолютной погрешности численного решения нелинейного дифференциального уравнения с использованием обеих функций def odein(),def oden() модуля scipy.integrate и адаптированного к Python методов Рунге—Кутта и Рунге—Кутта— Фельберга

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

Численный эксперимент по сравнению быстродействия численного решения СДУ при использовании функции ode с атрибутом dopri5 (метод Рунге – Кутты 5 порядка) и с использованием адаптированного к Python метода Рунге—Кутта— Фельберга

Сравнительный анализ проведём на примере модельной задачи, приведенной в [2]. Чтобы не повторять источник, приведу постановку и решение модельной задачи из [2].

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

где – радиус вектор движущегося тела, – вектор скорости тела, – коэффициент сопротивления, вектор силы веса тела массы m, g – ускорение свободного падения.

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

К системе следует добавить начальные условия: (h начальная высота), . Положим . Тогда соответствующая система ОДУ 1 – го порядка примет вид:

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

Flight time = 1.2316 Distance = 5.9829 Height =1.8542
Flight time = 1.1016 Distance = 4.3830 Height =1.5088
Flight time = 1.0197 Distance = 3.5265 Height =1.2912
Flight time = 0.9068 Distance = 2.5842 Height =1.0240
Время на модельную задачу: 0.454787

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

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6

Функция increment(f, t, y, tau) обеспечивает пятый порядок численного метода решения. Остальные особенности программы можно посмотреть в следующем листинге:

Время на модельную задачу: 0.259927

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

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

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

(11)

Для решения задачи (11) используем следующий алгоритм:

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

Введем обозначение для решения задачи Коши:

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

Введем обозначение для решения задачи Коши:

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

Введем обозначение для решения задачи Коши:

4. Общее решение краевой задачи (11) при помощи решений задач Коши записывается в виде линейной комбинации решений:

где p2, p3 — некоторые неизвестные параметры.

5. Для определения параметров p2, p3, используем краевые условия последних двух уравнений (11), то есть условия при x = b. Подставляя, получим систему линейных уравнений относительно неизвестных p2, p3:
(12)
Решая (12), получим соотношения для p2, p3.

По приведенному алгоритму с применением метода Рунге—Кутта—Фельберга получим следующую программу:

y0[0]= 0.0
y1[0]= 1.0
y2[0]= 0.7156448588231397
y3[0]= 1.324566562303714
y0[N-1]= 0.9900000000000007
y1[N-1]= 0.1747719838716767
y2[N-1]= 0.8
y3[N-1]= 0.5000000000000001
Время на модельную задачу: 0.070878

Вывод

Разработанная мною программа отличается от приведенной в [3] меньшей погрешностью, что подтверждает приведенный в начале статьи сравнительный анализ функции odeint с реализованным на Python метода Рунге—Кутта—Фельберга.

3. Н.М. Полякова, Е.В. Ширяева Python 3. Создание графического интерфейса пользователя (на примере решения методом пристрелки краевой задачи для линейных обыкновенных дифференциальных уравнений). Ростов-на-Дону 2017.

Видео:6.1 Численные методы решения задачи Коши для ОДУСкачать

Лекции по дисциплине «Численные методы» (стр. 1 )

Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6

Лекции по дисциплине

Тема 1. Теория погрешностей

§1. Источники и классификация погрешностей

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

Типы ошибок приближенных вычислений:

1. вычислительные ошибки;

2. ошибки округления;

3. проблемы устойчивости вычислительной схемы.

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

Определение. Численные методы это методы решения задач, сводящиеся к арифметическим и некоторым логическим действиям над ними, т. е. к действиям которые выполняет ЭВМ.

При численном решении математических и прикладных задач возможно появление на том или ином этапе погрешностей следующих типов:

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

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

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

Погрешности, соответствующие этим причинам называют:

ü Неустранимая погрешность;

ü Устранимая погрешность;

ü Вычислительная погрешность.

§2. Погрешность численного решения задачи

Введем некоторые переменные.

x* — точное значение вычисляемого параметра;

— значение этого параметра соответствующее принятому математическому описанию;

— решение задачи, получаемое при реализации численного метода в предположении отсутствия округлений;

— приближенное решение задачи, получаемое при реальных вычислениях;

— неустранимая погрешность;

— погрешность метода;

— вычислительная погрешность;

— полная погрешность;

Чтобы численный метод считался удачно выбранным, необходимо выполнение некоторых условий:

а) его погрешность должна быть в несколько раз меньше неустранимой погрешности (r0 0,001 Þ цифра неверная;

б) a = 27,017(6)8; «6»: единица ее разряда 0,0001; 0,003 > 0,0001 Þ цифра неверная;

в) a = 27,0(1)768; «1»: единица ее разряда 0,01; 0,003 0,00463.

Следовательно, приближенное значение числа А в пункте b) является более точным.

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

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

Тема 2. Решение нелинейного уравнения

§ 1. Общая постановка задачи

Пусть дано уравнение f(x) = 0, где f(x) — непрерывная функция.

Требуется вычислить действительные корни этого уравнения с заданной точностью.

I. Отделение корней уравнения.

Отделить корни — значит указать отрезки [ai, bi], на каждом из которых содержится ровно один корень уравнения.

а) Обозначим y = f(x) и построим график этой функции, корнями уравнения являются точки пересечения графика функции с осью (ох).

Введем обозначения y = g(x), y = h(x) и построим эти графики в одной системе координат.

Абсциссы точек пересечения и являются корнями уравнения f(x) = 0.

На выбранном отрезке [a, b] находится один корень уравнения f(x) = 0.

§2. Методы решения уравнений с одной переменной

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

а) Метод половинного деления (метод вилки)

Вилка это любой отрезок, на концах которого функция имеет различные знаки.

Пусть функция y = f(x) определена и непрерывна при всех x на отрезке [a,b] и имеет на концах этого отрезка значения разных знаков, т. е. f(а)×f(b) 0) Þ данный отрезок разбивается на два отрезка [a0, с0] и [с0, b0].

Найдем значение функции в этой точке f(с0). Если знак функции в т. С совпадает со знаком функции в т. b0, то в дальнейшем вместо f(b0) будем использовать f(с0). (Соответственно, если знак функции в т. С совпадает со знаком функции в т. а0, то вместо f(а0) будем использовать f(с0))

Таким образом из этих двух отрезков выбирают тот, на концах которого функция имеет значения разных знаков. Получается новый отрезок — [a1, b1].

Т. е. f(а1) × f(b1) 0 (в частности за x0 можно принять тот из концов отрезка [a, b], в котором соблюдено это условие).

Проведем в точке М0(x0, f(x0)) касательную к кривой f(x). За приближенное значение корня берется абсцисса точки пересечения этой касательной с осью (ох), т. е. точка М1(х1, 0).

Рассмотрим случай, когда f ’(x) и f ’’(x) имеют одинаковые знаки.

Пусть у = f(x). Напишем уравнение касательной в точке М0(x0, f(x0))

Подставим точку М1(x1, 0) в уравнение касательной:

0 = f(x0) + f ’(x0)(x1 – x0)

f ’(x0)(x1 – x0) = – f(x0)

Проведем теперь касательную в точке М1(x1, f(x1))

y = f(x1) + f ’(x1)(x – x1)

и подставим точку М2(х2, 0).

0 = f(x1) + f ’(x1)(x2 – x1)

f ’(x1)(x2 – x1) = – f(x1)

Полученная таким образом последовательность x0, x1, x2, … имеет своим пределом искомый корень.

Замечание: В случае, когда f ’(x) и f ’’(x) имеют разные знаки, формулы для нахождения xn+1 имеют тот же вид. (Доказать самостоятельно)

За исходную точку х0 следует выбирать тот конец отрезка [a, b], в котором знак функции совпадает со знаком f ’’(x).

Оценка погрешности (имеет место не только для метода касательных).

y = f(x) – дифференцируема на [a, b]

x* – точное значение корня уравнения f(x) = 0.

x* Î [a, b], сÎ [a, b], с ¹ x*

Рассмотрим отрезок с концами x* и с и к этому отрезку применим теорему Лагранжа:

Существует такая точка l из отрезка с концами x* и с для которой выполняется равенство

f(с) – f(x*) = f ’(l)(с – x*),

так как x* – точное значение корня уравнения f(x) = 0, то f(x*) = 0 Þ f(с) = f ’(l)(с – x*).

Т. к. точка с не является корнем этого уравнения, то можем записать f(с) ¹ 0 Þ f ’(l) ¹ 0. Выразим (с – x*):

Отсюда можем написать ,

где Þ Þ Þ

Þ – условие окончания вычислений.

Метод заключается в том, что дуга графика функции f(x) = 0 на отрезке [a,b] заменяется стягивающей ее хордой. В качестве приближенного значения корня берется точка пересечения хорды с осью (ох).

🎬 Видео

Численное решение задачи Коши методом ЭйлераСкачать

Численные методы. Лекция 9: численные методы решения дифференциальных уравненийСкачать

5 Численное решение дифференциальных уравнений Part 1Скачать

ЧМ-7. Численное интегрирование ОДУ. Часть 1/2Скачать

МЗЭ 2022 Численное решение дифференциальных уравнений. Неявный метод Эйлера. Ложкин С.А.Скачать

Обработка результатов эксперимента. 1. Классификация погрешностейСкачать

Метод ЭйлераСкачать

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

МЗЭ 2022 Численное решение дифференциальных уравнений Метод Эйлера Ложкин С. А.Скачать

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

2.2 Итерационные методы решения СЛАУ (Якоби, Зейделя, релаксации)Скачать

Метод Ньютона (метод касательных) Пример РешенияСкачать

«Фин. грамотность — от классиков до современности: советы литературных гениев и цифровые решения»Скачать

Численные методы. Вебинар. ( 03.09.21 ) Абсолютная и относительная погрешности.Скачать

Линейное неоднородное дифференциальное уравнение второго порядка с постоянными коэффициентамиСкачать

Численное решение дифференциальных уравнений (задачи Коши)Скачать