Метод Булирша-Штера с использованием рациональной экстраполяции для системы уравнений

2.3. Метод Булирша-Штера с использованием рациональной экстраполяции для системы уравнений

 

Метод Булирша-Штера (Bulirsch-Stoer Method) - это метод решения системы обыкновенных дифференциальных уравнений первого порядка с гладкими правыми частями. Гладкость правых частей является необходимой для работы метода. Если правые части вашей системы не являются гладкими или содержат разрывы, то лучше использовать метод Рунге-Кутта. В случае же гладкой системы метод Булирша-Штера позволяет добиться существенно большей точности, чем метод Рунге-Кутта.

 

Принцип работы метода

Основной идеей метода является вычисление состояния системы в точке x+h, как результата двух шагов длины h/2, четырех шагов длины h/4, восьми шагов длины h/8 и так далее с последующей экстраполяцией результатов. Метод строит рациональную интерполирующую функцию, которая в точке h/2 проходит через состояние системы после двух таких шагов, в точке h/4 проходит через состояние системы после четырех таких шагов, и т.д., а затем вычисляет значение этой функции в точке h = 0, проводя экстраполяцию.

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

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

2.4            Метод Адамса

Явная схема Адамса.

Рассмотренные выше методы являются явными одношаговыми (для нахождения последующего приближения используется лишь одно предыдущее). Приведённый ниже метод является многошаговым.

Пусть задана задача Коши:

(2.4.1)

Для точного решения (которое нам не известно) выполнено:

(2.4.2)

Предположим, нам известны приближенные значения  функции u(x) в k точках  (стартовые k точек, в частности, можно найти методом Эйлера или методом Рунге-Кутта того или иного порядка), тогда функцию f(x,u(x)) в (2.4.2) для приближенного вычисления интеграла можно заменить на интерполяционный полином  порядка k-1, построенный по k точкам , интеграл от которого считается явно и представляет собой линейную комбинацию значений c некоторыми множителями . Таким образом, мы получаем следующую рекуррентную процедуру вычисления приближенных значений  функции u(x) (являющимся точным решением задачи Коши) в точках  :

(2.4.3)

Описанная схема является k-шаговой явной формулой Адамса.

Неявная схема Адамса.

Пусть - интерполяционный полином порядка k, построенный по k+1 значению б одно из которых, именно , мы будем считать неизвестным. Модифицируем (2.4.3), заменив в нём  на полином более высокой степени , интеграл от которого выражается в виде линейной комбинации значений  с некоторыми новыми коэффициентами :

(2.4.4)

Формула (2.4.4) представляет собой неявную схему Адамса и является уравнением на , которое можно решать методом последовательных приближений. Естественно, что начальное приближение , должно быть разумно выбрано. Для этого удобно объединить явную и неявную схемы Адамса в одну, называемую «методом коррекции». Именно с помощью явной схемы определяется начальное приближение (прогноз), а затем по неявной схеме оно необходимое число раз (обычно один или два) корректируется методом последовательных приближений до достижения заданной точности (коррекция).

2.5.Метод Эйлера.

Решить дифференциальное уравнение у^/=f(x,y) численным методом - это значит для заданной последовательности аргументов х₀, х₁…, х_n и числа у₀, не определяя функцию у=F(x), найти такие значения у₁, у₂,…, у_n, что у_i=F(x_i)(i=1,2,…, n) и F(x₀)=y₀. (2.5.1)

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

У=F(x) получить таблицу значений этой функции для заданной последовательности аргументов. Величина h=x_k-x_k-1 называется шагом интегрирования.

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

Рассмотрим дифференциальное уравнение первого порядка (2.5.1)

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

x=x₀, y(x₀)=y₀  (2.5.2)

Требуется найти решение уравнения (2.5.1) на отрезке [а,b].

Разобьем отрезок [a, b] на n равных частей и получим последовательность х₀, х₁, х₂,…, х_n, где x_i=x₀+ih (i=0,1,…, n), а h=(b-a)/n-шаг интегрирования.

В методе Эйлера приближенные значения у(х_i)»y_i вычисляются последовательно по формулам у_i+hf(x_i, y_i) (i=0,1,2…).

