Как найти порядок аппроксимации

From Wikipedia, the free encyclopedia

In science, engineering, and other quantitative disciplines, order of approximation refers to formal or informal expressions for how accurate an approximation is.

Usage in science and engineering[edit]

In formal expressions, the ordinal number used before the word order refers to the highest power in the series expansion used in the approximation. The expressions: a zeroth-order approximation, a first-order approximation, a second-order approximation, and so forth are used as fixed phrases. The expression a zero-order approximation is also common. Cardinal numerals are occasionally used in expressions like an order-zero approximation, an order-one approximation, etc.

The omission of the word order leads to phrases that have less formal meaning. Phrases like first approximation or to a first approximation may refer to a roughly approximate value of a quantity.^[1][2] The phrase to a zeroth approximation indicates a wild guess.^[3] The expression order of approximation is sometimes informally used to mean the number of significant figures, in increasing order of accuracy, or to the order of magnitude. However, this may be confusing, as these formal expressions do not directly refer to the order of derivatives.

The choice of series expansion depends on the scientific method used to investigate a phenomenon. The expression order of approximation is expected to indicate progressively more refined approximations of a function in a specified interval. The choice of order of approximation depends on the research purpose. One may wish to simplify a known analytic expression to devise a new application or, on the contrary, try to fit a curve to data points. Higher order of approximation is not always more useful than the lower one. For example, if a quantity is constant within the whole interval, approximating it with a second-order Taylor series will not increase the accuracy.

In the case of a smooth function, the nth-order approximation is a polynomial of degree n, which is obtained by truncating the Taylor series to this degree. The formal usage of order of approximation corresponds to the omission of some terms of the series used in the expansion (usually the higher terms). This affects accuracy. The error usually varies within the interval. Thus the numbers zeroth, first, second etc. used formally in the above meaning do not directly give information about percent error or significant figures.

Zeroth-order[edit]

Zeroth-order approximation is the term scientists use for a first rough answer. Many simplifying assumptions are made, and when a number is needed, an order-of-magnitude answer (or zero significant figures) is often given. For example, you might say «the town has a few thousand residents», when it has 3,914 people in actuality. This is also sometimes referred to as an order-of-magnitude approximation. The zero of «zeroth-order» represents the fact that even the only number given, «a few», is itself loosely defined.

A zeroth-order approximation of a function (that is, mathematically determining a formula to fit multiple data points) will be constant, or a flat line with no slope: a polynomial of degree 0. For example,

could be – if data point accuracy were reported – an approximate fit to the data, obtained by simply averaging the x values and the y values. However, data points represent results of measurements and they do differ from points in Euclidean geometry. Thus quoting an average value containing three significant digits in the output with just one significant digit in the input data could be recognized as an example of false precision. With the implied accuracy of the data points of ±0.5, the zeroth order approximation could at best yield the result for y of ~3.7 ± 2.0 in the interval of x from −0.5 to 2.5, considering the standard deviation.

If the data points are reported as

the zeroth-order approximation results in

The accuracy of the result justifies an attempt to derive a multiplicative function for that average, for example,

One should be careful though, because the multiplicative function will be defined for the whole interval. If only three data points are available, one has no knowledge about the rest of the interval, which may be a large part of it. This means that y could have another component which equals 0 at the ends and in the middle of the interval. A number of functions having this property are known, for example y = sin πx. Taylor series is useful and helps predict an analytic solution, but the approximation alone does not provide conclusive evidence.

First-order[edit]

First-order approximation is the term scientists use for a slightly better answer.^[3] Some simplifying assumptions are made, and when a number is needed, an answer with only one significant figure is often given («the town has 4×10³, or four thousand, residents»). In the case of a first-order approximation, at least one number given is exact. In the zeroth-order example above, the quantity «a few» was given, but in the first-order example, the number «4» is given.

A first-order approximation of a function (that is, mathematically determining a formula to fit multiple data points) will be a linear approximation, straight line with a slope: a polynomial of degree 1. For example:

is an approximate fit to the data.
In this example there is a zeroth-order approximation that is the same as the first-order, but the method of getting there is different; i.e. a wild stab in the dark at a relationship happened to be as good as an «educated guess».

Second-order[edit]

Second-order approximation is the term scientists use for a decent-quality answer. Few simplifying assumptions are made, and when a number is needed, an answer with two or more significant figures («the town has 3.9×10³, or thirty-nine hundred, residents») is generally given. In mathematical finance, second-order approximations are known as convexity corrections. As in the examples above, the term «2nd order» refers to the number of exact numerals given for the imprecise quantity. In this case, «3» and «9» are given as the two successive levels of precision, instead of simply the «4» from the first order, or «a few» from the zeroth order found in the examples above.

A second-order approximation of a function (that is, mathematically determining a formula to fit multiple data points) will be a quadratic polynomial, geometrically, a parabola: a polynomial of degree 2. For example:

is an approximate fit to the data. In this case, with only three data points, a parabola is an exact fit based on the data provided. However, the data points for most of the interval are not available, which advises caution (see «zeroth order»).

Higher-order[edit]

While higher-order approximations exist and are crucial to a better understanding and description of reality, they are not typically referred to by number.

Continuing the above, a third-order approximation would be required to perfectly fit four data points, and so on. See polynomial interpolation.

Colloquial usage[edit]

These terms are also used colloquially by scientists and engineers to describe phenomena that can be neglected as not significant (e.g. «Of course the rotation of the Earth affects our experiment, but it’s such a high-order effect that we wouldn’t be able to measure it.» or «At these velocities, relativity is a fourth-order effect that we only worry about at the annual calibration.») In this usage, the ordinality of the approximation is not exact, but is used to emphasize its insignificance; the higher the number used, the less important the effect. The terminology, in this context, represents a high level of precision required to account for an effect which is inferred to be very small when compared to the overall subject matter. The higher the order, the more precision is required to measure the effect, and therefore the smallness of the effect in comparison to the overall measurement.

See also[edit]

• Linearization
• Perturbation theory
• Taylor series
• Chapman–Enskog method
• Big O notation

References[edit]

2.2. Порядок аппроксимации. Метод решения

Для
определения порядка аппроксимации
явной разностной схемы (6.2) подставим в
неё выражения (2.16)–(2.18), описывающие
разложение значений

в ряд Тейлора относительно точки

на разностной сетке:

Таким
образом, явная разностная схема (6.2)
аппроксимирует исходное дифференциальное
уравнение (6.1) с первым порядком и по
времени, и по координате:

Разностный
шаблон (см. рисунок), характеризующий
явную разностную схему (6.2), свидетельствует
о том, что она содержит одну неизвестную
величину – значение функции u
на (n + 1)-ом
шаге по времени. Выражая эту величину
из разностной схемы, получаем рекуррентное
соотношение

(6.4)

позволяющее
рассчитать все значения функции u
на (n + 1)-ом
шаге по времени (при известных значениях
функции u на n-ом
шаге), кроме значений
,
определяемых с помощью граничных
условий. Если заданы граничные условия
1-го рода, то значения

определяются непосредственно из их
разностной аппроксимации; если 2-го или
3-го рода, то – с помощью соотношений
(4.4a) и (4.4b).
Таким образом, алгоритм решения явной
разностной схемы (6.2) аналогичен алгоритму
решения явной разностной схемы (4.2),
аппроксимирующей дифференциальное
уравнение параболического типа, не
содержащее производную по координате
первого порядка.

3. Неявная разностная схема

3.1. Характеристика

Неявная
разностная схема для уравнения (6.1) имеет
вид:

(6.5)

Учитывая
порядок аппроксимации разностных
операторов, составляющих данную
разностную схему, легко видеть, что она
аппроксимирует дифференциальное
уравнение (6.1) с первым порядком и по
времени, и по координате:

Исследуем
устойчивость неявной разностной схемы
(6.5) с помощью спектрального метода. Для
этого отбрасываем член
,
наличие которого, как известно, не
оказывает влияния на устойчивость
разностной схемы, и представляем решение
в виде гармоники (3.7):

