Применение интерполяции для решения уравнений

Интерполяция применяется для решения уравнений вида

(2.17)

Если в области корня уравнения (2.17) вычислить его левую часть n +1 точке и результаты поместить в табл. 2.1, то для определения корня можно поменять местами столбцы таблицы и с помощью одного из алгоритмов интерполяции найти значение аргумента х, при котором функция f(x) принимает значение . Нахождение значений аргумента x по заданным значениям функции называется обратной интерполяци­ей.

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

где х^к — приближение к корню на k-й итерации; — заданная по­грешность определения корня; — предельная величина невязки.

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

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

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

Введем обозначение , тогда и квадратное уравнение принимает вид

где

-

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

На следующей итерации полином строим по точкам Окончание итерационного процесса наступит, когда выполнится одно из условий

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

 

Интерполяция сплайнами

Полиномиальная интерполяция не всегда дает удовлетворитель­ные результаты при аппроксимации зависимостей. Так, например, при представлении полиномами резонансных кривых колебательных си­стем большая погрешность возникает на концах ("крыльях") этих кри­вых. Несмотря на выполнение условий Лагранжа в узлах, интерполя­ционная функция может иметь значительное отклонение от аппрокси­мируемой кривой между узлами. При этом повышение степени интер­поляционного полинома приводит не к уменьшению, а к увеличению погрешности. Возникает так называемое явление волнистости [4].

Для проведения гладких кривых через узловые значения функции чертежники используют упругую металлическую линейку, совмещая ее с узловыми точками. Математическая теория подобной аппроксима­ции развита за последние тридцать лет [5, 6] и называется теорией сплайн-функций. Разработано и обширное программное обеспечение для применения сплайнов в различных областях науки и техники [7, 1].

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

(2.18)

Функцию и будем использовать для аппроксимации зависимости f(x), заданной в узлах ,...,х_п значениями .

Если в качестве функции выбрать полином, то в соответствии с уравнением (2.18) степень полинома должна быть не выше третьей. Этот полином называют кубическим сплайном, который на интервале записывают в виде

(2.19)

где и — коэффициенты сплайна, определяемые из дополни­тельных условий; i = 1,2,..., п — номер сплайна.

В отличие от полиномиальной интерполяции, когда вся аппрокси­мируемая зависимость описывается одним полиномом, при сплайно-вой интерполяции на каждом интервале строится отдельный полином третьей степени(2.19) со своими коэффициентами.

Коэффициенты сплайнов определяются из условий сшивания со­седних сплайнов в узловых точках:

1) Равенство значений сплайнов и аппроксимируемой функ­ции f(x) в узлах — условия Лагранжа

(2.20)

2) Непрерывность первой и второй производных от сплайнов в уз­лах

(2.21)

(2.22)

Кроме перечисленных условий необходимо задать условия на кон­цах, т.е. в точках и. В общем случае эти условия зависят от кон­кретной задачи. Часто используют условия свободных концов сплай­нов. Если линейка не закреплена в точках вне интервала [ , ], то там она описывается уравнением прямой, т.е. полиномом первой степени. Следовательно, исходя из условий (2.22), т.е. непрерывности вторых производных сплайнов на концах интервала, запишем соотношения:

(2.23)

- (2.24)

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

Получим алгоритм определения коэффициентов кубических сплай­нов из условий (2.20 — 2.24). Условия (2.20) в узлах x_i-1 и x _i после подстановки i -гo сплайна принимают вид

a_i = ƒ_{i-1 ,} (2.25)

a_i + b_i h_i + c_ih_i² + d_ih_i³ = ƒ_{i ,}(2.26)

где h_i = x_i — x_{i- 1}, 1 ≤ i ≤ n.

Продифференцируем дважды сплайн (2.19) по переменной x:

φ′_i(x) = b_i + 2c_i (x— x_i-1) + 3d_i (x— x_i-1)², (2.27)

φ″_i(x) = 2c_i + 6d_i (x— x_i-1). (2.28)

