Автоматический выбор шага. Правило Рунге

При выполнении практических расчетов с применением рассмотренных квадратурных формул оценку погрешности получаемых результатов можно осуществлять априорно, используя соотношения (4.3.15) – (4.3.18), (4.3.21), (4.3.22), (4.3.25) и (4.3.26). Эти соотношения позволяют получить (для заданного шага h между узлами равномерной сетки h_x) оценку сверху для погрешности аппроксимации квадратурной формулы с точностью до максимума модуля производной того или иного порядка подынтегральной функции на отрезке интегрирования. С другой стороны, основываясь на этих соотношениях можно определить максимально возможное значение h, при котором удовлетворяются заданные требования к точности результатов расчета. Однако с практической точки зрения получить достоверную оценку максимума модуля производной, даже если известно аналитическое задание подынтегральной функции, достаточно сложно. Если аналитическое задание неизвестно (что чаще всего и имеет место), а функция представлена лишь табличными значениями, то задача получения такой оценки, как правило, становится практически неразрешимой. Таким образом, возникает необходимость применения других способов оценки погрешности численного интегрирования, один из которых основан на использовании правила Рунге.

В отличие от априорной (до того как будут произведены вычисления) оценки погрешности на основании соотношений (4.3.15) – (4.3.18), (4.3.21), (4.3.22), (4.3.25) и (4.3.26) для предлагаемого подхода имеет место апостериорная оценка, которая получается в результате анализа уже произведенных вычислений.

Пусть I_k – точное, а J_k^(h) – приближенное значение определенного интеграла от функции f(x), полученное в результате расчета с шагом h по одной из рассмотренных ранее квадратурных формул на частичных отрезках разбиения [x_k, x_{k + 1}], k = 0, 1, 2, …, n – 1 отрезка [a, b]. Тогда будет справедливо приближенное равенство

I_k – J_k^(h) » q_kh^m, k = 0, 1, 2, …, n – 1, (4.3.27)

где m – порядок точности используемой квадратурной формулы.

Напомним, что m = 2 для квадратурных формул левых (4.3.9) и правых (4.3.11) прямоугольников, m = 3 - для квадратурной формулы центральных (4.3.13) прямоугольников и квадратурной формулы (4.3.19) трапеций, m = 5 - для квадратурной формулы (4.3.23) Симпсона.

Уменьшим шаг сетки в два раза и повторим расчет по той же именно???? квадратурной формуле. Обозначим результат расчета как J_k^(h/2). Для него справедливо аналогичное (4.3.27) приближенное равенство

I_k – J_k^(h/2) » q_k(h/2)^m, k = 0, 1, 2, …, n – 1. (4.3.28)

Из соотношений (4.3.27) и (4.3.28) имеем

(4.3.29)

Полученное приближенное равенство часто называют правилом Рунге.

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

(4.3.30)

Здесь I – точное значение определенного интеграла на отрезке [a, b]; J^(h) – приближенное значение этого интеграла, определенное с помощью составной квадратурной формулы.

Для каждого частичного отрезка разбиения [x_k, x_{k + 1}], k = 0, 1, 2, …, n – 1 дважды проведем вычисления по квадратурной формуле, имеющей порядок точности равный m, – с шагом h и h/2.

Согласно правилу Рунге (4.3.29) можно утверждать, что если
в результате для всех частичных отрезков разбиения будут выполняться неравенства

, (4.3.31)

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

Если для некоторых частичных отрезков неравенство (4.3.31) не будет выполнено, то для них необходимо повторять расчет, уменьшая каждый раз шаг h сетки в два раза до тех пор, пока не будет обеспечено его выполнение.

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

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

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

Основные понятия

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

Задача Коши для системы дифференциальных уравнений первого порядка формулируется следующим образом: найти функции yk = fk(x), k = 1, 2, …, n, удовлетворяющие системе дифференциальных уравнений = gk(x, y1, y2, …, yn), k = 1, 2, …, n (5.1.1) и начальным условиям yk(x0) = y0(k), k = 1, 2, …, n, (5.1.2) где y0(k), k = 1, 2, …, n – заданные числа.

Более подробно эту задачу можно сформулировать так: найти функции y_k = f_k(x), k = 1, 2, …, n, удовлетворяющие системе дифференциальных уравнений

(5.1.3)

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

y₁(x₀) = y₍₁₎₀, y₂(x₀) = y₍₂₎₀, …, y_{n – 1}(x₀) = y_{(n – 1)0}, y_n(x₀) = y_(n)0, (5.1.4)

где y_(1)0, y₍₂₎₀, …, y_{(n – 1)0,} y_(n)0 – заданные числа.

В дальнейшем предполагается, что условия теоремы существования и единственности для задачи (5.1.3), (5.1.4) выполнены.

Положим a:= x₀ и будем искать решение задачи Коши (5.1.3), (5.1.4) на отрезке [a, b]. При построении численного решения получим приближенное решение задачи

j_(1)h(x)»f₁(x); j_(2)h(x)»f₂(x);…;j_{(n – 1)h}(x)»f_{n – 1}(x); j_(n)h(x)»f_n(x), (5.1.5)

которое представлено значениями функций j_(k)h(x), k = 1, 2, …, n в узлах заданной на отрезке [a, b] одномерной сетки

h_x = {x_i/x_i = x_{i – 1} + h_i, h_i > 0, i = 1, 2, 3, …, m; x₀ = a, x_m = b}, (5.1.6)

где x_i, i = 0, 1, 2, …, m – узлы одномерной сетки; h_i, i = 1, 2, 3, …, m – ее шаг.