Далее,
упрощаем полученное выражение, деля
левую и правую его части на
:

Используя
зависимости (3.9), (3.10), получаем

Группируя
члены, содержащие ,
в левой части уравнения, выразим величину,
обратную :

При
этом необходимое условие устойчивости
разностных схем (3.8) также преобразуем
к виду:

(6.6)

Комплексный
вид полученного выражения свидетельствует
о том, что для устойчивости разностной
схемы (6.5) согласно условию (6.6) требуется,
чтобы величины, обратные собственным
числам оператора перехода, были
расположены вне или на границе круга
радиусом 1, центр которого находится в
начале координат комплексной плоскости.

Введём
следующие обозначения:

Следовательно,
величины, обратные собственным числам
оператора перехода, расположены на
комплексной плоскости на окружности с
центром в точке

и радиусом r (см.
рисунок). Данная окружность находится
вне круга, соответствующего условию
устойчивости, при любом значении r.
Это означает, что неявная разностная
схема (6.5) абсолютно устойчива.

3. Неявная разностная схема

3.2. Метод решения

Разностный
шаблон (см. рисунок), характеризующий
неявную разностную схему (6.5), свидетельствует
о том, что она содержит три неизвестные
величины – значения функции u
на (n + 1)-ом
шаге по времени. Следовательно, для
решения данной разностной схемы
необходимо использовать метод прогонки.

Приведём
разностную схему (6.5) к виду (4.10), удобному
для использования метода прогонки:

Следовательно,
коэффициенты, соответствующие уравнению
(4.10), имеют вид:

Легко
видеть, что для разностной схемы (6.5)
достаточное условие сходимости прогонки
(4.16) выполняется:

Алгоритм
решения, а также методики определения
прогоночных коэффициентов и решения
на правой границе аналогичны описанным
ранее.

Разностные схемы для решения задачи Коши

Рассмотрим основные принципы конструирования разностных схем.

1. Производная [math]left.{left(frac{dy}{dx}right)!}right|_{x=x_{i}}[/math] аппроксимируется на заранее выбранном шаблоне [math](x_{i-k+1}, x_{i-k+2},ldots,x_{i+1})[/math], а правая часть [math]f(x,y)[/math] обыкновенного дифференциального уравнения заменяется сеточным представлением [math]f(x_{i},widehat{y}_{i})[/math] в точке [math](x_{i},widehat{y}_{i})[/math]. Условно этот принцип называется аппроксимационным. С его помощью строятся одношаговые и многошаговые схемы в зависимости от шаблона, по которому аппроксимируется производная.

2. Дифференциальное уравнение заменяется разностным аналогом на одношаговом (двухточечном) шаблоне [math](x_{i},x_{i+1})[/math] (возможно использование дополнительных точек внутри шаблона) путем согласования с разложением функции [math]y=y(x)[/math] при [math]x=x_{i+1}[/math] no формуле Тейлора относительно точки [math]x=x_{i}[/math]. По этому принципу строятся одношаговые схемы, в частности семейство схем Рунге-Кутты.

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

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

---

Принцип аппроксимаций

Рассмотрим проблему нахождения численного решения задачи Коши (6.9): [math]frac{dy}{dx}= f(x,y),~ y(x_0)=y_0[/math]. Для аппроксимации производной [math]left.{left(frac{dy}{dx}right)!}right|_{x=x_{i}}[/math] используем формулу (5.4), записанную на двухточечном шаблоне [math]H_{2,i}= (x_{i},x_{i+1})colon[/math]

и формулу (5.10), записанную на трехточечном нерегулярном шаблоне [math]H_{3,i}= (x_{i-1}, x_{i}, x_{i+1})colon[/math]

где [math]H_{i}^{i+1}= h_{i}+h_{i+1},~ Pi_{i}^{i+1}= h_{i}h_{i+1},~ h_{i+1}=x_{i+1}-x_{i},~~ delta_{i+1}= frac{h_{i+1}}{h_{i}}[/math] — параметр нерегулярности.

При [math]delta_{i+1}=1[/math] шаблон становится регулярным, а среднее слагаемое с входящим в него [math]y_{i}[/math] — равным нулю. При этом формула (6.17) сводится к традиционной:

Аппроксимация формульных функций, входящих в правую часть (6.9), осуществляется путем перехода к их сеточному представлению.

---

Методика построения разностных схем с помощью аппроксимаций производной

1. Вводится в общем случае неравномерная сетка [math]Omega_n= (x_0, x_1,ldots, x_{i}, x_{i+1},ldots, x_{n})[/math]. Величина шага [math]h_{i+1}=x_{i+1}-x_{i}[/math] выражается через узловые точки.

2. Для получения рассматриваемых здесь двух схем (одношаговой и двух-шаговой) выбирается соответственно двухточечный (одношаговый) и трехточечный (двухшаговый) шаблоны [math]H_{2,i}= (x_{i},x_{i+1})[/math] и [math]H_{3,i}= (x_{i-1},x_{i},x_{i+1})[/math], на которых в точке [math]x_{i}[/math] производится аппроксимация производной [math]left.{left(frac{dy}{dx} right)!}right|_{x=x_{i}}[/math] по формулам. (6.16), (6.17) соответственно в левой и центральной точках [math]x_{i}[/math]. Далее заменяется правая часть уравнения (6.9) ее сеточным представлением, т.е. [math]f(x,y)to f(x_{i},y_{i})[/math], а вместо [math]y(x)[/math] рассматривается сеточная функция [math]widehat{y}_{i}approx y(x_{i})[/math], которая определяется только в точках сетки.

3. Выполняется подстановка аппроксимаций (6.16), (6.17) и [math]f(x_{i},y_{i})[/math] в дифференциальное уравнение:

[math]begin{gathered}frac{y_{i+1}-y_{i}}{h_{i+1}}+O(h_{i+1})= f(x_{i},y_{i});\ frac{1}{H_{i}^{i+1}}! left[-delta_{i+1}y_{i-1}+ frac{delta_{i+1}^2-1}{delta_{i+1}}y_{i}+ frac{y_{i+1}}{delta_{i+1}}right]+ O(h_{i}h_{i+1})= f(x_{i},y_{i}). end{gathered}[/math]

После отбрасывания остаточных слагаемых получаются разностные схемы: — явная схема Эйлера первого порядка (явный метод Эйлера):

– обобщенная на нерегулярный шаблон двухшаговая явная схема Эйлера второго порядка (где [math]Delta widehat{y}_{i}=widehat{y}_{i+1}-widehat{y}_{i}[/math]):

Обозначим схему (6.20) символами 2Y2B (суть этого обозначения поясняется далее). Анализ остаточного слагаемого [math]left(frac{h_{i}^2}{6}delta_{i+1}M_{3,i}right)[/math] обобщенной схемы Эйлера указывает на то, что при условном выборе шаблона, когда достигается соотношение [math]delta_{i+1}leqslant h_{i}[/math] (то есть [math]h_{i+1}leqslant h_{i}^2[/math]) , порядок аппроксимации этой схемы повышается на единицу без увеличения количества точек шаблона. Этот прием повышения точности может быть использован при конструировании вложенных алгоритмов путем измельчения шага.

При [math]delta_{i+1}=1[/math] (регулярный шаблон) схема (6.20) соответствует методу Эйлера-Коши:

Для начала расчетов по формулам (6.20),(6.21) требуется иметь две «разгонные» точки [math]widehat{y}_{0},,widehat{y}_{1}[/math]. Первая определяется известным начальным условием [math]widehat{y}_0=y_0[/math], а вторая может быть найдена с помощью другого метода, например, по формуле (6.19): [math]widehat{y}_1=y_0+h_1f(x_0,y_0)[/math].

Приведем геометрическую интерпретацию явного метода Эйлера (6.19), которая показана на рис. 6.5.

Пусть известна точка [math](x_{i},widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем к кривой [math]y=y(x)[/math] касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i}, widehat{y}_{i})[/math]. Ее уравнение имеет вид

