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

Метод Рунге — Кутты

Этот онлайн калькулятор реализует классический метод Рунге — Кутты (встречается также название метод Рунге — Кутта) четвертого порядка точности. Метод используется для решения дифференциальных уравнений первой степени с заданным начальным значением

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

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

и ввести правую часть уравнения f(x,y) в поле y’ калькулятора.

Также вам понадобится ввести начальное значение

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

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

Описание метода можно найти под калькулятором.

Методы численного дифференцирования и интегрирования

Общие понятия

Если функция [math]f(x)[/math] непрерывна на отрезке [math][a,b][/math] и известна ее первообразная [math]F(x)[/math] , то определенный интеграл от этой функции может быть вычислен по формуле Ньютона-Лейбница

где [math]F'(x)=f(x)[/math] . Однако во многих случаях возникают большие трудности, связанные с нахождением первообразной, или эта задача не может быть решена элементарными способами. Например, в элементарных функциях не выражается интеграл [math]\textstyle<\int\limits_<1>^ <2>\dfrac<\ln x>>[/math] .

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

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

Традиционно интегралы вычисляются с помощью квадратурных формул (явно), выражающихся линейными комбинациями значений функции [math]y_i=f(x_i)[/math] в узлах сетки [math]\Omega_n[/math] , функционально-дифференциальной квадратурной формулы Эйлера-Маклорена, формул Гаусса-Кристоффеля, Маркова и различных нестандартных формул. В данной главе кратко описаны как классические, так и новые методы численного дифференцирования и интегрирования, разработанные на основе аппарата интегрально-дифференциальных сплайнов. Приведен ряд обобщенных явных квадратурных формул функционально-дифференциального типа, а также описан новый способ вычисления интегралов на основе решения систем линейных алгебраических уравнений. Последний подход в отличие от способа вычисления интегралов по квадратурным формулам носит неявный характер и связан с построением неявных алгоритмов.

Особенность аппарата численного дифференцирования и интегрирования, изложенного в данной главе, состоит в том, что в некоторых аппроксимацион-ных операторах дифференцирования могут использоваться как значения сеточной функции в узлах, так и значения определенных интегралов на частичных отрезках [math][x_i,x_][/math] , а в аппроксимационных операторах интегрирования — значения функции и ее производных в узлах. В соответствии с этим далее выделяются классические постановки задач численного дифференцирования и интегрирования, формулируемые обычно на равномерной сетке [math](h_=\text)[/math] , и обобщенные постановки, учитывающие связи производных с интегралами в задаче численного дифференцирования и интегралов с производными в задаче численного интегрирования, формулируемые на неравномерной сетке. Классические постановки ниже нумеруются цифрой 1, а обобщенные — цифрой 2. Данные постановки относятся только к явным аппроксимационным формулам численного дифференцирования и интегрирования. Соответствующие неявные алгоритмы рассматриваются отдельно.

Постановка задачи численного дифференцирования

(h_= h=\text)[/math] заданы:
а) сеточная функция [math]y_i= f(x_i),

i=\overline<0,n>[/math] своими значениями [math]f_i=f(x_i)[/math] ;
б) точки [math]x_j[/math] сетки [math]\Omega_n[/math] , в которых требуется найти значения производных;
в) желаемый порядок [math]t[/math] точности (аппроксимации) относительно [math]h[/math] .

Требуется с заданным порядком точности (аппроксимации) вычислить значения производных [math]\Bigl.<\widehat^<(p)>(x)>\Bigr|_[/math] в точках [math]x_j[/math] сетки, где [math]p[/math] — порядок производной.

Иначе требуется получить аппроксимационный оператор [math]\widehat^<(p)>(x_j)[/math] , удовлетворяющий условию [math]\bigl|f^<(p)>(x_j)-\widehat^<(p)>(x_j)\bigr|\leqslant C\,h^t[/math] , где [math]C=\text[/math] , не зависящая от величины шага [math]h[/math] .

Задача 2. Пусть на отрезке [math][a,b][/math] в общем случае на неравномерной сетке

а) сеточная функция [math]y_i= f(x_i),

i=\overline<0,n>[/math] , своими значениями [math]f_i= f(x_i)[/math] и (или) значениями интегралов [math]\textstyle^= \int\limits_^> f(x)dx>[/math] на частичных отрезках [math][x_i,x_][/math] ;

б) точки [math]x_j[/math] сетки [math]\Omega_n[/math] , в которых требуется найти значения производных;

в) желаемый порядок [math]t[/math] точности (аппроксимации) относительно величины шага.

Требуется с заданным порядком точности (аппроксимации) вычислить значения производных [math]\Bigl.<\widehat^<(p)>(x)>\Bigr|_[/math] в точках [math]x_j[/math] сетки (получить аппроксимационный оператор), где [math]p[/math] — порядок производной.

