Приближение функций

ПРИБЛИЖЕНИЕ ФУНКЦИЙ 1 Основные понятия Приближение или аппроксимация функций – это приближенная замена одних функций другими. Приближение данной функции f  x  - это нахождение функции g  x  из некоторого определенного класса, в том или ином смысле близкой к функции f  x  , дающей ее приближенное представление. Существует много различных вариантов задачи о приближении функций в зависимости от того, какие функции приближающая используются функция для g x , приближения, как как понимается ищется близость приближаемой f  x  и приближающей g  x  функций. Многообразие методов обусловлено многообразием различных постановок задачи, продиктованных необходимостью изучения абстрактных моделей реальных проблем. В задачах приближения функций приближаемая функция обычно задается на дискретном множестве точек своей области определения. Это множество называют сеткой, а сами точки – узлами сетки. Значения приближаемой функции в узлах сетки называются ее узловыми значениями. Для функции, заданной на n  мерной вещественной сетке R , используется обозначение R : xi , f  xi , i  0,1,  , n , где xi  узлы сетки; f  xi   узловые значения приближаемой функции. 1 В качестве приближающих функций наиболее часто используются так называемые обобщенные полиномы Q  x  вида g  x   Q x   c00  x   c11  x     cm m  x . 0  x , 1 x , m  x   где базисная система функций приближающего обобщенного полинома; c1 , c2 , , cm  произвольные коэффициенты.   В простейшем случае, если базисная система функций  j  x  представляет собой последовательность целых неотрицательных степеней переменной x , то приближающая функция представляет собой обычный полином Q x   c0  c1x  c2 x 2    cm x m . В зависимости от характера критерия близости приближаемой f  x  и приближающей g  x  функций в проблеме приближения функций можно выделить два вида задач, имеющих наибольшую важность для приложений: 1) интерполирование (интерполяция) функций, когда требуется, чтобы в определенных точках (узлах интерполяции) совпадали значения функций f  x  и g  x ; 2) аппроксимирование (аппроксимация) функций, когда близость функций f  x  и g  x  оценивается по среднеквадратичному или максимальному абсолютному отклонению функции g  x  от функции f  x  в определенных точках (узлах аппроксимации). 2 Узлами интерполяции или аппроксимации в задачах приближения функций служат обычно узлы сетки, на которой задается приближаемая функция. 2 Интерполирование функций 2.1 Общая схема интерполяции функций Рассмотрим функцию f  x  , заданную на сетке R : xi , f  xi , i  0,1,  , n. , (2.1) В общей задаче интерполяции для функции f  x  отыскивается приближающая функция g  x  в виде обобщенного полинома g  x   Q x   c00  x   c11  x     cm m  x . (2.2) Условие интерполяции g  xi   f  xi , i  0,1,, n. (2.3) эквивалентно системе линейных уравнений порядка n  1  m  1 относительно произвольных коэффициентов c j , j  0,1,, m c00  xi   c11  xi     cm m  xi   f  xi , i  0,1,2, n. (2.4) Задача интерполяции состоит в нахождении коэффициентов c j , обеспечивающих выполнение равенств (2.3). Эта задача имеет однозначное решение, если выполняются условия: 3 0  x0  1  x0    n  x0  0  x1  1  x1    n  x1  m  n;    0,  0  xn  1  xn    n  xn  (2.5) где   определитель основной матрицы системы уравнений (2.4). Условие m  n обязательно, а неравенство   0 , очевидно, выполняется, если базисные функции линейно независимы и приближаемая функция задана на сетке без кратных узлов. 2.2 Интерполяция функций обычными Интерполяционный полином Лагранжа   Если базисная система функций  j  x  представляет собой совокупность степенных функций с основанием x неотрицательными показателями степени, полиномами. то и целыми соответствующая приближающая функция Q x   c0  c1x  c2 x 2    cm x m (2.6) называется интерполяционным полиномом. Если m  n , то система уравнений (2.4) для определения коэффициентов интерполяции в рассматриваемом случае будет иметь вид c0  c1x0  c2 x02    cn x0n  y0 ; c0  c1x1  c2 x12    cn x1n  y1;  c0  c1xn  c2 xn2    cn xnn  yn , 4 (2.7) где yi  f  xi   узловые значения приближаемой функции, заданной на сетке (2.1). Определитель  основной матрицы системы уравнений (2.7) представляет собой так называемый определитель Вандермонда 1 x0 x02  x0n  1 x1 x12  x1n   1 xn xn2  xnn  xq  x p . (2.8) 0 p  q  n Если приближаемая функция задана на сетке без кратных узлов, то определитель (2.8) не равен нулю и, следовательно, система (2.7) имеет единственное решение. Полином Ln  x   Q  x  , коэффициенты которого определяются из системы (2.7), называется интерполяционным полиномом Лагранжа. Этот полином можно представить в явном виде n Ln  x    n  i 0 j 0 ( j i ) x  x j  y  xi  x j  i  x  x0  x  xi 1   x  xi 1  x  xn  y . i i  0  xi  x0   xi  xi 1   xi  xi 1   xi  xn  n   5 (2.9) Рассмотрим пример построения интерполяционного полинома для функции, заданной на сетке i 1 2 3 xi -2,0 -1.0 1,0 yi 0,3 - 0,6 0,2 -1,2 По формуле (2.9) сформируем для заданной функции интерполяционный полином Лагранжа L3  x  (рис.1.1): L3  x   x  x j  y   x  x1  x  x2  x  x3  y    x  x  i  x  x  x  x  x  x  0 i j 1 2 3 i 0 j 0 3 3 ( j i)  x  x0  x  x2  x  x3  y   x  x0  x  x1  x  x3  y   x1  x0  x1  x2  x1  x3  1  x2  x0  x2  x1  x2  x3  2  x  x0  x  x1  x  x2  y   x  1 x  0 x  1 0,3    x3  x0  x3  x1  x3  x2  3  2  1 2  0 2  1  x  2 x  0 x  1  0,6   x  2 x  1 x  1 0,2    1  2 1  0 1  1 0  20  10  1  x  2 x  1 x  0  1,2  0,65 x3  1,1 x 2  0,35 x  0,2.  1  21  11  0  6 0.2 L(x) -0.2 -0.4 -0.6 -0.8 -1 -1.2 -2 -1.5 -1 -0.5 x 0.5 1 Рисунок 1.1 - График интерполирующего полинома Лагранжа для функции, заданной на сетке R : {[-2; 0.3],[-1;-0,6],[0;0,2],[1;-1,2]}. 7 3 Аппроксимация функций 3.1 Аппроксимация функций по методу наименьших квадратов Аппроксимация функций по методу наименьших квадратов – это стандартный метод приближения функций обобщенными полиномами, когда число m базисных функций меньше количества узлов n сетки, на которой задана приближаемая функция. Такая ситуация характеризуется как проблема избыточных данных и она типична для разнообразных задач, связанных с функциональной интерпретацией эмпирических данных. В таких задачах полагается, что измеряемые переменные связаны гипотетической функциональной зависимостью, которую требуется аналитически выразить, возможно, более простым способом. Аналогичная ситуация характерна для статистических задач, связанных с регрессионным анализом. В указанных задачах нельзя обеспечить интерполирующий характер приближающего обобщенного полинома, так как система уравнений для определения коэффициентов полинома оказывается “переопределенной”: число уравнений больше числа неизвестных. А такая задача точного решения, как правило, не имеет. Метод наименьших квадратов основан в большей степени на чистой математической логике, чем на индуктивном физическом рассуждении. Предложили этот метод независимо друг от друга Гаусс (1809) и Лежандр (1806), и вот уже двести лет он успешно используется в самых различных отраслях знания. Рассмотрим функцию f  x  , заданную на n  мерной сетке 8 xi , f  xi , i  1,, n, и приближающую функцию g  x  в виде m  компонентного обобщенного полинома Q  x  g  x   Q  x   c0 0  x   c11  x     cm m  x  при условии избыточности данных n  m. Критерий r близости функций f  x  и g  x  в методе наименьших квадратов принимается в виде суммы квадратов отклонений этих функций в узлах аппроксимации 2  m   r    c j j  xi   yi  ,   i 1  j 1  n (3.1) где yi  f  xi . r Величина представляет собой отклонение приближающей функции среднеквадратическое g x функции f  x  на множестве узлов аппроксимации Решением рассматриваемой совокупность c j , j  1,2,, m, значений задачи от приближаемой xi , i  1,2,, n. аппроксимации коэффициентов служит аппроксимации которой будет соответствовать минимальное значение критерия близости функций r 9 r c1 , c2 , , cm   min . (3.2) c j  Следовательно, задача отыскания совокупности коэффициентов c j , j  1,2,, m, наилучших в смысле принятого критерия близости функций f  x  и g  x  , будет заключаться в нахождении минимума r c1 , c2 ,, cm . функции Используя необходимые условия экстремума функции нескольких переменных, можно получить квадратную определения систему линейных значений уравнений порядка коэффициентов m для аппроксимации, удовлетворяющих условию (2.15) 2  m  r   0 c j j  xi   yi   0     ck ck i 1  j 1  n n m i 1 j 1   2 k  xi   (c j j  xi   yi )  0  (3.3) m n  n       j  xi  k  xi  c j   yi k  xi , k  1,2,, m.  j 1 i 1 i 1 Полученную систему уравнений можно представить в матричной форме WC  B, 10 (3.4)   n где W  wkj ; wkj    k  xi   j  xi  ; i 1 n B  bk  ; bk    k  xi  yi . i 1 Система уравнений (3.4) называется нормальной системой уравнений метода наименьших квадратов. В прикладных расчетах рационально используется мультипликативная факторизация основной матрицы W матричная и вектора B правых частей нормальной системы (3.4): W  ПT П ; B  ПTY , (3.5) где Y  вектор-столбец узловых значений приближаемой функции Y   y1, y2 ,, yn T , (3.6) а П  так называемая матрица планов порядка n  m , элементами i й строки которой служат значения базисных функций аппроксимирующего полинома в i  ом узле аппроксимации  1  x1   2  x1  m  x1      x    x   x   1 2 2 2 m 2  П .                  1  xn   2  xn  m  xn   11 (3.7) Рассмотрим пример построения по методу наименьших квадратов аппроксимирующего полинома второй степени для функции, заданной на сетке: i 1 2 3 4 xi -2,0 -1.0 1,0 yi 0,3 - 0,6 0,2 -1,2 Система базисных функций для аппроксимирующего полинома принимаем в виде 1  x   1;  2  x   x ; 3  x   x 2 . А сам полином определяется выражением Q x   c1  c2 x  c3 x 2 . Сформируем матрицу планов П в соответствии с выражением (3.7) 1  2  2 2  1  2 4  2  1  1 1  1  1  1  .  2 1 0 0  1 0 0   2  1 1 1 1 1 1     и вектор-столбец Y узловых значений аппроксимируемой функции 12  0,3   0,6 . Y   yi     0,2      1,2  По формулам (3.5) построим основную матрицу W и векторстолбец B правых частей нормальной системы уравнений рассматриваемой задачи: T 1  2 4 1  2 4 1  2 4 1 1 1 1  1  1 1  1  1 1   1  1 1 T         W     2 1 0 1   1 0 0 1 0 0 1 0 0  4 1 1         1 1 1  1 1 1  1 1 1   0,3  1 1 1   4 2 6  1    1,3   , 6     1,2  .   2 6  8 ; B   T Y   2  1 0 1       0,2     6  8 18   4 1 0 1    0,6  1 , 2   Расширенная матрица W нормальной системы уравнений имеет вид W  W  4  2 6  1,3  B    2 6  8  1,2  .    6  8 18  0,6 Найдем решение системы нормальных уравнений с помощью алгоритма Гаусса-Жордана: 13  4  2 6  1,3  1 W   2 6  8  1,2   0     6  8 18  0,6 0 1 0 1  0,51 1 0  0 1  1  0,37   0 1    0 0 4  0,5  0 0  0,5 1,5  0,325 5  5  1,85    5 9 1,35  0  0,385 0  0,495 ;  1  0,125  c1  0,385 ; c 2  0,495 ; c 3  0,125 . Построенный по методу наименьших квадратов аппроксимирующий полином для заданной функции будет иметь вид, рис.2.5 Q  x    0,385  0,495 x  0,125 x 2 . 0.2 Q (x) -0.2 -0.4 -0.6 -0.8 -1 -1.2 -2 -1.5 -1 -0.5 x 0.5 1 Рис. 2.5. График полинома Q  x  , аппроксимирующего по методу наименьших квадратов функцию, заданную на сетке R : {[-2; 0.3],[-1;-0,6],[0;0,2],[1;-1,2]}. 14 На рис. 2.6 для сравнения показаны совмещенные графики рассмотренных приближающих функций для функции, заданной на сетке R : {[-2; 0.3],[-1;-0,6],[0;0,2],[1;-1,2]}. SP(x) L(x), SP(x), Q(x), T(x) 0.2 Q(x) -0.2 T(x) -0.4 -0.6 L(x) -0.8 -1 -1.2 -2 -1.5 Рис. -1 -0.5 x 0.5 1 2.6. Совмещенные графики приближающих функций для функции, заданной L x   SP  x   Q x  на сетке R : {[-2; интерполирующий 0.3],[-1;-0,6],[0;0,2],[1;-1,2]}: полином интерполирующий аппроксимирующий Лагранжа; кубический по МНК сплайн; полином; T  x   аппроксимирующий полином наилучшего приближения. 15 ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ Основная задача Коши – это задача отыскания частного решения обыкновенного дифференциального уравнения порядка n   y n   f x, y , y,, y n 1 . которое удовлетворяет заданным начальным условиям: y  x0   y0 , где y x0   y0 ,, y n 1  x0   y0n 1 , x0 , y0 , y0 ,, y0n 1  фиксированные начальные значения аргумента, искомой функции и ее первых (n  1)  ой производных. Аналогичная задача может решаться для системы обыкновенных дифференциальных уравнений или уравнений с частными современной практике производными. Наиболее употребительным в приближенным методом решения задачи Коши для обыкновенных дифференциальных уравнений является метод Рунге-Кутта четвертого порядка. Исторически одним из первых и наиболее широко известных методов численного интегрирования дифференциальных уравнений является метод Эйлера для решения дифференциального уравнения первого порядка. Общий анализ показывает, что метод Эйлера по существу может рассматриваться как один из методов Рунге-Кутта. 16 Существует два широких класса приближенных методов решения задачи Коши для обыкновенных дифференциальных уравнений. Это одноступенчатые и многоступенчатые методы. В одноступенчатых методах на каждом шаге вычислительного процесса используется информация об интегральной кривой решаемого уравнения только в одной ее точке, а именно в текущей точке. Практически наиболее удобными методами этого класса являются методы Рунге-Кутта. Эти методы являются прямыми – без локальных итераций на отдельных шагах. Общий недостаток таких методов состоит в том, что при их использовании трудно оценить допускаемую ошибку. В то же время эти методы очень легко позволяют менять величину шага интегрирования. В многоступенчатых методах на каждом шаге общего вычислительного процесса осуществляются локальные итерации для обеспечения достаточной точности решения. Большинство методов этого класса называются методами прогноза и коррекции. Достоинством многоступенчатых методов является возможность осуществлять хорошие оценки точности вычислений. Общий недостаток таких методов – это необходимость временно обращаться к методам Рунге-Кутта при любом изменении величины шага интегрирования. Для сравнительной оценки точности различных методов численного интегрирования дифференциальных уравнений в качестве эталона обычно принимается разложение искомого решения в ряд Тейлора. 17 1. Метод Тейлора Этот метод, будучи теоретически пригодным для решения любых дифференциальных уравнений, с вычислительной точки зрения не представляет почти никакого практического интереса. Однако как средство для сравнительной оценки различных численных методов решения задачи Коши метод Тейлора играет универсальную роль. Рассмотрим дифференциальное уравнение первого порядка y  f  x, y  (1.1) и запишем разложение функции y  x  в ряд Тейлора в окрестности точки x  xk y  x   yk  yk  x  xk   ( j) где yk  yk  x  xk 2  yk  x  xk 3  , (1.2) 2! 3!  j -я производная функции y  x  в точке x  xk . Если для задачи Коши дифференциального уравнения (1.1) с начальным условием y  x0   y0 найдено приближенное решение для k точек x1 , x2 ,, xk с шагом h , то приближенное решение для точки xk 1  xk  h , используя формулу (1.2), можно представить в виде h2 h3 yk 1  yk  h yk  yk  yk   . 2! 3! 18 (1.3) Чем больше членов ряда (1.3) будет удерживаться для yk 1 , тем точнее будет приближение. вычисления величины Практическое неудобство формулы (3.1) состоит в необходимости вычислять производные различных порядков функции Выражение для первой производной в точке x  xk yx  . имеем из уравнения (1.1) yk  f  xk , yk . Общее выражение для второй производной можно получить, дифференцируя по x f  x, y  как сложную функцию y  где f x  f ; x   y ( x ) f  x, y   f  x, y   fx  f  f y , x y x fy  (1.4) f . y Тогда выражение (1.3), удерживая три члена разложения, представим в виде yk 1  yk  hf  xk , yk   h2   f x  xk , yk   f  xk , yk   f y  xk , yk   O ( h3 ), 2 (1.5)   где слагаемое O h 3 означает, что в следующие члены ряда величина h входит в степени не ниже третьей. Иначе говоря, если для нахождения приближенного значения 19 yk 1 использовать формулу (1.5) без слагаемого O ( h3 ) , то ошибка ограничения будет 3 приблизительно равна Kh , где K  некоторая постоянная. Решение задачи Коши с помощью ряда Тейлора является одноступенчатым методом, так как для вычисления решения yk 1 в очередной точке требуется информация только об одной предыдущей точке  xk , yk . Практическое неудобство и вычислительная неэффективность этого метода заключается в том, что при необходимости получить лучшее приближение с меньшей ошибкой ограничения придется вычислять производные высших порядков, характеризующиеся сложными выражениями, в особенности для структурно сложных функций. Так в рассматриваемом случае для разложения (1.5) с 4 остаточным членом вида O ( h ) необходимо будет использовать выражение для производной третьего порядка функции f  x, y  y  f xx  2 ff xy  f 2 f yy  f x f y  ff y2 . где f xx 2 f  2; x f xy 2 f  ; xy f yy 2 f  2. y Производные более высоких порядков будут иметь еще более сложные выражения. 2. Метод Эйлера Метод Эйлера - это приближенный метод решения задачи Коши для обыкновенного дифференциального уравнения первого порядка, 20 представленного в стандартной форме, разрешенной относительно производной y  f  x, y . (2.1) Начальное условие для такого дифференциального имеет вид y  x0   y0 . Интегральной кривой (2.2) дифференциального уравнения (2.1) называется график функции y  x  , являющейся решением этого уравнением. В геометрической интерпретации решение дифференциального уравнения (2.1) – это интегральная кривая L , проходящая через точку M 0  x0 , y0  . В методе Эйлера на оси абсцисс назначается координатная сетка в виде системы равноотстоящих точек с шагом h xi  x0  ih, i  0,1,2, (2.3) а искомая интегральная кривая L приближенно заменяется ломаной L  M 0 M 1M 2  с вершинами M i  xi , yi  , прямолинейные звенья которой M i , M i 1 имеют угловой коэффициент ki ki  yi 1  yi  f  xi , yi  , i  0,1,2, . h 21 (2.4) Ломаная L называется ломаной Эйлера. В соответствии с выражением (2.4) звенья M i , M i 1 ломаной Эйлера в каждой вершине M i имеют направления, которые определяются угловым коэффициентом yi  f  xi , yi  касательной к интегральной кривой уравнения (2.1) в точке M i . В связи с этим метод Эйлера называют также методом построения частного решения дифференциального уравнения первого порядка по схеме касательной аппроксимации искомой интегральной кривой. Действительно, через точку M 0  x0 , y0  , определяемую начальным условием (2.2), можно провести касательную к искомой интегральной кривой, проходящей через эту точку. Уравнение этой касательной имеет вид y  y0  f ( x0 , y0 )  x  x0 . (2.5) Перемещаясь по этой касательной на один шаг по оси x , можно попасть в точку M 1  x1 , y1  с координатами x1  x0  h ; y1  y0  f ( x0 , y0 ) h . При достаточно малой величине шага интегральной кривой точка M1 будет (2.6) h и нежесткой достаточно близко располагаться относительно искомой интегральной кривой. Через точку M 1 можно провести прямую с угловым коэффициентом, равным угловому коэффициенту f  x1, y1  квазикасательной к искомой интегральной кривой в точке с абсциссой x1 22 y  y1  f ( x1, y1 )  x  x1 . (2.7) Перемещаясь по этой прямой на один шаг по оси x , можно попасть в точку M 2  x2 , y2  с координатами x2  x0  2h ; y2  y1  f ( x1, y1 ) h , (2.8) которая будет располагаться вблизи соответствующей точки (с абсциссой x2 ) искомой интегральной кривой. Проделанные два перемещения по прямым, имеющим направления, определяемые касательными к искомой интегральной уравнения (2.1), соответствуют двум первым звеньям кривой M 0 M1 и M 1M 2 ломаной Эйлера. Таким образом, построению ломаной Эйлера соответствует процесс последовательного перемещения по отрезкам прямых M i M i 1 , угловые коэффициенты которых определяется как значения f  xi , yi  правой части дифференциального уравнения (2.1) в начальной точке M i каждого отрезка. Рис. 2.1 иллюстрирует геометрическое представление i -го шага метода Эйлера. 23 Рис. 2.1. Геометрическое представление метода Эйлера Приближенные значения yi решения y  x  уравнения (2.1) в узлах xi расчетной сетки в соответствии с зависимостями (2.4) – (2.6) определяются по следующей рекуррентной формуле: yi 1  yi  ki h, i  0,1,2,, (2.9) где ki  f  xi , yi . Приближенное решение дифференциального уравнения (2.1), определенное на дискретном множестве полученных точек xi , yi  , можно представить в виде непрерывной функции с помощью одной из схем интерполяции. В соответствии с рекуррентной формулой (2.9) метод Эйлера согласуется с рядом Тейлора (1.3) по точности приближения до члена 24 первого порядка относительно показателя степени шага интегрирования h . По функциональной схеме метод Эйлера является одноступенчатым и прямым, что позволяет считать его по общей классификации методом Рунге-Кутта. Оценка точности методов Рунге-Кутта осуществляется на основе сравнения с разложением Тейлора. Если приближения в методе осуществляются с точностью p до слагаемого пропорционального h , то степень p называется порядком метода. Поэтому классический метод Эйлера можно рассматривать как простейший метод Рунге-Кутта первого порядка. 25