[math]y=widehat{y}_{i}+ y'(x_{i})(x-x_{i})= widehat{y}_{i}+ f(x_{i}, widehat{y}_{i})(x-x_{i}).[/math]

Следующая точка [math]widehat{y}_{i+1}[/math] приближенного решения получается при [math]x=x_{i+1}[/math] (в результате нахождения точки пересечения прямых [math]L[/math] и [math]x_{i+1}=text{const}[/math]): [math]widehat{y}_{i+1}= widehat{y}_{i}+ h_{i+1}f(x_{i},widehat{y}_{i})[/math]. Это соотношение совпадает с формулой (6.19).

Схема (6.19) иногда называется схемой ломаных (рис. 6.6), так как очередное значение [math]widehat{y}_{i+1}[/math] получается в результате проведения касательной к интегральной кривой в точке [math](x_{i},widehat{y}_{i})[/math]. После выполнения каждого шага метода Эйлера приближенное решение переходит с одной интегральной кривой на другую, что обусловлено погрешностью численного решения. Большая погрешность схемы связана с тем, что порядок этой схемы равен единице, и поэтому ломаная, определяемая по (6.19), довольно значительно отклоняется от точного решения [math]y=y(x)[/math]. В связи с этим данная схема используется, как правило, в качестве вспомогательной при конструировании схем повышенного порядка.

Для анализа устойчивости явного метода Эйлера применим формулу (6.19) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате разностное уравнение будет иметь вид [math]widehat{y}_{i+1}= widehat{y}_{i}+ h,mu,widehat{y}_{i}[/math] или [math]widehat{y}_{i+1}-(1+h,mu)widehat{y}_{i}=0[/math]. Корень соответствующего характеристического уравнения [math]lambda_1=1+h,mu[/math]. Условие (6.13) принимает форму [math]|1+h,mu|<1[/math]. Отсюда следует, что область устойчивости представляется внутренней частью круга единичного радиуса с центром в точке [math](-1;0)[/math] (рис. 6.7). Если [math]mu[/math] действительная константа, то можно указать интервал устойчивости [math]-2<h,mu<0[/math], то есть [math]h,muin (-2;0)[/math], а при [math]mu<0[/math] имеем [math]0<h<-frac{2}{mu}[/math]. Поэтому явный метод Эйлера является ограниченно устойчивым с критическим шагом [math]h_{text{kr}}=-frac{2}{mu}[/math].

Дадим геометрическую интерпретацию метода Эйлера-Коши (6.21) применительно к решению задачи (6.9) (рис. 6.8). Пусть известны две точки [math](x_{i-1},widehat{y}_{i-1}),, (x_{i}, widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Проведем через точку [math](x_{i},widehat{y}_{i})[/math] касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i},widehat{y}_{i})[/math]. После этого через точку [math](x_{i-1}, widehat{y}_{i-1})[/math] проведем прямую [math]L_1[/math], параллельную [math]L[/math]. Ее уравнение имеет вид [math]y=widehat{y}_{i-1}+ f(x_{i},widehat{y}_{i})(x-x_{i-1})[/math]. Следующая точка [math]widehat{y}_{i+1}[/math], приближенного решения получается при [math]x=x_{i+1}colon[/math] [math]widehat{y}_{i+1}= widehat{y}_{i-1}+2h f(x_{i}, widehat{y}_{i})[/math]. Это соотношение совпадает с формулой (6.21).

Замечание. Все описанные здесь и далее методы легко переносятся на системы обыкновенных дифференциальных уравнений (6.7). Под [math]y[/math] и [math]f(x,y)[/math] следует подразумевать векторы [math]Y=(y_1,ldots,y_n)^T[/math] и [math]F(x,Y)= bigl(f_1(x,Y), ldots, f_n(x,Y)bigr)^T[/math]. Например, для системы двух дифференциальных уравнений

[math]begin{cases}y’=f(x,y,z),& y(x_0)=y_0,\ z’=g(x,y,z),& z(x_0)=z_0 end{cases}[/math]

явный метод Эйлера (6.19) записывается в форме