Из условий непрерывности производных (2.21) и (2.22) при пере­ходе в точке x_iот i-го к (i₁)сплайну с учетом выражений (2.27) и (2.28) получим следующие соотношения:

b_i + 2c_i h_i + 3d_ih_i² = b_{i+1 ,}(2.29)

c_i + 3d_ih_i = c_i + 1. (2.30)

И, наконец, из граничных условий (2.23) и (2.24) на основании вы­ражения для второй производной (2.28) получим, что

c₁ = 0, (2.31)

c_n + 3d_n h_n = 0. (2.32)

Соотношения (2.25), (2.26), (2.29) — (2.32) представляют собой полную систему линейных алгебраических уравнений относительно ко­эффициентов сплайнов a_i, b_i, c_i и d_i. Но прежде чем решать эту систе­му, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов c_i.

Из уравнений (2.30) коэффициенты d_iвыразим через коэффициен­ты c_i:

c_i+1 —c _i

d_i = ————. (2.33)

3 h_i

Объединяя уравнения (2.25) и (2.26) с уравнением (2.33), предста­вим коэффициенты b_iтакже через коэффициенты c_i:

f _i — f _i-1(c _i+1 + 2c _i)h

b_i = ———— - —————. (2.34)

h_i 3

После подстановки выражений (2.33) и (2.34) в соотношение (2.29) получим уравнение, в которое входят только неизвестные коэффици­енты c_i. Для симметричности записи в полученном уравнении умень­шим значение индекса i на единицу:

h_i-1 c_i-1+2(h_i-1+h_i)c_i + h_ic_i+1 = 3 (2.35)

где 2≤ i ≤ n.

При i = n, учитывая условие свободного конца сплайна, в уравне­ний (2.35) следует положить

c_n+1 = 0. (2.36)

Таким образом, п — 1 уравнение вида (2.35) вместе с условиями (2.31) и (2.36) образует систему линейных алгебраических уравнений для определения коэффициентов c_i. Коэффициенты d _i и b _i вычисляют­ся после нахождения c_i по формулам (2.33) и (2.34), коэффициенты a_i равны значениям аппроксимируемой функции в узлах в соответствии с формулой (2.25).

В каждое из уравнений типа (2.35) входит только три неизвестных с последовательными значениями индексов c_i-1, c_i, c_i+1. Следователь­но, матрица системы линейных алгебраических уравнений относитель­но c_i является трехдиагональной, т.е. имеет отличные от нуля элементы только на главной и двух примыкающих к ней диагоналях. Для реше­ния систем с трехдиагональной матрицей наиболее эффективно при­менять так называемый метод прогонки, являющийся частным случаем метода Гаусса.

Рассмотрим алгоритм метода прогонки применительно к решаемой задаче. Для сокращения записи уравнение (2.35) представим в виде

h_i-1 c_i-1+s_i c_i + h_ic_i+1 = r_i, (2.37)

где s_i = 2 (h_i-1 - h_i ),

r_i = 3 . (2.38)

Так же, как и метод Гаусса, метод прогонки разделяется на два эта­па — прямой и обратный ходы. В процессе прямого хода метода про­гонки вычисляют значения коэффициентов линейной связи каждого предыдущего неизвестного c_i с последующим c_i+1.

Из уравнения (2.37) при i = 2 с учетом граничного условия (2.31) установим связь коэффициента с₂ с коэффициентом с₃:

 

c₂ = k₂ - l₂с_3, (2.39)

где k₂ и l₂ — прогоночные коэффициенты,

 

. r ₂ h₂

k₂ = —, l₂ = —.

S ₂ S₂

Затем, подставляя выражение (2.39) в уравнение (2.37) при i = 3, получим линейную связь коэффициента с₃ с коэффициентом с₄:

c₃ = k₃ - l₃с_4,

Поступая аналогичным образом для любых соседних коэффициентов с номерами i и i + 1, можно установить линейную связь в виде

