4.7 Аппроксимация с помощью взвешенных невязок

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

Пусть требуется построить аппроксимацию функции f(x) на отрезке [a, b]. Аппроксимирующая функция g(x) берется в виде линейной комбинации

, (4.9)

где {g_k} – некоторые параметры, вычисляемые таким образом, чтобы получить хорошее приближение, {N_k , k=1, 2, 3, …} – полная система линейно независимых базисных функций. Условие полноты этой системы означает, что она способна приближать g(x) к f(x) с любой наперед заданной точностью при M.

Вводится функция R – погрешность аппроксимации или невязка, определяемая в каждой точке промежутка [a, b] как разность f и g:

, axb. (4.10)

Понятно, что для хорошего приближения невязка должна как можно меньше отличаться от нуля, и задача формулируется следующим образом: подобрать коэффициенты из (4.9) так, чтобы R(x)0 для всех x[a, b]. Чтобы уменьшить эту невязку (по модулю) на всем указанном промежутке, потребуем равенства нулю соответствующего числа интегралов от невязки, взятых с различными весами, т.е.

, i=1, 2, …, M, (4.11)

где {W_i , i=1, 2, 3, …} – множество линейно независимых весовых функций. Тогда общее требование сходимости R0 при M можно записать в форме условия выполнения равенств (4.11) для всех i при M. В функциональном анализе показывается, что это действительно будет верно лишь при условии R0 во всех точках отрезка [a, b].

С учетом выражения (4.10) система уравнений метода взвешенных невязок сводится к системе линейных алгебраических уравнений для неизвестных коэффициентов

; i=1, 2, …, M, (4.12)

или в общем виде –

, (4.12*)

где

; ; ; 1i, kM.

Таким образом, если известна аппроксимируемая функция f и выбраны подходящие системы базисных и весовых функций, то, решая систему (4.12), можно получить коэффициенты {g_k} в аппроксимации (4.9).

На практике могут быть использованы различные системы весовых функций {W_i , i=1, 2, …}, ведущие к разным частным методам аппроксимации посредством взвешенных невязок. Ниже приводятся некоторые из наиболее употребительных таких систем.

1) Поточечная коллокация. Весовые функции заданы формулой

,

где (x–x_i) – дельта-функция Дирака, по определению обладающая свойствами:

(x–x_i)=0 при xx_i; (x–x_i)= при x=x_i; _.

Такой выбор весовых функций эквивалентен тому, что невязка полагается равной нулю (т.е. f=g) в ряде фиксированных точек {x_i}. А последнее требование означает, что в данном случае реализуется интерполяция.

2) Коллокация по отрезкам. Пусть на отрезке [a, b] введена система точек {x_i}: x₁=a<x₂<x₃<…<x_M<x_M+1=b. Если весовые функции определены по правилу

то уравнения метода взвешенных невязок равносильны требованию равенства нулю интегралов от невязки по каждому из M промежутков [x_i, x_i+1]. Другими словами, площади под кривыми y=f(x) и y=g(x) на каждом таком промежутке должны быть равны. Матрица S и вектор f в (4.12*) имеют элементы

; .

3) Метод Галеркина. В этом наиболее популярном варианте метода взвешенных невязок в качестве весовых функций используются сами базисные функции, т.е. W_i=N_i. Матричные элементы системы метода взвешенных невязок и элементы её вектора правых частей имеют вид

; .

Отметим, что здесь матрица S симметричная и положительно определенная, что дает этому методу дополнительные вычислительные преимущества.

Особая простота уравнений метода Галеркина возникает, когда функции {N_k , k=1, 2, …} на [a, b] образуют ортогональную систему, т.е. выполняется соотношение

, d0,

(_ik – символ Кронекера; _ik=0 при ik и _ik=1 при i=k ). Например, рассматривая аппроксимацию на отрезке [0, l] базисными функциями { , k=1, 2, …}, получим для коэффициентов системы:

.

Следовательно, матрица S системы (4.12*) имеет диагональную форму, что сразу позволяет записать решение

, k=1, 2, …, M.

Полученное выражение тождественно формуле для коэффициентов разложения функций по синусам кратных дуг из теории рядов Фурье. Вообще говоря, известные из математики ортогональные разложения можно рассматривать как аппроксимацию по Галёркину с помощью взвешенных невязок при М.

4) Метод моментов. Один из очевидных способов выбора весовых функций состоит в использовании системы {W_i=x^i–1, i=1, 2, 3, …}. Тогда метод взвешенных невязок приводит в требованию, чтобы площадь под кривой невязки R и её различные моменты относительно начала координат на отрезке [a, b] были равны нулю, т.е.

, , , …, .

Коэффициенты системы уравнений метода выражаются через интегралы

; .