[math]left{! begin{aligned}&widehat{y}_{i+1}= widehat{y}_{i}+ h_{i+1}cdot f(x,widehat{y}_{i},widehat{z}_{i}), &~ & widehat{y}_0=y_0,\ &widehat{z}_{i+1}= widehat{z}_{i}+ h_{i+1}cdot g(x,widehat{y}_{i},widehat{z}_{i}), &~ & widehat{z}_0=z_0. end{aligned}right.[/math]

Пример 6.3. Исследовать устойчивость явного метода Эйлера на примере решения уравнения [math]Ty’+y=1,[/math] [math]y(0)=y_0[/math], где [math]T[/math] — действительное положительное число. Найти приближенное решение задачи Коши [math]0,!1y’+y=1,~y(0)=0[/math] на отрезке [math][0;1][/math] явным методом Эйлера.

Решение

Замечание. Порядок точности метода, как правило, определяется порядком аппроксимации схем. Например, для схемы (6.19) первого порядка, записанной в виде, соответствующем исходному уравнению (6.9), получаем (первое слагаемое аппроксимирует производную):

[math]frac{y(x_{i+1})-y(x_{i})}{h_{i+1}}-f(x_{i},y_{i})= frac{y_{i}+ h_{i+1}y’_{i}+ O(h_{i+1}^2)-y_{i}}{h_{i+1}}-f_{i}= y’_{i}-f_{i}+ O(h_{i+1})=0,[/math]

то есть [math]y’_{i}=f_{i}+O(h_{i+1})[/math]. Здесь были использованы формула Тейлора и свойства, приведенные ранее.

Для двухшаговой схемы (6.21) аналогично получается второй порядок аппроксимации. Действительно, воспользуемся частным случаем формулы Тейлора (см. (В. 18)) при

[math]x_{0}=x_{i},quad x_{i+1}=x_{i}+h,quad x_{i-1}=x_{i}-h_{i},quad f(x)= y(x)[/math]

имеем

[math]y_{i+1}= y_{i}+ hy’_{i}+ frac{h^2}{2}y»_{i}+ O(h^3),qquad y_{i-1}= y_{i}-hy’_{i}+ frac{h^2}{2}y»_{i}+ O(h^3).[/math]

Подставляя полученные соотношения в схему (6.21), записанную в виде [math]frac{y_{i+1}-y_{i-1}}{2h}-f(x_{i},y_{i})=0[/math], соответствующем дифференциальному уравнению (6.9), с учетом свойств, приведенных ранее, имеем

[math]frac{y_{i}+ hy’_{i}+ dfrac{h^2}{2}y»_{i}+ O(h^3)-y_{i}+ hy’_{i}-dfrac{h^2}{2}y»_{i}-O(h^3)}{2h}-f_{i}= y’_{i}+ O(h^2)-f_{i}=0[/math], то есть [math]y’_{i}=f_{i}+O(h^2)[/math].

Также легко проверить, что и схема (6.20) на нерегулярном шаблоне имеет второй порядок при безусловной аппроксимации. Если же [math]h_{i+1}leqslant h_{i}^2[/math] (условная аппроксимация), то этот порядок равен трем.

---

Интегрально-интерполяционный принцип

В данном подходе явные или неявные методы соответствующего порядка получаются путем интегрирования уравнения (6.9) вдоль решения [math]y(x)[/math]. Так, на двухточечном шаблоне

Функция [math]f(x,y(x))[/math] в (6.22) заменяется интерполяционным многочленом соответствующей степени [math]p[/math] ([math]p=0[/math] для метода первого порядка, [math]p=1[/math] для метода второго порядка, [math]p=2[/math] для метода третьего порядка и т.д.). Подчеркнем, что замена [math]f(x,y)[/math] интерполяционным многочленом в (6.22) становится возможной благодаря тому, что [math]f(x,y)[/math] рассматривается не на всей плоскости [math](x,y)[/math], а только вдоль интегральной кривой [math]y(x)[/math]. Для вычисления интеграла в правой части могут быть использованы известные квадратурные формулы прямоугольников, трапеций, парабол (Симпсона) и др.

Примем [math]x_{i+1}[/math] в качестве узловой точки многочлена [math]L_{0}(x)[/math], то есть [math]L_0(x)= f(x_{i+1},widehat{y}_{i+1})[/math]. Тогда из (6.22) следует формула неявного метода Эйлера:

Принимая для интерполяционного многочлена [math]L_1(x)[/math] в качестве узловых точек [math]x_{i}[/math] и [math]x_{i+1}[/math] и используя квадратурную формулу трапеций (5.46), получаем неявную одношаговую схему второго порядка (метод трапеций):

где [math]f_{i}= f(x_{i},widehat{y}_{i+1})[/math]. Формулы (6.23),(6.24) имеют вид (6.11) при [math]k=1[/math]. Подчеркнем, что свойство неявности схем (6.23),(6.24) обусловлено наличием искомой величины [math]widehat{y}_{i+1}[/math] в левой и правой частях в общем случае нелинейных уравнений.

При реализации алгоритма решения задачи Коши неизвестное значение [math]widehat{y}_{i+1}[/math] вычисляется одним из методов решения нелинейных уравнений. Применение метода Ньютона связано с записью уравнения в форме

и с дифференцированием функции [math]F(widehat{y}_{i+1})[/math], что увеличивает время расчетов из-за возможной сложности производных или увеличивает погрешность в силу неточного задания правых частей.

Как правило, используется метод простых итераций:

Таким образом, уравнение (6.24) решается с помощью формулы (рекуррентного соотношения):

При применении методов Ньютона и простых итераций вначале задается или находится нулевое приближение решения по формуле [math]y_{i+1}^{,(0)}= widehat{y}_{i}[/math], (так называемый «постоянный» прогноз) или явным методом Эйлера (6.19): [math]widehat{y}_{i+1}^{,(0)}= widehat{y}_{i}+ h_{i+1} f(x_{i}, widehat{y}_{i})[/math].

Итерации завершаются при выполнении условия окончания [math]bigl|widehat{y}_{i+1}^{,(k+1)}-widehat{y}_{i+1}^{,(k)}bigr| leqslant varepsilon[/math], где [math]varepsilon[/math] — малое положительное число.

Для сокращения описания и удобства ссылок получаемые здесь схемы обозначаются буквенно-цифровыми символами. Первая цифра в нем обозначает «шаговость» схемы, следующая за ней буква «Я» или буквы «НЯ» — явный или неявный характер, затем цифра указывает порядок точности схемы и за ней буква — возможную модификацию. Таким образом, схема (6.24) обозначается 1НЯ2, т.е. она является одношаговой, неявной и имеет второй порядок точности. Последняя буква указывается только для тех схем, которые имеют модификации.

Приведем геометрическую интерпретацию неявного метода Эйлера (6.23) (рис.6.10). Пусть известна точка [math](x_{i}, widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем прямую, тангенс угла наклона которой равен [math]y'(x_{i+1})= f(x_{i+1}, widehat{y}_{i+1})[/math], т.е. следующая точка [math]widehat{y}_{i+1}[/math] приближенного решения берется на той кривой из семейства интегральных кривых, касательная к которой в точке [math](x_{i+1},widehat{y}_{i+1})[/math] проходит через точку [math](x_{i},widehat{y}_{i})[/math]. Поэтому [math]frac{widehat{y}_{i+1}-widehat{y}_{i}}{h_{i+1}}= f(x_{i+1}, widehat{y}_{i+1})[/math]. Полученное соотношение совпадает с формулой (6.23).

Исследуем устойчивость неявного метода Эйлера и метода трапеций. Применим формулу (6.23) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате получим разностное уравнение [math]widehat{y}_{i+1}= widehat{y}_{i}+hcdotmucdot widehat{y}_{i+1}[/math] или [math]widehat{y}_{i+1}(1-hcdotmu)-widehat{y}_{i}=0[/math]. Корень характеристического уравнения [math]lambda_1= frac{1}{1-hcdotmu}[/math]. Условие устойчивости (6.13) имеет вид [math]left|frac{1}{1-hcdotmu}right|<1[/math] или [math]|1-hcdotmu|>1[/math]. Область устойчивости представляется внешней частью круга единичного радиуса с центром в точке [math](1;0)[/math]. Из геометрической интерпретации (рис. 6.11) следует, что неявный метод Эйлера обладает свойством A-устойчивости (см. рис.6.4,а), поскольку выделенная на рис. 6.11 область включает левую полуплоскость.

Применим формулу (6.24) метода трапеций при [math]h_{i+1}=h[/math] к решению задачи (6.12). Получим разностное уравнение

[math]widehat{y}_{i+1}= widehat{y}_{i}+ frac{h}{2} (mu widehat{y}_{i}+ mu widehat{y}_{i+1})[/math] или [math]widehat{y}_{i+1}! left(1-frac{hmu}{2}right)-left(1+ frac{hmu}{2}right)! widehat{y}_{i}=0[/math].

Корень характеристического уравнения [math]lambda_1= left(1+ frac{hmu}{2}right),colon left(1-frac{hmu}{2}right)[/math]. Условие устойчивости (6.13) принимает вид [math]left|left(1+ frac{hmu}{2}right),colon left(1-frac{hmu}{2}right)right|<1[/math] или [math]|2+mmu|< |2-hmu|[/math]. Полученное неравенство выполняется в области, изображенной на рис. 6.4,а, т.е. метод трапеций является A-устойчивым.

Продолжим рассмотрение интегрально-интерполяционного метода. Использование квадратурной формулы прямоугольников дает модифицированный метод Эйлера второго порядка:

[math]widehat{y}_{i+1!not{phantom{|}},,2}= widehat{y}_{i}+ frac{h_{i+1}}{2} f(x_{i}, widehat{y}_{i}),quad i=overline{0,n-1},[/math]

Дадим геометрическую интерпретацию модифицированного метода Эйлера второго порядка (рис. 6.12). Пусть известна точка [math](x_{i},widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Через эту точку проведем касательную [math]L[/math], тангенс угла наклона которой равен [math]y'(x_{i})= f(x_{i},widehat{y}_{i})[/math]. При [math]x=x_{i}+ frac{h_{i+1}}{2}[/math] получаем точку [math]widehat{y}_{i+1!not{phantom{|}},,2}[/math]. Далее проведем прямую [math]L_1[/math] с тангенсом угла наклона [math]f! left(x_{i}+ frac{h_{i+1}}{2}, widehat{y}_{i+1!not{phantom{|}},,2}right)[/math], а затем прямую [math]L_2[/math], параллельную [math]L_1[/math] и проходящую через точку [math](x_{i}, widehat{y}_{i})[/math]. Ее уравнение имеет вид [math]y=widehat{y}_{i}+ f! left(x_{i}+ frac{h_{i+1}}{2}, widehat{y}_{i+1!not{phantom{|}},,2}right)!cdot (x-x_{i})[/math]. Следующая точка [math]widehat{y}_{i+1}[/math] приближенного решения находится при [math]x=x_{i+1}[/math] (это соотношение совпадает с формулой (6.28)):

[math]widehat{y}_{i+1}= widehat{y}_{i}+ h_{i+1} f! left(x_{i}+ frac{h_{i+1}}{2}, widehat{y}_{i+1!not{phantom{|}},,2}right)!.[/math]

Исследуем устойчивость модифицированного метода Эйлера второго порядка. Применим формулы (6.28) при [math]h_{i+1}=h[/math] к решению задачи (6.12). В результате получим разностное уравнение:

[math]widehat{y}_{i+1}= widehat{y}_{i}+ hmu left(widehat{y}_{i}+ frac{h}{2}mu widehat{y}_{i}right)[/math] или [math]widehat{y}_{i+1}-widehat{y}_{i} left(1+ hmu+ frac{(hmu)^2}{2}right)=0[/math].

Корень соответствующего характеристического уравнения [math]lambda_1= 1+hmu+ frac{(hmu)^2}{2}[/math]. Условие устойчивости (6.13) имеет вид [math]left|1+hmu+ frac{(hmu)^2}{2}right|<1[/math]. Если [math]mu[/math] — действительное число, то, раскрывая модуль, получаем [math]-1<1+hmu+ frac{(hmu)^2}{2}<1[/math]. Неравенство [math]2+hmu+ frac{(hmu)^2}{2}>0[/math] выполняется для любых [math]hmu[/math], а неравенство [math]hmu+ frac{(hmu)^2}{2}<0[/math], то есть [math]hmucdot (2+hmu)<0[/math], справедливо при [math]-2<hmu<0[/math]. Таким образом, интервал устойчивости [math]hmuin (-2;0)[/math] модифицированного метода Эйлера совпадает с интервалом устойчивости явного метода Эйлера (6.19).

Для использования регулярного трехточечного шаблона [math](x_{i-1},x_{i},x_{i+1})[/math] заменим (6.22) аналогичным выражением: [math]textstyle{intlimits_{y_{i-1}}^{y_{i+1}} dy= intlimits_{x_{i-1}}^{x_{i+1}} f(x,y(x))dx}[/math]. Если функцию [math]f(x,y(x))[/math], входящую в правую часть, заменить многочленом Лагранжа второй степени и применить формулу (5.47), то получается неявная схема парабол (Симпсона) четвертого порядка (2НЯ4):

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

В соответствии с классификацией, схема (6.29) относится к классу неявных схем Адамса. Для получения других схем Адамса следует переписать (6.22) в форме

и заменить [math]y'(x)[/math] соответствующим интерполяционным полиномом, построенным на регулярном шаблоне с [math]h=text{const}[/math].

В силу формулы (4.36) второго интерполяционного многочлена Ньютона получаем (где [math]widehat{q}= frac{x-x_{i}}{h}[/math] — фаза интерполяции)

[math]y'(widehat{q})= y’_{i}+ widehat{q} Delta y’_{i-1}+ frac{Delta^2 y’_{i-2}}{2!} widehat{q} (widehat{q}+1)+ frac{Delta^3 y’_{i-3}}{3!} widehat{q} (widehat{q}+1)(widehat{q}+2)+ldots,[/math]

или

[math]y'(widehat{q})= y’_{i}+ widehat{q} Delta y’_{i-1}+ frac{widehat{q}^2+ widehat{q}}{2} Delta^2 y’_{i-2}+ frac{widehat{q}^3+ 3widehat{q}^2+ 2widehat{q}}{6} Delta^3 y’_{i-3}+ldots[/math]

Подставляя полученное выражение в правую часть (6.30) и учитывая, что [math]dx= hcdot dwidehat{q}[/math], вычисляем [math]textstyle{intlimits_{0}^{1} y'(widehat{q})cdot h, dwidehat{q}}[/math], поскольку при [math]x=x_{i}[/math] и [math]x=x_{i+1}[/math] справедливо [math]widehat{q}=0[/math] и [math]widehat{q}=1[/math] соответственно. В результате имеем

[math]widehat{y}_{i+1}-widehat{y}_{i}= hcdot y’_{i}+ frac{1}{2}hcdot Delta y’_{i-1}+ frac{5}{12}hcdot Delta^2 y’_{i-2}+ frac{3}{8}hcdot Delta^3 y’_{i-3}+ldots[/math]

или с учетом (6.9) общую формулу семейства явных схем Адамса:

где [math]eta_{i}= hcdot f(x_{i}, widehat{y}_{i})= hcdot f_{i}[/math], а конечные разности вычисляются с помощью табл. 6.2. Подчеркнем, что здесь интерполяционный многочлен записан на шаблоне, не включающем точку [math]x_{i+1}[/math].

Если взять в правой части (6.31) два слагаемых, то получится явный метод Эйлера (6.19), если три слагаемых, то явный метод, описываемый формулой

[math]widehat{y}_{i+1}= widehat{y}_{i}+ hcdot f_{i}+ frac{1}{2} (hcdot f_{i}-hcdot f_{i-1})= widehat{y}_{i}+ frac{h}{2} (3f_{i}-f_{i-1}).[/math]

а если четыре слагаемых, то явный метод вида

[math]begin{aligned}widehat{y}_{i+1}&= widehat{y}_{i}+ frac{h}{2} (3f_{i}-f_{i-1})+ frac{5}{12} bigl[(hcdot f_{i}-hcdot f_{i-1})-(hcdot f_{i-1}-hcdot f_{i-2})bigr]=\ &=widehat{y}_{i}+ frac{h}{12} bigl(23f_{i}-16f_{i-1}+ 5f_{i-2}bigr) end{aligned}[/math]

В результате имеем схемы Адамса-Бэшфорта

– второго порядка (2Я2В):

– третьего порядка (ЗЯЗА):

– четвертого порядка (4Я4):

Для начала расчетов по формуле (6.32) требуются две «разгонные» точки: [math]widehat{y}_{0},, widehat{y}_{1}[/math], по формуле (6.33) — три «разгонные» точки: [math]widehat{y}_{0},, widehat{y}_{1},, widehat{y}_{2}[/math], по формуле (6.34) — четыре «разгонные» точки: [math]widehat{y}_{0},, widehat{y}_{1},, widehat{y}_{2},, widehat{y}_{3}[/math]. Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Далее схема (6.32) обобщается на нерегулярный шаблон.

Если для построения интерполяционного многочлена использовать шаблоны, включающие точку [math]x_{i+1}[/math], то получится следующая формула семейства неявных схем Адамса:

где [math]eta_{i+1}= hcdot f(x_{i+1}, widehat{y}_{i+1})= hcdot f_{i+1}[/math], а конечные разности вычисляются с помощью табл. 6.3.

Если взять в правой части (6.35) два слагаемых, то получится неявный метод Эйлера (6.23), если три слагаемых, — метод трапеций (6.24):

[math]widehat{y}_{i+1}= widehat{y}_{i}+ hcdot f_{i+1}-frac{1}{2} (hcdot f_{i+1}-hcdot f_{i})= widehat{y}_{i}+ frac{h}{2} bigl[f_{i}+ f(x_{i+1}, widehat{y}_{i+1})bigr],[/math]

если четыре слагаемых, — неявный метод, описываемый соотношением (и т.д.)

[math]begin{aligned}widehat{y}_{i+1}&= widehat{y}_{i}+ frac{h}{2} (f_{i+1}+f_{i+1})-frac{1}{12} bigl[(hcdot f_{i+1}-hcdot f_{i})-(hcdot f_{i}-hcdot f_{i-1})bigr]=\ &=widehat{y}_{i}+ frac{h}{12} bigl[-f_{i-1}+ 8f_{i}+ 5f(x_{i+1}, widehat{y}_{i+1})bigr]. end{aligned}[/math]

В результате получаются неявные схемы Адамса-Мултона:

– первого порядка (6.23);

– второго порядка (6.24);

– третьего порядка (2НЯЗ):

– четвертого порядка (ЗНЯ4):

Для расчетов по формуле (6.36) требуются две «разгонные» точки: [math]widehat{y}_{0},, widehat{y}_{1}[/math], по формуле (6.37) — три «разгонные» точки: [math]widehat{y}_{0},, widehat{y}_{1},, widehat{y}_{2}[/math]. Их необходимо вычислить с порядком точности не меньше порядка точности схемы. Чтобы найти искомое значение [math]widehat{y}_{i+1}[/math] с помощью формул (6.36),(6.37), так же как в неявном методе Эйлера и методе трапеций, требуется решить в общем случае нелинейное уравнение.

Замечание. Аналогичные методы могут быть получены из интегрального уравнения вида [math]textstyle{intlimits_{y_{i-k}}^{y_{i+1}} dy= intlimits_{x_{i-k}}^{x_{i+1}} f(x,y(x)),dx,~ k geqslant 2}[/math].

---

Принцип согласования с разложением по формуле Тейлора

Рассмотрим задачу Коши (6.9) в общем случае на неравномерной сетке [math]Omega_{n}[/math]. Для получения формулы нахождения последующего значения [math]widehat{y}_{i+1}[/math] на шаблоне [math](x_{i}, x_{i+1})[/math] разложим решение [math]y(x)[/math] при [math]x=x_{i+1}[/math] относительно точки [math]x_{i}[/math] по формуле Тейлора на промежутке [math]x_{i} leqslant x leqslant x_{i+1}colon[/math]

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

[math]y’=f(x,y),quad y»= frac{d}{dx}f(x,y)= f_x+f_ycdot y’= f_x+fcdot f_y,quad ldots[/math]

Если функция [math]f(x,y)[/math] имеет p-е непрерывные производные по обоим аргументам, то в разложении (6.38) можно учесть слагаемые до [math]O(h_{i+1}^{p+1})[/math]. Чем больше членов разложения используется для вычисления, тем точнее будет приближение и тем выше порядок схемы. Однако рассчитывать производные путем дифференцирования правой части уравнения (6.9) невыгодно. Поэтому Рунге и Кутта независимо друг от друга предложили более удобную и рациональную форму представления численного решения, при использовании которой дифференцировать функцию [math]f(x,y)[/math] не требуется.

В общем случае методы Рунге-Кутты имеют следующую структуру (где [math][/math]):

[math]widehat{y}_{i+1}= widehat{y}_{i}+ h_{i+1}cdot bigl[b_1K_{1,i}+ b_2K_{2,i}+ ldots+ b_sK_{s,i}bigr],quad widehat{y}_{0}= y_0,quad i=overline{0,n-1},[/math]

где

где [math]s[/math] — число стадий (этапов), [math]K_{s,i}[/math] — значения коэффициентов схемы Рунге-Кутты, вычисленные на основе правой части уравнения (6.9), [math]c_{j},~ j= overline{s,s};[/math] [math]a_{ell,m},~ ell= overline{2,s};[/math] [math]m=overline{1,s-1};[/math] [math]b_k,~ k=overline{1,s}[/math]. Первый индекс в обозначениях коэффициентов является порядковым номером, а второй соответствует индексу точки [math]x_{i}[/math] — началу отрезка [math][x_{i},x_{i+1}][/math], на котором производится расчет.

В некоторых методах кроме вычисления приближенного решения [math]widehat{y}_{i+1}[/math] определяется еще дополнительное значение [math]widetilde{y}_{i+1}[/math] по формуле

порядок которого, как правило, на единицу больше или меньше обеспечиваемого выражением для [math]widehat{y}_{i+1}[/math]. Величина [math]|widehat{y}_{i+1}-widetilde{y}_{i+1}|[/math] служит для учета погрешности и управления величиной шага.

Коэффициенты схемы подбираются так, чтобы выражение для [math]widehat{y}_{i+1}[/math] согласовывалось с формулой Тейлора (6.38) вплоть до членов степени [math]h_{i+1}^{p}[/math], где степень [math]p[/math] меняется в зависимости от метода и соответствует порядку схемы.

Наибольшее распространение в вычислительной практике нашел метод Рунге-Кутты четвертого порядка:

где [math]K_{1,i}= f_{i}= f(x_{i},widehat{y}_{i}),~~ K_{2,i}= f! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2}K_{1,i}right)!,[/math]
[math]K_{3,i}= f! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2} K_{2,i}right)!,~~ K_{4,i}= f bigl(x_{i}+ h_{i+1},, widehat{y}_{i}+ h_[i+1}cdot K_{3,i}bigr).[/math]