Для простоты будем полагать, что сетка h_x равномерная:
h_i = h, i = 1, 2, 3, …, m. Таким образом, в результате численного решения задачи (5.1.3) - (5.1.4) получим значения j_(k)h(x₀), j_(k)h(x₁), j_(k)h(x₂), …, j_(k)h(x_n) функций j_(k)h(x), k = 1, 2, …, n. Эти результаты удобно представлять в табличной форме (табл. 5.1.1).

Таблица 5.1.1

Как известно, к задаче (5.1.3) - (5.1.4) сводится и задача Коши для обыкновенного дифференциального уравнения n-го порядка:

найти функцию y = f(x), удовлетворяющую дифференциальному уравнению

y⁽ⁿ⁾ = g(x, y, y^¢, y⁽²⁾, …, y^{(n – 2)}, y^{(n – 1)}) (5.1.7)

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

(5.1.8)

где – заданные числа (здесь y^(k) - производная k-го порядка функции y = f(x)).

Сведение задачи (5.1.7) - (5.1.8) к решению системы дифференциальных уравнений вида (5.1.3), (5.1.4) осуществляется посредством переобозначения:

y₁:= y^¢; y₂:= y⁽²⁾;…; y_{n – 1}:= y^{(n – 1)}; y_n:= y⁽ⁿ⁾. (5.1.9)

В результате вместо дифференциального уравнения n-го порядка (5.1.7) может быть записана система n дифференциальных уравнений первого порядка

(5.1.10)

а начальные условия (5.1.8) переписываются в обозначениях (5.1.9), что приводит к задаче вида (5.1.3), (5.1.4).

Как частный случай (n = 1) задачи (5.1.3), (5.1.4) запишем задачу Коши для дифференциального уравнения первого порядка, которая может быть сформулирована следующим образом:

найти функцию y = f(x), удовлетворяющую дифференциальному уравнению

y¢ = g(x, y) (5.1.11)

и начальному условию

y(x₀) = y₀, (5.1.12)

где y₀ – заданное число.

Предположим, как и ранее, что условия теоремы существования
и единственности решения задачи (5.1.11), (5.1.12) выполнены. Положим a:= x₀ и будем искать решение задачи Коши (5.1.11), (5.1.12) на отрез-
ке [a, b]. В результате численного решения поставленной задачи получим приближенное решение

j_h(x) » f(x), (5.1.13)

которое представлено значениями j_h(x) в узлах заданной на отрезке [a, b] одномерной равномерной сетки h_x (5.1.6) с шагом h.

Таким образом, в результате численного решения задачи (5.1.11), (5.1.12) получаем значения j_h(x₀), j_h(x₁), j_h(x₂), …, j_h(x_m) функции j_h(x) – приближенное решение задачи.

В последующих параграфах данной главы будут рассмотрены конкретные численные методы решения применительно к задаче (5.1.11), (5.1.12) – задаче Коши для одного обыкновенного дифференциального уравнения первого порядка. Вместе с тем рассматриваемые методы могут применяться и для решения задачи Коши для системы дифференциальных уравнений первого порядка (5.1.3), (5.1.4).

Рассматриваемые численные методы решения задачи Коши (5.1.11), (5.1.12) можно разделить на две группы:

– многошаговые разностные методы. К ним относятся метод Эйлера (§ 5.2) и явный 4-шаговый метод Адамса (§ 5.3);

– методы Рунге – Кутта (в § 5.4 будет рассмотрен метод Рунге – Кутта четвертого порядка).

Эти методы имеют некоторые общие свойства и характеристики.

Как следует из вышесказанного, j_h(x) – приближенное решение, полученное в результате численного решения задачи Коши, представлено значениями в узлах одномерной равномерной сетки h_x с шагом h - j_h(x₀), j_h(x₁),j_h(x₂), …, j_h(x_m). Расчет значений приближенного решения j_h(x) в узлах сетки осуществляется для всех рассматриваемых численных методов решения задачи Коши по правилу

(5.1.14)

Таким образом, отличие одного численного метода от другого заключается лишь в способе построения Dj_h(x_k), k = 0, 1, 2, …, m – 1.

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

Определение 5.1.1

Определение 5.1.2

Следует отметить, что при стремлении шага h сетки к нулю общее количество узлов m равномерной сетки h_x на отрезке [a, b] стремится к бесконечности.

Определение 5.1.3

Метод Эйлера

Пусть на отрезке [a, b] требуется найти решение задачи Коши для обыкновенного дифференциального уравнения первого порядка

y¢ = g(x, y), (5.2.1)

y(x₀) = y₀, (5.2.2)

где y₀ – заданное число.

Зададим на отрезке [a, b] равномерную сетку

h_x = {x_i / x_i = x_{i – 1} + h, h > 0, i = 1, 2, 3, …, m; x₀ = a, x_m = b}, (5.2.3)

где x_i, i = 0, 1, 2, …, m – узлы одномерной сетки h_x; h – шаг сетки.

Приближенное решение j_h(x) задачи (5.2.1) - (5.2.2) будем строить по правилу

(5.2.4)

Организация вычислений при реализации метода Эйлера предельно проста. Промежуточные результаты удобно размещать в таблице, аналогичной табл. 5.2.1.

Таблица 5.2.1

В заключение сделаем несколько замечаний относительно метода Эйлера (5.2.4).

1. Рассмотренный метод является явным одношаговым разностным. Это один из наиболее простых численных методов решения задачи Коши.

2. Сходится на отрезке [a, b] в соответствии с определением и имеет первый порядок точности.

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

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

Рассмотрим систему из двух уравнений

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

Как и ранее, будем искать решение поставленной задачи Коши на отрезке [a, b], задав на нем равномерную сетку h_x с шагом h.

Пусть функции j_{h y}(x) » y(x) и j_{h z}(x) » z(x) обозначают приближенное решение задачи, тогда метод Эйлера запишется в виде