Что касается термина “момент”, то помимо известных величин из механики (например, момент силы), можно также проследить аналогию с теорией вероятностей, где вводится понятие начального момента k-го порядка случайной величины X:

.

Замечания. 1) В ряде случаев целесообразно вместо (4.9) брать аппроксимацию вида

, (4.13)

где (x) – любая гладкая функция, принимающая на концах отрезка [a, b] те же значения, что и f(x), т.е. (a)=f(a), (b)=f(b); при этом базисные функции {N_k} должны равняться нулю при x=a и x=b. Тогда в системе уравнений метода взвешенных невязок (4.12*) изменится вычисление только элементов вектора правых частей, а именно

, i=1, 2, …, M.

В качестве  обычно берут линейную функцию

. (4.14)

Данный способ аппроксимации позволяет лучше учесть поведение функции вблизи границ отрезка. К примеру, предположим, что на [a, b] производится разложение с помощью усеченных синус-рядов Фурье

.

Заметив, что все базисные функции, участвующие в этом разложении, обращаются в ноль при x=a и x=b, получаем g(a)=g(b)=0. При этом сама аппроксимируемая функция f(x) в этих точках может быть и не равной нулю. Тогда около значений a и b будем иметь заметную ошибку аппроксимации. Однако, выбрав аппроксимацию в виде

,

обеспечим и g(a)=f(a), и g(b)=f(b) и тем самым лучшее соответствие вблизи граничных точек.

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

В заключение рассмотрим примеры решения задач на основе метода взвешенных невязок в системе Maple.

Пример. Аппроксимировать функцию на промежутке [0, ].

Решение. Выбор весовых функций осуществим по Галеркину, число параметров возьмем M=4. Составляем программу:

> restart;

> f:=x->int(exp(-2*sin(t)),t=0..x);

> N:=[1,x,x^2,x^3,x^4,x^5]; # базисные функции

> W:=N; # весовые функции как в методе Галеркина

> m:=4; # число параметров аппроксимации

> g:=x->sum(c[k]*N[k],k=1..m); # аппроксимация в виде суммы (4.9)

> a:=0: b:=Pi: diap:=x=a..b; # диапазон изменения переменной x

> for i from 1 to m do eqns[i]:=(int(g(x)*W[i],diap)=evalf(Int(f(x)*W[i],diap))); od; # формирование системы уравнений метода взвешенных невязок

> r:=solve({eqns[j]$j=1..m},{c[j]$j=1..m}); #решение системы

{c₄ = 0.087577824, c₃ = -0.41270077, c₂ =0.75908203, c₁ = 0.023952682}

> plot([subs(r,g(x)),f(x)],x=a..b); # графики аппроксимирующей и исходной функций

Рис. 4.10.

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

.

Чем меньше эта величина, тем лучше аппроксимация. Её вычисление в Maple оформляем в виде строки:

> delta:=evalf(int((subs(r,g(x))-f(x))^2,diap))^0.5;

Если в представленном тексте программы поменять число параметров m на 5, то величина  не изменится, но с шестью параметрами она резко уменьшается: =0.000724. Это свидетельствует о сходимости метода.

Может возникнуть вопрос: зачем нужен метод взвешенных невязок, если можно функцию f(x) разложить в степенной ряд и его частичную сумму использовать в качестве аппроксимирующей функции? Для сравнения выведем первые четыре члена разложения нашей функции в ряд в окрестности точки x=/2 (середина отрезка).

> h:=convert(series(1.*f(x),x=Pi/2,4),polynom);

> plot([h,f(x)],diap,thickness=3);

Рис. 4.11.

Хорошо видно, что приближение тем хуже, чем дальше от точки разложения x=/2. Отсюда вывод: аппроксимация методом взвешенных невязок вносит более равномерную ошибку на заданном промежутке и, как правило, дает лучший результат, нежели с помощью разложения в степенной ряд (при одинаковой степени аппроксимирующего полинома).

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

> a:=0;b:=Pi; # начало и конец промежутка изменения x

> W:=[Dirac(x-('k'-1)*(b-a)/('m'-1))$'k'=1..m]; # дельта-функции Дирака

И ещё, для того чтобы Maple корректно вычислял интегралы с дельта-функциями Дирака, необходимо расширить пределы интегрирования, например, задать

> diap:=x=a-1..b+1;

Коллокацию по отрезкам получим, если задать

> a:=0;b:=Pi;

> W:=[piecewise(x>a+('k'-1)*(b-a)/m and

x<a+'k'*(b-a)/m,1,0)$'k'=1..m]; # единичные запаздывающие импульсы

Пример. Получить аппроксимацию решения задачи Коши для дифференциального уравнения , , на отрезке [0.5, 5].

Решение. Попробуем решить это уравнение аналитически