Схема (6.41) является четырехчленной, первый коэффициент [math]K_{1,i}[/math] относится к точке [math]x_{i}[/math], второй и третий — к средней точке [math]x_{i}+ frac{h_{i+1}}{2}[/math] четвертый — к точке [math]x_{i+1}[/math]. Для этой схемы выбираются следующие параметры, входящие в (6.39):

[math]s=4,quad c_2=c_3=frac{1}{2};quad c_4=1;quad a_{2,1}= frac{1}{2};quad a_{3,1}= a_{4,1}= a_{4,2}=0;quad a_{3,2}=frac{1}{2};quad a_{4,3}=1.[/math]

Приведем геометрическую интерпретацию метода Рунге-Кутты четвертого порядка (рис. 6.13). Пусть известна точка [math](x_{i},widehat{y}_{i})[/math] на искомой интегральной кривой [math]y=y(x)[/math]. Прямая [math]AE[/math], проходящая через точку [math](x_{i}, widehat{y}_{i})[/math] и определяющая величину [math]widehat{y}_{i+1}[/math], имеет угловой коэффициент. Этот коэффициент рассчитывается по угловым коэффициентам касательных в точках [math]A,B,C,D[/math], проведенным к проходящим через эти точки интегральным кривым (штриховые линии). При этом коэффициент [math]K_{1,i}[/math] рассчитывается в точке [math]A[/math], коэффициент [math]K_{2,i}[/math] — в точке [math]C[/math], коэффициент [math]K_{3,i}[/math] — в точке [math]B[/math], а коэффициент [math]K_{4,i}[/math] — в точке [math]D[/math].