c_i = k_i - l_iс_i+1 (2.40)

В процессе выполнения прямого хода метода прогонки необходимо вы­числить значения всех прогоночных коэффициентов k _i и l_i, для кото­рых получим рекуррентные соотношения. Подставим формулу связи

(i -1)-го и i -го коэффициентов

c_i-1 = k_i-1 - l_i-1с_i

в уравнение (2.37), в результате получим

r_i - h_i-1k_i-1h_i

c_i = ――――― - ————— c_i+1.

s_i - h_i-1l_i-1 s_i - h_i-1 l_i-1

Сравнение последнего соотношения с формулой (2.40) позволяет за­писать рекуррентные соотношения для определения прогоночных ко­эффициентов k_i и l_i:

r_i - h_i-1k_i-1h_i

k_i = ――――――, l_i = ——————. (2.41)

s_i - h_i-1l_i-1 s_i - h_i-1 - l_i-1

 

Учитывая граничное условие (2.31) и соотношение

c₁ = k₁ - l₁с_2,

а также полагая c₂ ≠0, задаем начальные коэффициенты k₁ = 0 и l₁ = 0. Затем по формуле (2.41) вычислим все п пар прогоночных ко­эффициентов k _i и l_i..

На основании соотношения c_n = k_n - l _nc_n+1 и граничного условия (2.36) получим, что

c_n = k_n. (2.42)

далее последовательно применим формулу (2.40) при i = n — 1, n — 2,...., 2 и вычислим значения искомых величин c_n-1, c_n-2,…, c_{2 .}

Эта процедура называется обратным ходом метода прогонки.

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

После определения всех коэффициентов с_i другие коэффициенты сплайна вычисляются по формулам (2.25), (2.33) и (2.34). Аппрокси­мирующую функцию φ(x) можно рассчитать с помощью соотношения (2.19) в любой точке x на интервале [ х₀, х_n ].

Глава 3

Метод наименьших квадратов

Общие положения

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

Обозначим узлы исходной таблицы данных через x _j, где 0 ≤ j ≤ n — номер узла. Считаем известными значения эксперимен­тальных данных в узловых точках ƒ(x_j) = ƒ_j. Введем непрерывную функцию φ (х) для аппроксимации дискретной зависимости f(x _j). В узлах функции φ (х) и f (x)будут отличаться на величину ε _j = φ (x _j) — - f(x_j). Отклонения ε _j могут принимать положительные и отрица­тельные значения. Чтобы не учитывать знаки, возведем каждое откло­нение в квадрат и просуммируем квадраты отклонений по всем узлам:

n n

Q = ∑ ε _j ² = ∑[φ(x _j) - ƒ(x _j)]² (3.1)

j=0 j=0

Метод построения аппроксимирующей функции φ (х)из условия минимума величины Qполучил название метода наименьших квадра­тов (МНК).

Наиболее распространен способ выбора функции φ (х)в виде ли­нейной комбинации

φ (х) = c₀ φ ₀ (x)+ c₁ φ ₁ (x) +... + c_m φ _m(x), (3.2)

φ ₀ (х), φ ₁ (x),..., φ _m (х) — базисные функции; m ≤ n;

c₀(x), c₁ (x),.. c_m (x)- коэффициенты, определяемые при минимизации величины Q.

Математически условия минимума квадратов отклонений Q запи­шем, приравнивая к нулю частные производные от Qпо коэффициентам 0 ≤ k ≤ m:

n

∂ Q⁄∂ c₀ = 2∑ [c₀ φ ₀ (x _j) + c₁ φ ₁(x _j) +…+

j=0

+ c_m φ _m (x _j) - ƒ_j] φ ₀ (x _j) = 0,

n

∂ Q⁄∂ c₁ = 2∑ [c₀ φ ₀ (x_j)+c₁φ₁(x_j) +…+ (3.3 )

j=0

+ c_m φ _m (x _j) - ƒ_j] φ ₁ (x _j) = 0,