Постановка задачи численного интегрирования

Задача 1. Пусть на отрезке [math][a,b][/math] на равномерной сетке [math]\Omega_n

а) сеточная функция [math]y_i= f(x_i),

i=\overline<0,n>[/math] , своими значениями [math]f_i= f(x_i)[/math] или сеточное представление формульной функции [math]y=f(x)[/math] ;

б) желаемый порядок [math]t[/math] точности (аппроксимации) относительно величины шага [math]h[/math] .

Требуется с заданным порядком точности вычислить значение интеграла

Иначе требуется получить аппроксимационный оператор интегрирования [math]\widehat_a^b[/math] , удовлетворяющий условию [math]\bigl|\widehat_a^b-I_a^b\bigr| \leqslant C\,h^t[/math] , где [math]C=\text[/math] , не зависящая от [math]h[/math] .

Задача 2. Пусть на отрезке [math][a,b][/math] в общем случае на неравномерной сетке

а) сеточная функция [math]y_i=f(x_i)

i=\overline<0,n>[/math] , своими значениями [math]f_i= f(x_i)[/math] и возможно значениями производных [math]f^<(p)>(x)[/math] ;

б) желаемый порядок [math]t[/math] точности (аппроксимации) относительно шага.

Требуется с заданным порядком точности вычислить значение интеграла (получить аппроксимационный оператор)

Отметим, что символом » ^ » здесь и далее обозначаются операторы дифференцирования и интегрирования.

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

1) разложения функции [math]f(x)[/math] по формуле Тейлора;

2) разложения первообразной функции [math]F(x)[/math] по формуле Тейлора;

3) замены заданной функции интерполяционными многочленами и последующим дифференцированием или интегрированием этих многочленов. Дополнительные возможности построения аппроксимационных операторов предоставляют интегрально-дифференциальные сплайны (ИД-сплайны).

Общая интегрально-функционально-дифференциальная формула для аппрок-симационного оператора дифференцирования [math]\widehat^<(p)>(x_j)[/math] может быть записана в форме

где [math]a_(h),\, b_(h),\, c_(h),\,\ldots[/math] — коэффициенты, индексы при [math]h[/math] для упрощения записи не указаны.

Общая функционально-дифференциальная формула для аппроксимационного оператора интегрирования [math]\widehat_^[/math] может быть записана в виде

где [math]a_(h),\, b_(h),\, c_(h),\,\ldots[/math] — коэффициенты.

Возможность такого комбинированного представления следует из разложения первообразной [math]F(x)[/math] по формуле Тейлора относительно точки [math]x_i[/math] и последующего выражения из него производной некоторого порядка или определенного интеграла.

В формулы (5.1), (5.2), очевидно, входят суммы нескольких групп линейных комбинаций. Так, вторая сумма в (5.1) и первая сумма в (5.2) соответствуют функциональным комбинациям, первая сумма в (5.1) соответствует интегральной части суммы. Дифференциальным комбинациям соответствуют третья и последующие части в (5.1) и вторая, третья и т.д. суммы в (5.2).

Если в (5.1) присутствуют только интегральные комбинации, формула численного дифференцирования называется интегральной, если только функциональные комбинации — функциональной (или точечной), а если только дифференциальные комбинации — дифференциальной. Если же в (5.1) отсутствуют интегральные комбинации, формула называется функционально-дифференциальной, если отсутствуют функциональные комбинации — интегрально-дифференциальной, если отсутствуют дифференциальные комбинации — интегрально-функциональной.

Аналогичная классификация справедлива и для формулы (5.2) численного интегрирования.

Метод Рунге уточнения результатов численного дифференцирования и интегрирования

Этот метод является апостериорным, и он основан на получении поправки, выражаемой через два приближенных результата [math]z_[/math] и [math]z_h[/math] , получающихся на двух равномерных сетках с шагами [math]h[/math] и [math]kh[/math] , где [math]k[/math] — коэффициент, характеризующий уменьшение шага [math](k или его увеличение 1)»>[math](k>1)[/math] при переходе с одной сетки на другую. Наиболее распространен способ генерации поправки в итерационном процессе с [math]k=1\!\!\not<\phantom<|>>\,2[/math] при необходимости снижения погрешности и [math]k=2[/math] при возможном ее увеличении.

Общая формула уточненного результата, получающегося на двух сетках с [math]h[/math] и [math]kh[/math] , записывается в виде

где [math]p[/math] — порядок точности метода, по которому рассчитаны [math]z_[/math] и [math]z_h[/math] . Второе слагаемое в правой части (5.79) представляет собой поправку или апостериорную оценку точности результата, получаемого на сетке с параметром [math]kh[/math] . Если заданная точность не достигается, то при [math]k процесс продолжается путем дробления шага и перехода на новую, более мелкую сетку. Если же указанное.слагаемое намного меньше заданной точности, то в алгоритме необходимо перейти на укрупненную сетку, полагая, например, [math]k=2[/math] .