Замечания

1. Схемы Рунге-Кутты легко обобщаются на системы дифференциальных уравнений. Так, применительно к системе двух уравнений

[math]begin{cases}y'(x)= f(x,y,z),& y(x_0)=y_0,\ z'(x)= g(x,y,z),& z(x_0)=z_0 end{cases}[/math]

схема Рунге-Кутты четвертого порядка имеет вид [math](i=overline{0,n-1})[/math]

где

[math]begin{aligned}& K_{1,i}= f(x_{i}, widehat{y}_{i}, widehat{z}_{i}); &quad & K_{2,i}= f! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2} K_{1,i},, widehat{z}_{i}+ frac{h_{i+1}}{2} q_{1,i}right)!;\[2pt] & qK_{1,i}= g(x_{i}, widehat{y}_{i}, widehat{z}_{i}); &quad & q_{2,i}= g! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2} K_{1,i},, widehat{z}_{i}+ frac{h_{i+1}}{2} q_{1,i}right)!;end{aligned}[/math]

[math]begin{aligned}& K_{3,i}= f! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2} K_{2,i},, widehat{z}_{i}+ frac{h_{i+1}}{2} q_{2,i}right)!;\[2pt] & q_{3,i}= g! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2} K_{2,i},, widehat{z}_{i}+ frac{h_{i+1}}{2} q_{2,i}right)!;\[2pt] & K_{4,i}= f bigl(x_{i}+h_{i+1},, widehat{y}_{i}+ h_{i+1}K_{3,i},, widehat{z}_{i}+ h_{i+1} q_{3,i}bigr);\[2pt] & q_{4,i}= g bigl(x_{i}+h_{i+1},, widehat{y}_{i}+ h_{i+1}K_{3,i},, widehat{z}_{i}+ h_{i+1} q_{3,i}bigr). end{aligned}[/math]

2. Соотношения модифицированного метода Эйлера (6.28) можно записать в форме схемы Рунге-Кутты:

[math]K_{1,i}= f(x,widehat{y}_{i}),quad K_{2,i}= f! left(x_{i}+ frac{h_{i+1}}{2},, widehat{y}_{i}+ frac{h_{i+1}}{2}K_{1,i}right)!.[/math]

Покажем, что эта схема второго порядка. Разложение (6.38) с учетом равенства [math]y»= f_x+ f_{i}cdot f_{y}[/math], где [math]f_{i}= f(x_{i},y_{i})[/math], частные производные [math]f_x, f_y[/math] и вторая производная [math]y»[/math] вычисляются в точке [math](x_{i}, y_{i})[/math], принимает вид

Разложение функции [math]f(x,y)[/math] относительно точки [math](x_{i},y_{i})[/math] можно записать следующим образом:

где частные производные [math]f_x,f_y[/math] вычисляются в точке [math](x_{i},y_{i})[/math]. Тогда при [math]x=x_{i}+frac{h_{i+1}}{2},~ y=y_{i}+frac{h_{i+1}}{2}f_{i}[/math] получаем