…………………………………….

n

∂ Q⁄∂ c_m = 2∑ [c₀ φ ₀ (x _j) + c₁ φ ₁(x _j) +…+

j=0

+ c_m φ _m (x _j) - ƒ_j] φ _m (x _j) = 0,

Из системы линейных алгебраических уравнений (3.3) определя­ются все коэффициенты с‌‍_k. Система (3.3) называется системой нор­мальных уравнений. Матрица этой системы имеет следующий вид:

(3.4)

и называется матрицей Грамма. Элементы матрицы Грамма являются скалярными произведениями базисных функций:

n

(φ _j, φ _k) = ∑ φ _j (x _j) φ _k(x _j). (3.5)

j=0

 

Расширенная матрица системы уравнений (3.5) получится добав­лением справа к матрице Грамма столбца свободных членов

(3.6)

,

где скалярные произведения, являющиеся элементами столбца, опре­деляются аналогично (3.5):

(3.7)

n

(φ _j, ƒ) = ∑ φ _j (x _j) ƒ _j

j=0

Отметим основные особенности матрицы Грамма, полезные при про­граммной реализации алгоритмов МНК:

1) Матрица симметрична, т.е. a‍‍ _{i‍‍‍‍‍‍ j} = a‍‍_{‍‍‍ j i}, что позволяет сократить объем вычислений при заполнении матрицы.

2) Матрица является положительно определенной, следовательно, при решении системы нормальных уравнений методом исключе­ния Гаусса можно отказаться от процедуры выбора главного эле­мента.

3) Определитель матрицы будет отличен от нуля, если в качестве базиса выбраны линейно независимые функции φ‌_k (x), при этом система (3.3) имеет единственное решение

При обработке экспериментальных данных, определенных с погреш­ностью е в каждой узловой точке, обычно начинают с аппроксима­ции функцией φ‌_k (x), представимой одной-двумя базисными функци­ями. После определения коэффициентов c‍_kвычисляют величину Q по формуле(3.1). Если получится, что √Q > e, то необходимо расширить базис добавлением новых функций φ‌_k (x).Расширение базиса необхо­димо осуществлять до тех пор, пока не выполнится условие √Q ≈ e.

Выбор конкретных базисных функций зависит от свойств аппрок­симируемой функции ƒ(x), таких как периодичность, экспоненциаль­ный или логарифмический характер, свойства симметрии, наличие ас-симптотики и т.д.[8]

Степенной базис

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

φ‌₀ (x) = х ⁰ = 1, φ‌₁ (x) = х¹ = х,..., φ‌_m (х) = х^m. (3.8)

В этом случае также, как и при интерполяции, мы будем аппрок­симировать экспериментальную зависимость полиномом. Однако сте­пень полинома m выбираем обычно m << n(при лагранжевой интер­поляции m = n). Аппроксимирующая кривая в МНК не проходит че­рез значения исходной функции в узлах, но проведена из условия наи­меньшего суммарного квадратичного отклонения. Экспериментальные данные "сглаживаются" с помощью функции φ‌(x). Если же выбрать m = n, то на основании единственности интерполяционного полинома получим функцию φ(x), совпадающую с каноническим интерполяци­онным полиномом степени n, аппроксимирующая кривая пройдет че­рез все экспериментальные точки и величина Qбудет равна нулю. По­следнее обстоятельство используется для отладки и тестирования про­грамм, использующих МНК.

Запишем расширенную матрицу системы нормальных уравнений для базиса (3.8):

(3.9)

Нетрудно видеть, что для формирования расширенной матрицы (3.9) достаточно вычислить элементы первой строки и двух последних столб­цов, остальные элементы не являются "оригинальными" и заполняются с помощью циклического заполнения.

Для решения систем уравнений с матрицей Грамма разработаны методы сингулярного разложения [9]. Если же m ≤ 4 ÷ 5, то такие системы можно решать и более простым методом исключения Гаусса.