> restart;

> p:=dsolve({diff(x^2*diff(y(x),x),x)+(x^2-1)*y(x)=0,y(0.5)=0.,D(y)(0.5)=2.5},y(x));

Maple справился с этой задачей, но решение выдал в виде громоздкого выражения через специальные функции (Бесселя). Работать с таким выражением крайне неудобно. Вот здесь и приходит на помощь аппроксимация. Сначала введем функцию для хранения точного решения

> u:=unapply(subs(p,y(x)),x);

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

> a:=0.5: b:=5: m:=4: # т.е. взято четыре параметра

> N:=[1,x,x^2,x^3,x^4,x^5,x^6]: W:=N: # весовые функции совпадают с базисными

> g:=x->sum(c[k]*N[k],k=1..m): # стандартная аппроксимация

> for i from 1 to m do eqns[i]:=(int(g(x)*W[i],x=a..b)=

evalf(Int(u(x)*W[i],x=a..b)));od;

> r:=solve({eqns[j]$j=1..m},{c[j]$j=1..m});

> subs(r,g(x));

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

> plot([subs(r,g(x)),u(x)],x=a..b,thickness=2);

Рис. 4.12.

Графики почти совпадают. Погрешность аппроксимации:

> delta:=evalf(Int((subs(r,g(x))-u(x))^2,x=a..b,method=_CCquad))^0.5;

Если метод сходится, то теоретически, увеличивая число параметров M (в программе это переменная m), должны получать более точное приближение, а значит, уменьшающуюся величину . Составим таблицу.

M	Аппроксимирующая функция g(x)	Погрешность 
4	–0.6622+1.923x–0.7123x^2+0.0688x^3	0.03905
5	–0.8572+2.359x–1.014x^2+0.1496x^3–0.0073x^4	0.02450
6	–1.064+2.955x–1.600x^2+0.4037x^3–0.0573x^4+0.0036x^5	0.01765
7	–1.392+4.112x–3.084x^2+1.311x^3–0.342x^4+0.048x^5–0.0027x^6	0.00979

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

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

Взяв в рассматриваемом примере аппроксимацию, удовлетворяющую граничным условиям на концах промежутка [a, b],

> g:=(u(b)-u(a))*(x-a)/(b-a)+u(a)+ sum(c[k]*N[k]*(x-a)*(x-b),k=1..m);

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

Разобранный пример касался случая, когда Maple был способен найти аналитическое решение дифференциального уравнения. Но что делать, когда такого решения нет и результат выдается в численном виде? Оказывается, часто удаётся построить аппроксимацию методом взвешенных невязок и в этом случае.

Пример. Аппроксимировать решение задачи Коши , при 0x4.

Решение

> restart;

Решение ищем в численном варианте (с параметром type=numeric):

> p:=dsolve({D(y)(x)+y(x)=sin(x*y(x)),y(0)=0.5}, y(x),type=numeric);

Решение хранится в переменной p, но его формат весьма необычный. Это видно из примера вывода значения функции в некоторой точке:

> p(2.3);

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

> u:=proc(t) global i;local w; w:=subs(p(t)[2],y(x)); RETURN(w) end;

Чтобы проверить, снова вычислим:

> u(2.3)

Цель достигнута. Теперь реализация самого метода:

> m:=5: a:=0: b:=4.:

Выбираем системы базисных и весовых функций

> N:=x->[sin(Pi*'k'*(x-a)/(b-a))$'k'=1..m];

> W:=x->[1,x,x^2,x^3,x^4,x^5,x^6];

Обратите внимание, здесь в отличие от предыдущих примеров массивы M и N введены как вектор-функции. Строим аппроксимацию с совпадением граничных значений (см. (4.13)–(4.14)):

> g:=x->(u(b)-u(a))*(x-a)/(b-a)+u(a)+ sum(c[k]*N(x)[k],k=1..m);

> for i from 1 to m do

eqns[i]:=int(g(x)*W(x)[i],x=a..b)=

evalf(Int(W(t)[i]*'u'(t),a..b));

od; # система уравнений метода

> r:=solve({eqns[j]$j=1..m},{c[j]$j=1..m}); # решение системы

> evalf[6](subs(r,g(x))); # получена аппроксимирующая функция

⁺

> gr1:=plots[odeplot](p,a..b): # графики функций

> gr2:=plot([subs(r,g(x))], x=a..b,color=blue,thickness=2):

> plots[display](gr1,gr2);

Рис. 4.13.

Соответствие имеет место. Для количественной оценки погрешности вычисляем:

> delta:=evalf(Int((subs(r,g(t))-'u'(t))^2, a..b,method=_Gquad))^0.5;

:=0.009961375

Проверку сходимости проведите самостоятельно.