[math]f! left(x_{i}+frac{h_{i+1}}{2},, y_{i}+frac{h_{i+1}}{2}f_{i}right)= f_{i}+ frac{h_{i+1}}{2}f_x+ frac{h_{i+1}}{2}f_{i}cdot f_y+ O(h_{i+1}^2).[/math]

Следовательно, формулу (6.43) можно переписать в виде

[math]widehat{y}_{i+1}= widehat{y}_{i}+ h_{i+1} f! left(x_{i}+frac{h_{i+1}}{2},, y_{i}+frac{h_{i+1}}{2}f(x_{i}, widehat{y}_{i})right)= widehat{y}_{i}+ h_{i+1}f_{i}+ frac{h_{i+1}^2}{2}(f_x+ f_{i}cdot f_y)+ O(h_{i+1}^3).[/math]

Сравнивая последнее соотношение с (6.44), можно сделать вывод о том, что схема модифицированного метода Эйлера согласуется с разложением по формуле Тейлора вплоть до членов степени [math]h_{i+1}^2[/math] и поэтому является методом Рунге-Кутты второго порядка.

3. Коэффициенты семейства методов Рунге-Кутты (6.39),(6.40) удобно записывать в виде табл. 6.7. Табл. 6.8 соответствует классическому методу Рунге-Кутты четвертого порядка (6.41), а в табл. 6.9 приведены коэффициенты другого метода четвертого порядка (правило [math]3!!not{phantom{|}},8[/math]). Табл. 6.10 соответствует методу Рунге–Кутты (6.43) второго порядка, табл. 6.11 и 6.12 — методам третьего порядка. В табл. 6.13 даны коэффициенты метода Фельберга четвертого порядка, где формула для [math]widetilde{y}_{i+1}[/math] имеет пятый порядок, а в табл. 6.14 — семистадийного метода Батчера шестого порядка. В табл. 6.15 приведены коэффициенты метода Рунге-Кутты-Мерсона четвертого порядка, где формула для [math]widetilde{y}_{i+1}[/math] имеет пятый порядок для систем линейных дифференциальных уравнений с постоянными коэффициентами и третий порядок для нелинейных систем. Параметры метода Рунге-Кутты-Дормана-Принса пятого порядка приведены в табл. 6.16. Метод Фвльберга седьмого порядка и метод Дормана-Принса восьмого порядка изложены.

4. Для s-стадийных методов Рунге-Кутты справедливо следующее условие устойчивости, определяемое для задачи (6.12):

[math]left|1+ hmu+ frac{(hmu)^2}{2}+ frac{(hmu)^3}{6}+ frac{(hmu)^4}{24}+ ldots+frac{(hmu)^s}{s!}right|<1.[/math]

Если [math]mu[/math] — действительное число, можно указать интервалы устойчивости: [math]hmuin (-2;0)[/math] при [math]s=1,!2;~ hmuin(-2,!51;0)[/math] при [math]s=3;~ hmuin (-2,!78;0)[/math] при [math]s=4[/math] и т.д.

Пример 6.6. Для задачи Коши [math]y’=x+y,~ y(0)=1,~ xin[0;0,!2][/math] выполнить два шага методом Рунге-Кутты второго порядка (6.43) и один (первый шаг) методом четвертого порядка (6.41). Шаг интефирования принять постоянным и равным 0,1. Сравнить оба численных решения с точным [math]y(x)= 2e^x-x-1[/math].

Решение

В заключение данного раздела отметим преимущества и недостатки одношаговых и многошаговых схем:

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

2. В многошаговых схемах, напротив, учитываются значения искомой функции, рассчитанные на предыдущих шагах интегрирования, что повышает порядок точности этих схем без использования дополнительных точек на отрезке [math][x_{i},x_{i+1}][/math]. Однако эти схемы требуют усложнения расчетных алгоритмов при реализации переменного шага.

Замечание. В явных одношаговых методах обеспечить заданную точность на каждом очередном шаге из текущей точки [math](x_{i}, widehat{y}_{i})[/math] можно по следующей методике.

1. Задать начальную величину шага [math]h^{(0)}= x_{i+1}^{(0)}-x_{i}[/math] и рассчитать значение [math]widehat{y}_{i+1}^{(0)}[/math] в точке [math]x_{i+1}^{(0)}[/math].

2. Вычислить величину [math]widehat{y}_{i+1}^{(1)}[/math] в точке [math]x_{i+1}^{(0)}[/math] с половинным значением шага [math]h^{(1)}= frac{1}{2}h^{(0)}[/math] (до достижения точки [math]x_{i+1}^{(0)}[/math] выполняются два шага).

3. Если [math]bigl|widehat{y}_{i+1}^{(1)}-widehat{y}_{i+1}^{(0)}bigr|< varepsilon[/math], то перейти к нахождению приближенного решения в точке [math]x_{i+2}^{(0)}= x_{i+1}^{(0)}+ h^{(0)}[/math]. Иначе осуществить еще одно дробление шага, т.е. вычислить [math]h^{(2)}= frac{1}{2}h^{(1)}[/math], и найти [math]widehat{y}_{i+1}^{(2)}[/math] в точке [math]x_{i+1}^{(1)}= x_{i}+ h^{(1)}[/math].

4. Сравнить решения, полученные с шагом [math]h^{(1)}[/math] и [math]h^{(2)}[/math] в точке [math]x_{i+1}^{(1)}[/math]. Если желаемая точность достигнута, расчеты продолжить с шагом [math]h^{(1)}[/math], а если нет, то процесс дробления шага продолжить до достижения заданной точности.

---

Принцип аналогий

Данный метод следует из обобщения интегрально-интерполяционного принципа. Он основан на применении интерполяционных интегрально-дифференциальных сплайн-функций, преобразования которых — параметрические соотношения, аппроксимационные формулы для производных и квадратурные формулы — используются для получения разностных схем решения задачи Коши. Приведем некоторые из этих соотношений, соответствующих задаче аппроксимации некоторой произвольной функции [math]Phi(x)[/math] многочленами второй степени:

а) параметрическое соотношение, связывающее интерполируемую функцию [math]Phi(x_{i})[/math] и ее производную [math]overline{m}_{i}= Pho'(x_{i})colon[/math]

б) формула аппроксимации первой производной в точке [math]x_{i}[/math] на трехточечном нерегулярном шаблоне [math](x_{i-1}, x_{i},x_{i+1})colon[/math]

где [math]H_{si}^{p(i+1)}= scdot h_{i}+ pcdot h_{i+1};~ s,,p[/math] — целые числа [math]left(H_{si}^{p(i+1)}= h_{i}(s+ pcdot delta_{i+1}),~ delta_{i+1}= frac{h_{i+1}}{h_{i}}right)[/math];

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

Эта формула при [math]h_{i+1}= h=text{const}[/math] упрощается:

[math]I_{i}^{i+1}= frac{h}{12} (5Phi_{i-2}-16Phi_{i-1}+ 23Phi_{i})quad left(frac{3}{8} h^4M_{3,i}right)!;[/math]

г) двухинтервальная (трехточечная) обобщенная на нерегулярный шаблон квадратурная формула парабол четвертого порядка на отрезке [math][x_{i-1},x_{i+1}][/math], выражающая сумму «удельных» значений интегралов через линейную комбинацию значений [math]Phi_{i-1}, Phi_{i}, Phi_{i+1}colon[/math]

д) одноинтервальная (трехточечная) функциональная квадратурная формула четвертого порядка на отрезке [math][x_{i-1},x_{i+1}][/math]

которая при [math]h_{i+1}= h=text{const}[/math] преобразуется к более простому виду:

[math]I_{i}^{i+1}= frac{h}{12} bigl(-Phi_{i-1}+ 8Phi_{i}+ 5Phi_{i+1}bigr).[/math]