Данную методику рассмотрим применительно к расчету производных на основе сгущающихся сеток [math](k=1\!\!\not<\phantom<|>>\,2)[/math] , когда для первоначального вычисления производных принимается одна из формул (5.10). Схема уточнения приведена в табл. 5.5, где численные значения производных обозначены символом [math]Z[/math] с двумя нижними индексами, первый из которых (т>\) указывает на ступени уточнения ( [math]m=1[/math] — первая, [math]m=2[/math] — вторая и т.д.), а второй — [math]k[/math] — на порядковые номера значений производных, полученных при различных шагах на каждой ступени уточнения. Начиная с третьей колонки таблицы помещаются значения [math]Z_<0,0>,Z_<0,1>,\ldots[/math] , полученные по некоторой исходной аппроксимационной формуле при различных шагах разбиения сетки:

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

где значение [math]2m[/math] соответствует порядку аппроксимации [math]p[/math] .

Порядок уточнения численных результатов

Пусть по заданной сеточной функции [math]y_i= f(x_i),

i=\overline<0,n>[/math] , в некоторой фиксированной точке [math]x_j=x_i[/math] производная вначале рассчитывается по формуле (5.10) второго порядка аппроксимации

Тогда уточнение до порядка [math]O(h^6)[/math] производится следующим образом.

1. По формуле (5.81) рассчитываются значения [math]Z_<0,0>,\, Z_<0,1>,\, Z_<0,2>[/math] производных в заданной точке с использованием таблицы значений исходной сеточной функции. При этом для вычисления [math]Z_<0,k>,

k=0,1,2[/math] , берутся значения [math]f(x_i)[/math] в тех точках сетки, которые соответствуют шагам [math]h_0,\, h_1[/math] и [math]h_2[/math] . Все величины, входящие в третью колонку таблицы, имеют порядок аппроксимации [math]O(h^2)[/math] , что указано в нижней строке табл. 5.5.

2. Реализуется первая ступень процесса уточнения по формуле (5.80), в которой следует принять [math]m=1[/math] (знаменатель второго слагаемого формулы (5.80) равен 3), а [math]p=2m=2[/math] . При этом вычисляются [math]Z_<1,0>[/math] и [math]Z_<1,1>[/math] . Порядок аппроксимации результатов, получаемых в первой ступени уточнения, равен четырем.

3. Реализуется вторая ступень процесса уточнения по формуле (5.80), в которой следует принять [math]m=2[/math] (знаменатель второго слагаемого формулы (5.80) равен 15), а [math]p=2m=4[/math] . При этом вычисляется окончательное значение [math]Z_<2,0>[/math] , порядок аппроксимации которого равен шести. Процесс уточнения до заданного порядка [math]p=6[/math] закончен.

Замечание. В случае, когда необходимо обеспечить порядок аппроксимации, превышающий полученный, например восьмой, производятся дополнительные расчеты значений [math]Z_<0,k>[/math] при более мелких шагах и расчеты [math]Z_<1,2>,\, Z_<2,1>[/math] по первой и второй ступени уточнения. После этого рассчитывается [math]Z_<3,0>[/math] с порядком [math]O(h^8)[/math] (см. табл.5.5). Данный процесс может быть продолжен и далее.

Пример 5.13. Для сеточного представления функции [math]y=\sin x,

x\in[0; 2\pi\!\!\not<\phantom<|>>\,5][/math] , заданной на сетке с шагом [math]h= \pi\!\!\not<\phantom<|>>\,10[/math] , на основе применения формулы (5.10) второго порядка и метода Рунге найти первую производную [math]\Bigl.\Bigr|_<5>>[/math] с порядком [math]O(h^4)[/math] .

1. По формуле (5.81) рассчитываются значения [math]Z_<0,0>[/math] и [math]Z_<0,1>[/math] производных в точке [math]x_2=\frac<\pi><5>[/math] . При этом [math]Z_<0,0>[/math] получается при [math]h_0= \frac<2\pi><10>[/math] , а значение [math]Z_<0,1>[/math] при [math]h_1= \frac<\pi><10>[/math] . Записывая формулу (5.81) для указанных шагов, получаем

Исходная сеточная функция и результаты расчетов производной в точке [math]x_2[/math] помещены в табл. 5.6.

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

2. Реализуется первая (в данном случае и последняя) ступень процесса уточнения по формуле (5.80), в которой [math]m=1[/math] . При этом принимается [math]f’_<2,\frac<2>>= Z_<0,1>,

Полученный результат соответствует порядку аппроксимации [math]O(h^4)[/math] .

Замечание. Для уточнения результатов численного интегрирования метод Рунге реализуется аналогично.