1.10.3. Распространение ошибок в начальных данных при решении обыкновенных дифференциальных уравнений.

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

, ,

решением которой является функция y(x) = y₀ exp(lx), для анализа распространения ошибок при использовании метода Эйлера.

Приближенное начальное значение y₀ отличается от точного начального значения y(x₀) либо из-за погрешностей в задании начальных данных, либо из-за ошибок округления в представлении значения y(x₀) в ЭВМ, либо по обеим причинам сразу. Поэтому в реальном начальном значении y₀ содержится ошибка e₀ = y(x₀) - y₀. Необходимо также отметить, что применение метода Эйлера в данном случае эквивалентно использованию выражения 1 + hl для аппроксимации функции exp(lx).

Выражение для глобальной ошибки E_n+1 на (n+1)-ом шаге интегрирования запишется так:

.

Таким образом, глобальная ошибка состоит из двух компонент. Первая компонента представляет собой ошибку, вытекающую из аппроксимации значения exp(lx) выражением 1 + hl, вторая компонента отражает распространение ошибки e₀ в начальном условии. Если ç1 + hlç > 1, то вторая компонента растет при увеличении n и в конечном счете станет доминирующей частью глобальной ошибки E_n+1. Если l < 0, то, для того чтобы ограничить распространение ошибки e₀, необходимо обеспечить выполнение условия

или .

Данное условие на величину шага интегрирования h называется условием устойчивости, поскольку в случае его нарушения приближенное решение будет носить неустойчивый характер. Таким образом, условие устойчивости накладывает ограничение на максимальную величину шага интегрирования h_max, значение которого зависит как от параметров задачи, так и от применяемого метода. Условие устойчивости h < 2/çlç получено для метода Эйлера.

Для методов Адамса (3.41’) это условие записывается в виде

и приводит к таким ограничениям на шаг интегрирования:

для метода второго порядка аппроксимации hçlç < 1,

для метода третьего порядка аппроксимации hçlç < 0.545,

для метода четвертого порядка аппроксимации hçlç < 0.3,

для метода пятого порядка аппроксимации hçlç < 0.163.

Для классического метода Рунге-Кутта (2.22) четвертого порядка аппроксимации hçlç < 2.785.

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

. (4.8)

Если матрица системы дифференциальных уравнений имеет большие по модулю отрицательные собственные значения l_j, то ограничение (4.8) на шаг h является на больших интервалах интегрирования слишком обременительным, т.к. оно требует использования очень малого шага на протяжении всего интервала интегрирования. Если матрица системы имеет большой разброс собственных значений, например, в случае двух уравнений çl₁ç?çl₂ç и l₁ < 0, l₂ < 0, то первая составляющая точного решения

,

где e₁, e₂ — собственные вектора матрицы, очень быстро затухает на отрезке, длина которого имеет порядок t = 1/çl₁ç, а затем становится ничтожно малой. Именно на этом отрезке она вносит свой вклад в решение y(x).

Таким образом, на интервале интегрирования выделяются два промежутка с разным характером поведения решения. На первом промежутке, называемом пограничным слоем, для воспроизведения быстро изменяющегося решения с приемлемой точностью необходим шаг интегрирования, удовлетворяющий условию h ? 1/çl₁ç. На втором промежутке, где решение характеризуется малой скоростью изменения, шаг нельзя увеличить, т.к. при нарушении условия (4.8) соответствующая составляющая решения y_n не будет затухать.

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

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

Если же l > 0, то неравенство ç1 + hlç > 1 выполняется для любого h > 0, т.е. на каждом шаге процесса интегрирования наблюдается увеличение локальных ошибок. Однако на этот раз неустойчивость не столь очевидна, поскольку доминирующей является первая компонента глобальной ошибки E_n+1: каким бы малым ни был выбран шаг h, разность (exp[(n + 1)hl] —(1 + hl)ⁿ⁺¹) при достаточно большом n всегда будет превосходить (1 + hl)ⁿ⁺¹. Такого рода задачи называют плохо обусловленными.

Если аналогичным образом провести анализ распространения ошибок на общий случай y’ = f(x, y), то тогда имеют место следующие утверждения:

а) если f/y < 0, то влияние локальных ошибок уменьшается при значениях h, удовлетворяющих условиям устойчивости;

б) если f/y > 0, то влияние локальных ошибок увеличивается независимо от того, насколько малым выбран шаг h.

Во многих задачах знак f/y меняется на интервале интегрирования, т.е. в процессе интегрирования локальные ошибки то увеличиваются, то уменьшаются. Имеет смысл выполнять интегрирование так, чтобы на каждом шаге знать знак f/y и тем самым по крайней мере контролировать ситуацию. В простейшем случае y’ = ly, l > 0, лучше интегрировать в направлении уменьшения x, т.е. с отрицательным шагом.

21.Гантмахер Ф.Р. Теория матриц. — М.: Гл. ред. физ.-мат. лит.

изд-ва Наука, 1966. —576 с.