Анализ остаточного слагаемого, получающегося при оценке погрешности квадратурной формулы (6.50) (см. мажоранту в скобках справа от данной формулы), указывает на то, что при сгущающейся вправо сетке (условная аппроксимация интеграла) порядок точности формулы может быть повышен на единицу без увеличения количества точек шаблона. Действительно, при [math]delta_{i+1}< h_{i}[/math] легко получить приближенное условие повышения порядка: [math]h_{i+1} leqslant h_{i}^2[/math]. Таким образом, в данном случае порядок повышается за счет значительного уменьшения шага, такого, что [math]h_{i+1}[/math] не превышает квадрата шага [math]h_{i}[/math]. Этот факт может быть использован при решении задачи Коши на нерегулярном шаблоне;

е) две одноинтервальные (двухточечные) квадратурные формулы функционально-дифференциального типа:

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

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

Так, заменяя в (6.46) [math]DeltaPhi_{k}[/math] на [math]Delta y_{k}[/math], а [math]overline{m}_{k}[/math] на [math]f_{k}~ (k=i,i-1)[/math] и разрешая относительно [math]widehat{y}_{i+1}[/math], получаем двухшаговую явную схему второго порядка (2Я2А):

которая при [math]h_{i+1}= h= text{const}[/math] преобразуется к виду [math]widehat{y}_{i+1}= 2widehat{y}_{i}-widehat{y}_{i-1}+ h(f_{i}-f_{i-1})[/math].

Комбинация схем (6.54) и (6.24) дает обобщенную на нерегулярный шаблон схему Адамса—Бэшфорта (2Я2В):

которая при [math]h_{i+1}= h= text{const}[/math] преобразуется к виду (6.32): [math]widehat{y}_{i+1}= widehat{y}_{i}+ frac{h}{2}(3f_{i}-f_{i-1})[/math].

Аппроксимационная формула (6.47) определяет неявную двухшаговую схему второго порядка (2НЯ2):

которая при [math]h_{i+1}= h= text{const}[/math] принимает вид

Замечание. Схема (6.57) является одним из методов Гира, получаемых на основе формул дифференцирования назад, связывающих значение производной в точке [math]x_{i+1}[/math] со значениями функции в предыдущих точках [math]x_{i},ldots, x_{i-k+1}[/math] (см. рис.6.2). Схема (6.57) является схемой второго порядка, а неявный метод Эйлера — схемой первого порядка. Приведем более точные многошаговые схемы, записанные на регулярных шаблонах:

– схема третьего порядка:

– схема четвертого порядка:

[math]widehat{y}_{i+1}= frac{48}{25} widehat{y}_{i}-frac{16}{25} widehat{y}_{i-1}+ frac{16}{25} widehat{y}_{i-2}-frac{3}{25}widehat{y}_{i-3}+ frac{12}{25}h, f(x_{i+1}, widehat{y}_{i+1});[/math]

– схема пятого порядка:

[math]widehat{y}_{i+1}= frac{300}{137} widehat{y}_{i}-frac{300}{137} widehat{y}_{i-1}+ frac{200}{137} widehat{y}_{i-2}-frac{75}{137}widehat{y}_{i-3}+ frac{12}{137} widehat{y}_{i-4}+ frac{60}{137}h, f(x_{i+1}, widehat{y}_{i+1});[/math]

– схема шестого порядка:

[math]widehat{y}_{i+1}= frac{360}{147} widehat{y}_{i}-frac{450}{147} widehat{y}_{i-1}+ frac{400}{147} widehat{y}_{i-2}-frac{225}{147}widehat{y}_{i-3}+ frac{17}{147} widehat{y}_{i-4}-frac{10}{147} widehat{y}_{i-5}+ frac{60}{147}h, f(x_{i+1}, widehat{y}_{i+1}).[/math]

Известно, что схема (6.58) является A(α)-устойчивой с [math]alphacong 68^{circ}[/math].

Все рассматриваемые ниже схемы получаются из квадратурных формул (6.48)-(6.52), в которых вместо разности первообразных [math]F_{i+1}-F_{i}= I_{i}^{i+1}[/math] подставляется [math]Delta widehat{y}_{i}[/math], а вместо [math]Phi_{k}(k=i-2,i-1,i)[/math] и [math]Phi’_{k}[/math] подставляются [math]f_{k}[/math] и [math]f’_{k}[/math]. При этом порядки точности схем на отрезке [math][a,b][/math] соответствуют порядкам точности вычисления интегралов на всем отрезке интегрирования.

Преобразуя таким образом квадратурную формулу (6.48), получаем явную трехшаговую схему третьего порядка (ЗЯЗ):

При [math]h_{i+1}= h= text{const}[/math] она является традиционной, получающейся интегрально-интерполяционным методом (6.33):

[math]widehat{y}_{i+1}= widehat{y}_{i}+ frac{h}{12}bigl(5f_{i-2}-16 f_{i-1}+ 23f_{i}bigr).[/math]

Из обобщенной формулы (6.49) следует первая двухшаговая неявная схема третьего порядка (2НЯЗА), которая при [math]h_{i+1}= h= text{const}[/math] преобразуется к виду (6.29)

Одноинтервальная (трехточечная) квадратурная формула (6.50) определяет вторую двухшаговую неявную схему третьего порядка (2НЯЗБ)

Эта схема в соответствии с вышеприведенной оценкой точности интеграла при [math]delta_{i+1} leqslant h_{i}[/math] имеет не третий, а четвертый порядок точности. При [math]h=text{const}[/math] эта схема преобразуется к виду (6.36):

[math]widehat{y}_{i+1}= widehat{y}_{i}+ frac{h}{12} bigl(-f_{i-1}+ 8f_{i}+ 5f(x_{i+1}, widehat{y}_{i+1})bigr).[/math]

И наконец, последние приводимые здесь схемы следуют из квадратурных формул (6.51), (6.52) (неявные одношаговые схемы) и параметрического соотношения (6.53) (неявная двухшаговая схема):

[math]begin{aligned}widehat{y}_{i+1}&= widehat{y}_{i}+ frac{h_{i+1}}{3} bigl[f_{i}+ 2f(x_{i+1}, widehat{y}_{i+1})bigr]-frac{h_{i+1}^2}{6}f'(x_{i+1}, widehat{y}_{i+1});\[2pt] widehat{y}_{i+1}&= widehat{y}_{i}+ h_{i+1}f_{i}+ frac{h_{i+1}^2}{6} bigl[2f’_{i}+ f'(x_{i+1}, widehat{y}_{i+1})bigr];\[2pt] widehat{y}_{i+1}&= widehat{y}_{i}+ delta_{i+1} Delta widehat{y}_{i+1}+ frac{h_{i+1}}{6} bigl[h_{i} f’_{i-1}+ 2H_{i}^{i+1}f’_{i}+ h_{i+1} f'(x_{i+1}, widehat{y}_{i+1})bigr], end{aligned}[/math]

где [math]f’=f_{x}+f_{y}cdot f[/math]. Эти схемы обозначаются 1НЯЗА, 1НЯЗБ, 2НЯЗД соответственно.

При [math]h_{i+1}= h= text{const}[/math] последняя схема преобразуется к виду

[math]widehat{y}_{i+1}=-widehat{y}_{i-1}+ 2widehat{y}_{i}+ frac{h^2}{6} bigl[f’_{i-1}+ 4f’_{i}+ f'(x_{i+1}, widehat{y}_{i+1})bigr].[/math]

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

---

Методика решения задачи Коши методами сеток

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

2. Задать сетку [math]Omega_{n}[/math] или выбрать ее параметр ([math]h[/math], если сетка равномерная; [math]h_{i+1},~ i=overline{0,n-1}[/math], если сетка неравномерная). Если выбранный метод является ограниченно устойчивым, то для задания параметра сетки оценить величину критического шага.

3. Если выбранный метод многошаговый, найти требуемое число «разгонных» точек.

4. Найти приближенное решение задачи Коши во всех точках сетки [math]Omega_{n}[/math]. Если выбран явный метод, каждая следующая точка [math]widehat{y}_{i+1}[/math] рассчитывается по явной формуле, а если неявный, то с помощью решения в общем случае нелинейного алгебраического уравнения (6.11) одним из численных методов.

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

Математический форум (помощь с решением задач, обсуждение вопросов по математике).

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