При этом искомая интегральная кривая у=у(х), проходящая через точку М₀(х₀, у₀), заменяется ломаной М₀М₁М₂… с вершинами М_i(x_i, y_i) (i=0,1,2,…); каждое звено М_iM_i+1 этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения (2.5.1), которая проходит через точку М_i. Если правая часть уравнения (2.5.1) в некотором прямоугольнике R{|x-x₀|£a, |y-y₀|£b}удовлетворяет условиям:
|f(x, y₁)- f(x, y₂)| £ N|y₁-y₂| (N=const), (2.5.3)

|df/dx|=|df/dx+f(df/dy)| £ M (M=const),

то имеет место следующая оценка погрешности:

|y(x_n)-y_n| £ hM/2N[(1+hN)ⁿ-1], (2.5.4)

где у(х_n)-значение точного решения уравнения (2.5.1) при х=х_n, а у_n- приближенное значение, полученное на n-ом шаге.

Формула (13) имеет в основном теоретическое применение. На практике иногда оказывается более удобным двойной просчет: сначала расчет ведется с шагом h, затем шаг дробят и повторный расчет ведется с шагом h/2. Погрешность более точного значения у_n^* оценивается формулой

|y_n-y(x_n)|»|y_n^*-y_n|. (2.5.5)

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

Модифицированный метод Эйлера

Рассмотрим дифференциальное уравнение (2.5.1) y^/=f(x,y) с начальным условием y(x₀)=y₀. Разобьем наш участок интегрирования на n равных частей. На малом участинтегральную кривую заменим прямой  линией.

 Рис.1 Метод Эйлера в графическом видa

 Получаем точку М_к(х_к,у_к). Через М_к проводим касательную: у=у_к=f(x_k,y_k)(x-x_k). Делим отрезок (х_к,х_к1) пополам:

x_Nk^/=x_k+h/2=x_{k+1/2 (2.5.6)}

y_Nk^/=y_k+f(x_k,y_k)h/2=y_k+y_k+1/2

Получаем точку N_k^/. В этой точке строим следующую касательную:

y(x_k+1/2)=f(x_k+1/2, y_k+1/2)=α_{k (2.5.7)}

Из точки М_к проводим прямую с угловым коэффициентом α_к и определяем точку пересечения этой прямой с прямой Х_к1. Получаем точку М_к^/. В качестве у_к+1 принимаем ординату точки М_к^/. Тогда:

у_к+1=у_к+α_кh

x_k+1=x_k+h

(2.5.8) α_k=f(x_k+h/2, y_k+f(x_k,Y_k)h/2)

y_k=y_k-1+f(x_k-1,y_k-1)h

(2.5.8)-рекурентные формулы метода Эйлера.

Сначала вычисляют вспомогательные значения искомой функции у_к+1/2 в точках х_к+1/2, затем находят значение правой части уравнения (11) в средней точке y^/_k+1/2=f(x_k+1/2, y_k+1/2) и определяют у_к+1.

Для оценки погрешности в точке х_к проводят вычисления у_к с шагом h, затем с шагом 2h и берут 1/3 разницы этих значений:

| у_к^*-у(х_к)|=1/3(y_k^*-y_k), (2.5.9)

где у(х)-точное решение дифференциального уравнения.

 Таким образом, методом Эйлера можно решать уравнения любых порядков. Например, чтобы решить уравнение второго порядка y^//=f(y^/,y,x) c начальными условиями y^/(x₀)=y^/₀, y(x₀)=y₀, выполняется замена:

y^/=z (2.5.10)

z^/=f(x,y,z)

Тем самым преобразуются начальные условия: y(x₀)=y₀, z(x₀)=z₀, z₀=y^/₀. (2.5.12)

3.Описание алгоритмов решения задачи

3.1.Описание переменных.

Наименование	Тип	Описание
Входные данные
Xi	double	Начальное значение (x) интервала вычисления
Xkon	double	Конечное значение (x) интервала вычисления
n	integer	Количество шагов
Yi	double	Начальное значение y
kx	double	Коэффициент при переменой x
ky	double	Коэффициент при переменной y
Выходные данные
h	double	Фиксированное приращение аргумента (x)
res	double	Расчётное значение уравнение y’=F(x,y) в точке (x)
Промежуточные
i	integer	Счётчик цикла
Yprom	double	Промежуточное значение y в точке Xprom
Xprom	double	Промежуточное значение x при h/2
a	double	Решение уравнения в точках f(Xprom,Yprom)
f1	double	Функция f(x,y)