Численные методы решения дифференциальных уравнений. Численное решение обыкновенных дифференциальных уравнений Что значит решить дифференциальное уравнение численным методом
Лабораторная работа 1
Численные методы решения
обыкновенных дифференциальных уравнений (4 часа)
При решении многих физических и геометрических задач приходится искать неизвестную функцию по данному соотношению между неизвестной функцией, ее производными и независимыми переменными. Такое соотношение называется дифференциальным уравнением , а отыскание функции, удовлетворяющей дифференциальному уравнению, называется решением дифференциального уравнения.
Обыкновенным дифференциальным уравнением называется равенство
, (1)в котором
- независимая переменная, изменяющаяся в некотором отрезке , а - неизвестная функция y ( x ) и ее первые n производные. называется порядком уравнения .Задача заключается в нахождении функции y, удовлетворяющей равенству (1). Более того, не оговаривая это отдельно, будем предполагать, что искомое решение обладает той или иной степенью гладкости, необходимой для построения и «законного» применения того или иного метода.
Различают два типа обыкновенных дифференциальных уравнений
Уравнения без начальных условий
Уравнения с начальными условиями.
Уравнения без начальных условий - это уравнение вида (1).
Уравнение с начальными условиями - это уравнение вида (1), в котором требуется найти такую функцию
, которая при некотором удовлетворяет следующим условиям: ,т.е. в точке
функция и ее первые производных принимают наперед заданные значения.Задачи Коши
При изучении способов решения дифференциальных уравнений приближенными методами основной задачей считается задача Коши.
Рассмотрим наиболее популярный метод решения задачи Коши – метод Рунге-Кутта. Этот метод позволяет строить формулы расчета приближенного решения практически любого порядка точности.
Выведем формулы метода Рунге-Кутта второго порядка точности. Для этого решение представим куском ряда Тейлора, отбрасывая члены с порядком выше второго. Тогда приближенное значение искомой функции в точке x 1 можно записать в виде:
(2)Вторую производную y "( x 0 ) можно выразить через производную функции f ( x , y ) , однако в методе Рунге-Кутта вместо производной используют разность
соответственно подбирая значения параметров
Тогда (2) можно переписать в виде:
y 1 = y 0 + h [ β f ( x 0 , y 0 ) + α f ( x 0 + γh , y 0 + δh )], (3)
где α , β , γ и δ – некоторые параметры.
Рассматривая правую часть (3) как функцию аргумента h , разложим ее по степеням h :
y 1 = y 0 +( α + β ) h f ( x 0 , y 0 ) + αh 2 [ γ f x ( x 0 , y 0 ) + δ f y ( x 0 , y 0 )],
и выберем параметры α , β , γ и δ так, чтобы это разложение было близко к (2). Отсюда следует, что
α + β =1, αγ =0,5, α δ =0,5 f ( x 0 , y 0 ).
С помощью этих уравнений выразим β , γ и δ через параметры α , получим
y 1 = y 0 + h [(1 - α ) f ( x 0 , y 0 ) + α f ( x 0 +, y 0 + f ( x 0 , y 0 )], (4)
0 < α ≤ 1.
Теперь, если вместо (x 0 , y 0 ) в (4) подставить (x 1 , y 1 ), получим формулу для вычисления y 2 – приближенного значения искомой функции в точке x 2 .
В общем случае метод Рунге-Кутта применяется на произвольном разбиении отрезка [ x 0 , X ] на n частей, т.е. с переменным шагом
x 0 , x 1 , …,x n ; h i = x i+1 – x i , x n = X. (5)
Параметры α выбирают равными 1 или 0,5. Запишем окончательно расчетные формулы метода Рунге-Кутта второго порядка с переменным шагом для α =1:
y i+1 =y i +h i f(x i + , y i + f(x i , y i)), (6.1)
i = 0, 1,…, n -1.
и α =0,5:
y i+1 =y i + , (6.2)
i = 0, 1,…, n -1.
Наиболее употребляемые формулы метода Рунге-Кутта – формулы четвертого порядка точности:
y i+1 =y i + (k 1 + 2k 2 + 2k 3 + k 4),
k 1 =f(x i , y i), k 2 = f(x i + , y i + k 1), (7)
k 3 = f(x i + , y i + k 2), k 4 = f(x i +h, y i +hk 3).
Для метода Рунге-Кутта применимо правило Рунге для оценки погрешности. Пусть y ( x ; h ) – приближенное значение решения в точке x , полученное по формулам (6.1), (6.2) или (7) с шагом h , а p – порядок точности соответствующей формулы. Тогда погрешность R ( h ) значения y ( x ; h ) можно оценить, используя приближенное значение y ( x ; 2 h ) решения в точке x , полученное с шагом 2 h :
(8)где p =2 для формул (6.1) и (6.2) и p =4 для (7).
Численное решение дифференциальных уравнений
Многие задачи науки и техники сводятся к решению обыкновенных дифференциальных уравнений (ОДУ). ОДУ называются такие уравнения, которые содержат одну или несколько производных от искомой функции. В общем виде ОДУ можно записать следующим образом:
Где x – независимая переменная, - i-ая производная от искомой функции. n - порядок уравнения. Общее решение ОДУ n–го порядка содержит n произвольных постоянных , т.е. общее решение имеет вид .
Для выделения единственного решения необходимо задать n дополнительных условий. В зависимости от способа задания дополнительных условий существуют два различных типа задач: задача Коши и краевая задача. Если дополнительные условия задаются в одной точке, то такая задача называется задачей Коши. Дополнительные условия в задаче Коши называются начальными условиями. Если же дополнительные условия задаются в более чем одной точке, т.е. при различных значениях независимой переменной, то такая задача называется краевой. Сами дополнительные условия называются краевыми или граничными.
Ясно, что при n=1 можно говорить только о задачи Коши.
Примеры постановки задачи Коши :
Примеры краевых задач :
Решить такие задачи аналитически удается лишь для некоторых специальных типов уравнений.
Численные методы решения задачи Коши для ОДУ первого порядка
Постановка задачи . Найти решение ОДУ первого порядка
На отрезке при условии
При нахождении приближенного решения будем считать, что вычисления проводятся с расчетным шагом , расчетными узлами служат точки промежутка [x 0 , x n ].
Целью является построение таблицы
x i |
x n |
|||
y i |
y n |
т.е. ищутся приближенные значения y в узлах сетки.
Интегрируя уравнение на отрезке , получим
Вполне естественным (но не единственным) путем получения численного решения является замена в нем интеграла какой–либо квадратурной формулой численного интегрирования. Если воспользоваться простейшей формулой левых прямоугольников первого порядка
,
то получим явную формулу Эйлера :
Порядок расчетов:
Зная , находим , затем т.д.
Геометрическая интерпретация метода Эйлера :
Пользуясь тем, что в точке x 0 известно решение y (x 0) = y 0 и значение его производной , можно записать уравнение касательной к графику искомой функции в точке :. При достаточно малом шаге h ордината этой касательной, полученная подстановкой в правую часть значения , должна мало отличаться от ординаты y (x 1) решенияy (x ) задачи Коши. Следовательно, точка пересечения касательной с прямой x = x 1 может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую , которая приближенно отражает поведение касательной к в точке . Подставляя сюда (т.е. пересечение с прямой x = x 2), получим приближенное значение y (x ) в точке x 2: и т.д. В итоге для i –й точки получим формулу Эйлера.
Явный метод Эйлера имеет первый порядок точности или аппроксимации.
Если использовать формулу правых прямоугольников: , то придем к методу
Этот метод называют неявным методом Эйлера , поскольку для вычисления неизвестного значения по известному значению требуется решать уравнение, в общем случае нелинейное.
Неявный метод Эйлера имеет первый порядок точности или аппроксимации.
В данном методе вычисление состоит из двух этапов:
Данная схема называется еще методом предиктор – корректор (предсказывающее – исправляющее). На первом этапе приближенное значение предсказывается с невысокой точностью (h), а на втором этапе это предсказание исправляется, так что результирующее значение имеет второй порядок точности.
Методы Рунге – Кутта: идея построения явных методов Рунге–Кутты p –го порядка заключается в получении приближений к значениям y (x i +1) по формуле вида
…………………………………………….
Здесь a n , b nj , p n , – некоторые фиксированные числа (параметры).
При построения методов Рунге–Кутты параметры функции (a n , b nj , p n ) подбирают таким образом, чтобы получить нужный порядок аппроксимации.
Схема Рунге – Кутта четвертого порядка точности :
Пример . Решить задачу Коши:
Рассмотреть три метода: явный метод Эйлера, модифицированный метод Эйлера, метод Рунге – Кутта.
Точное решение:
Расчетные формулы по явному методу Эйлера для данного примера:
Расчетные формулы модифицированного метода Эйлера:
Расчетные формулы метода Рунге – Кутта:
y1 – метод Эйлера, y2 – модифицированный метод Эйлера, y3 – метод Рунге Кутта.
Видно, что самым точным является метод Рунге – Кутта.
Численные методы решения систем ОДУ первого порядка
Рассмотренные методы могут быть использованы также для решения систем дифференциальных уравнений первого порядка.
Покажем это для случая системы двух уравнений первого порядка:
Явный метод Эйлера:
Модифицированный метод Эйлера:
Схема Рунге – Кутта четвертого порядка точности:
К решению систем уравнений ОДУ сводятся также задачи Коши для уравнений высших порядков. Например, рассмотрим задачу Коши для уравнения второго порядка
Введем вторую неизвестную функцию . Тогда задача Коши заменяется следующей:
Т.е. в терминах предыдущей задачи: .
Пример. Найти решение задачи Коши :
На отрезке .
Точное решение:
Действительно:
Решим задачу явным методом Эйлера, модифицированным методом Эйлера и Рунге – Кутта с шагом h=0.2.
Введем функцию .
Тогда получим следующую задачу Коши для системы двух ОДУ первого порядка:
Явный метод Эйлера:
Модифицированный метод Эйлера:
Метод Рунге – Кутта:
Схема Эйлера:
Модифицированный метод Эйлера:
Схема Рунге - Кутта:
Max(y-y теор)=4*10 -5
Метод конечных разностей решения краевых задач для ОДУ
Постановка задачи : найти решение линейного дифференциального уравнения
удовлетворяющего краевым условиям:. (2)
Теорема. Пусть . Тогда существует единственное решение поставленной задачи.
К данной задаче сводится, например, задача об определении прогибов балки, которая на концах опирается шарнирно.
Основные этапы метода конечных разностей:
1) область непрерывного изменения аргумента () заменяется дискретным множеством точек, называемых узлами: .
2) Искомая функция непрерывного аргумента x, приближенно заменяется функцией дискретного аргумента на заданной сетке, т.е. . Функция называется сеточной.
3) Исходное дифференциальное уравнение заменяется разностным уравнением относительно сеточной функции. Такая замена называется разностной аппроксимацией.
Таким образом, решение дифференциального уравнения сводится к отысканию значений сеточной функции в узлах сетки, которые находятся из решения алгебраических уравнений.
Аппроксимация производных.
Для аппроксимации (замены) первой производной можно воспользоваться формулами:
- правая разностная производная,
- левая разностная производная,
Центральная разностная производная.
т.е., возможно множество способов аппроксимации производной.
Все эти определения следуют из понятия производной как предела: .
Опираясь на разностную аппроксимацию первой производной можно построить разностную аппроксимацию второй производной:
Аналогично можно получить аппроксимации производных более высокого порядка.
Определение. Погрешностью аппроксимации n- ой производной называется разность: .
Для определения порядка аппроксимации используется разложение в ряд Тейлора.
Рассмотрим правую разностную аппроксимацию первой производной:
Т.е. правая разностная производная имеет первый по h порядок аппроксимации.
Аналогично и для левой разностной производной.
Центральная разностная производная имеет второй порядок аппроксимации .
Аппроксимация второй производной по формуле (3) также имеет второй порядок аппроксимации.
Для того чтобы аппроксимировать дифференциальное уравнение необходимо в нем заменить все производные их аппроксимациями. Рассмотрим задачу (1), (2) и заменим в(1) производные:
В результате получим:
(4)
Порядок аппроксимации исходной задачи равен 2, т.к. вторая и первая производные заменены с порядком 2, а остальные – точно.
Итак, вместо дифференциальных уравнений (1), (2) получена система линейных уравнений для определения в узлах сетки.
Схему можно представить в виде:
т.е., получили систему линейных уравнений с матрицей:
Данная матрица является трехдиагональной, т.е. все элементы, которые расположены не на главной диагонали и двух прилегающих к ней диагоналях равны нулю.
Решая полученную систему уравнений, мы получим решение исходной задачи.
Дифференциальные уравнения - это уравнения, в которые неизвестная функция входит под знаком производной. Основная задача теории дифференциальных уравнений -- изучение функций, являющихся решениями таких уравнений.
Дифференциальные уравнения можно разделить на обыкновенные дифференциальные уравнения, в которых неизвестные функции являются функциями одной переменной, и на дифференциальные уравнения в частных производных, в которых неизвестные функции являются функциями двух и большего числа переменных.
Теория дифференциальных уравнений в частных производных более сложная и рассматривается в более полных или специальных курсах математики.
Изучение дифференциальных уравнений начнем с наиболее простого уравнения--уравнения первого порядка.
Уравнение вида
F(x,y,y") = 0,(1)
где х -- независимая переменная; у -- искомая функция; у" -- ее производная, называется дифференциальным уравнением первого порядка.
Если уравнение (1) можно разрешить относительно у", то оно принимает вид
и называется уравнением первого порядка, разрешенным относительно производной.
В некоторых случаях уравнение (2) удобно записать в виде f (х, у) dх - dy = 0, являющемся частным случаем более общего уравнения
P(x,y)dx+Q(x,y)dy=O,(3)
где Р(х,у) и Q(х,у) -- известные функции. Уравнение в симметричной форме (3) удобно тем, что переменные х и у в нем равноправны, т. е. каждую из них можно рассматривать как функцию другой.
Дадим два основных определения общего и частного решения уравнения.
Общим решением уравнения (2) в некоторой области G плоскости Оху называется функция у=ц(х,С), зависящая от х и произвольной постоянной С, если она является решением уравнения (2) при любом значении постоянной С, и если при любых начальных условиях y x=x0 =y 0 таких, что (x 0 ;y 0)=G, существует единственное значение постоянной С = С 0 такое, что функция у=ц(х,С 0) удовлетворяет данным начальным условиям у=ц(х 0 ,С).
Частным решением уравнения (2) в области G называется функция у=ц(х,С 0), которая получается из общего решения у=ц(х,С) при определенном значении постоянной С=С 0 .
Геометрически общее решение у=ц(х,С) представляет собой семейство интегральных кривых на плоскости Оху, зависящее от одной произвольной постоянной С, а частное решение у=ц(х,С 0) -- одну интегральную кривую этого семейства, проходящую через заданную точку (х 0 ; у 0).
Приближенное решение дифференциальных уравнений первого порядка методом Эйлера. Суть этого метода состоит в том, что искомая интегральная кривая, являющаяся графиком частного решения, приближенно заменяется ломаной. Пусть даны дифференциальное уравнение
и начальные условия y |x=x0 =y 0 .
Найдем приближенно решение уравнения на отрезке [х 0 ,b], удовлетворяющее заданным начальным условиям.
Разобьем отрезок [х 0 ,b] точками х 0 <х 1 ,<х 2 <...<х n =b на n равных частей. Пусть х 1 --х 0 =х 2 -- x 1 = ... =x n -- x n-1 = ?x. Обозначим через y i приближенные значения искомого решения в точках х i (i=1, 2, ..., n). Проведем через точки разбиения х i - прямые, параллельные оси Оу, и последовательно проделаем следующие однотипные операции.
Подставим значения х 0 и у 0 в правую часть уравнения y"=f(x,y) и вычислим угловой коэффициент y"=f(x 0 ,y 0) касательной к интегральной кривой в точке (х 0 ;у 0). Для нахождения приближенного значения у 1 искомого решения заменяем на отрезке [х 0 ,x 1 ,] интегральную кривую отрезком ее касательной в точке (х 0 ;у 0). При этом получаем
y 1 - y 0 =f(x 0 ;y 0)(x 1 - x 0),
откуда, так как х 0 , х 1 , у 0 известны, находим
y1 = y0+f(x0;y0)(x1 - x0).
Подставляя значения х 1 и y 1 , в правую часть уравнения y"=f(x,y), вычисляем угловой коэффициент y"=f(x 1 ,y 1) касательной к интегральной кривой в точке (х 1 ;y 1). Далее, заменяя на отрезке интегральную кривую отрезком касательной, находим приближенное значение решения у 2 в точке х 2:
y 2 = y 1 +f(x 1 ;y 1)(x 2 - x 1)
В этом равенстве известными являются х 1 , у 1 , х 2 , а у 2 выражается через них.
Аналогично находим
y 3 = y 2 +f(x 2 ;y 2) ?x, …, y n = y n-1 +f(x n-1 ;y n-1) ?x
Таким образом, приближенно построена искомая интегральная кривая в виде ломаной и получены приближенные значения y i искомого решения в точках х i . При этом значения у i вычисляются по формуле
y i = y i-1 +f(x i-1 ;y i-1) ?x (i=1,2, …, n).
Формула и является основной расчетной формулой метода Эйлера. Ее точность тем выше, чем меньше разность?x.
Метод Эйлера относиться к численным методам, дающим решение в виде таблицы приближенных значений искомой функции у(х). Он является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других методов.
Степень точности метода Эйлера, вообще говоря, невелика. Существуют гораздо более точные методы приближенного решения дифференциальных уравнений.
Основные вопросы, рассматриваемые на лекции:
1. Постановка задачи
2. Метод Эйлера
3. Методы Рунге-Кутта
4. Многошаговые методы
5. Решение краевой задачи для линейного дифференциального уравнения 2 порядка
6. Численное решение дифференциальных уравнений в частных производных
1. Постановка задачи
Простейшим обыкновенным дифференциальным уравнением (ОДУ) является уравнение первого порядка, разрешённое относительно производной: y " = f (x, y) (1). Основная задача, связанная с этим уравнением известна как задача Коши: найти решение уравнения (1) в виде функции y (x), удовлетворяющей начальному условию: y (x0) = y0 (2).
ДУ n-ого порядка y (n) = f (x, y, y",:, y(n-1)), для которого задача Коши состоит в нахождении решения y = y(x), удовлетворяющего начальным условиям:
y (x0) = y0 , y" (x0) = y"0 , :, y(n-1)(x0) = y(n-1)0 , где y0 , y"0 , :, y(n-1)0 - заданные числа, можно свести к системе ДУ первого порядка.
· Метод Эйлера
В основе метода Эйлера лежит идея графического построения решения ДУ, однако этот же метод даёт одновременно и численную форму искомой функции. Пусть дано уравнение (1) с начальным условием (2).
Получение таблицы значений искомой функции y (x) по методу Эйлера заключается в циклическом применении формулы: , i = 0, 1, :, n. Для геометрического построения ломаной Эйлера (см. рис.) выберем полюс A(-1,0) и на оси ординат отложим отрезок PL=f(x0 , y0) (точка P - это начало координат). Очевидно, что угловой коэффициент луча AL будет равен f(x0 , y0), поэтому чтобы получить первое звено ломаной Эйлера достаточно из точки М провести прямую MM1 параллельно лучу AL до пересечения с прямой х = х1 в некоторой точке М1(х1, у1). Приняв точку М1(х1, у1) за исходную откладываем на оси Оу отрезок PN = f (x1, y1) и через точку М1 проводим прямую М1М2 | | AN до пересечения в точке М2(х2, у2) с прямой х = х2 и т.д.
Недостатки метода: малая точность, систематическое накопление ошибок.
· Методы Рунге-Кутта
Основная идея метода: вместо использования в рабочих формулах частных производных функции f (x, y) использовать лишь саму эту функцию, но на каждом шаге вычислять её значения в нескольких точках. Для этого будем искать решение уравнения (1) в виде:
Меняя α, β, r, q, будем получать различные варианты методов Рунге-Кутта.
При q=1 получаем формулу Эйлера.
При q=2 и r1=r2=½ получаем, что α, β= 1 и, следовательно, имеем формулу: , которая называется усовершенствованный метод Эйлера-Коши.
При q=2 и r1=0, r2=1 получаем, что α, β = ½ и, следовательно, имеем формулу: - второй усовершенствованный метод Эйлера-Коши.
При q=3 и q=4 также существуют целые семейства формул Рунге-Кутта. На практике они применяются наиболее часто, т.к. не наращивают ошибок.
Рассмотрим схему решения дифференциального уравнения методом Рунге-Кутта 4 порядка точности. Расчёты при использовании этого метода ведутся по формулам:
Их удобно вносить в следующую таблицу:
x | y | y" = f (x,y) | k=h · f(x,y) | Δy |
x0 | y0 | f(x0,y0) | k1(0) | k1(0) |
x0 + ½·h | y0 + ½·k1(0) | f(x0 + ½·h, y0 + ½·k1(0)) | k2(0) | 2k2(0) |
x0 + ½·h | y0 + ½·k2(0) | f(x0 + ½·h, y0 + ½·k2(0)) | k3(0) | 2k3(0) |
x0 + h | y0 + k3(0) | f(x0 + h, y0 + k3(0)) | k4(0) | k4(0) |
Δy0 = Σ / 6 | ||||
x1 | y1 = y0 + Δy0 | f(x1,y1) | k1(1) | k1(1) |
x1 + ½·h | y1 + ½·k1(1) | f(x1 + ½·h, y1 + ½·k1(1)) | k2(1) | 2k2(1) |
x1 + ½·h | y1 + ½·k2(1) | f(x1 + ½·h, y1 + ½·k2(1)) | k3(1) | 2k3(1) |
x1 + h | y1 + k3(1) | f(x1 + h, y1 + k3(1)) | k4(1) | k4(1) |
Δy1 = Σ / 6 | ||||
x2 | y2 = y1 + Δy1 | и т.д. | до получения всех искомых | значений y |
· Многошаговые методы
Рассмотренные выше методы - это так называемые методы пошагового интегрирования дифференциального уравнения. Они характерны тем, что значение решения на следующем шаге ищется с использованием решения, полученного лишь на одном предыдущем шаге. Это так называемые одношаговые методы.
Основная идея же многошаговых методов заключается в использовании при вычислении значения решения на следующем шаге нескольких предыдущих значений решения. Также эти методы носят название m-шаговых по числу m используемых для расчётов предыдущих значений решения.
В общем случае для определения приближённого решения yi+1 m-шаговые разностные схемы записываются таким образом (m 1):
Рассмотрим конкретные формулы, реализующие простейшие явный и неявный методы Адамса.
Явный метод Адамса 2 порядка (2-шаговый явный метод Адамса)
Имеем a0 = 0, m = 2.
Таким образом, - расчётные формулы явного метода Адамса 2-ого порядка.
При i = 1 имеем неизвестное y1, которое будем находить по методу Рунге-Кутта при q = 2 илиq = 4.
При i = 2, 3, : все необходимые значения известны.
Неявный метод Адамса 1 порядка
Имеем: a0 0, m = 1.
Таким образом, - расчётные формулы неявного метода Адамса 1-ого порядка.
Основная проблема неявных схем заключается в следующем: yi+1 входит и в правую и в левую часть представленного равенства, поэтому имеем уравнение для поиска значения yi+1. Данное уравнение является нелинейным и записано в форме, подходящей для итерационного решения, поэтому будем использовать метод простой итерации для его решения:
Если шаг h выбран удачно, то итерационный процесс быстро сходится.
Данный метод также не является самостартующимся. Так для вычисления y1 надо знать y1(0). Его можно найти по методу Эйлера.