2 семестр / ПОСОБИЕ_ВычМат

51

3.3 Аппроксимация функций

Задача аппроксимации функций относится к задачам среднеквадратичного дискретного приближения, при котором мерой отклонения функции ϕ(x) от заданной функции f(x) на множестве точек (xi, yi) (i = 1, …, n) является величина S, равная сумме квадратов разностей между значениями заданной и аппроксимирующей функций в данных точках:

n
S = å(( yi −ϕ( xi ))2 .	(11)

i=1

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

После того, как выбрано семейство функций g(x, α), осуществляется переход к решению второй задачи – определению наилучшего параметра α = (α0, ..., αm). В соответствии с мерой отклонения (11) для выбора наилучших параметров применим метод наименьших квадратов, который требует минимизации значения S.

Пусть в качестве семейства функций g(x, α) выбраны полиномы m-ой степени

m
Pm(x)= α0 + α1x +...+ αm xm = åα j x j .	(12)

j =0

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

n	m
S = å( yi	− åa j xij )2 → min.	(13)
i=1	j=0

Таким образом, для отыскания аппроксимирующего полинома необходимо определить m+1 его коэффициентов α0, α1,..., αm так, чтобы выполнился сформулированный критерий минимальности суммы (13). Следовательно, в точке α с необходимостью должны выполняться соотношения:

	∂ S	= 0;	∂ S	= 0,...,	∂ S		= 0 .	(14)
	∂ a		∂ a		∂ a	m
	0		1

Считая коэффициенты αi , i = 0, 1, 2,…, m, непрерывными величинами, выполним дифференцирование по αi функции (13) и, подставляя результаты в соотношения (14), получаем систему линейных уравнений m+1-го порядка относительно неизвестных α0, α1,..., αm:

						52
		¶ s			n	m					ü
					= 2 å ( yi	- åα j xij )			= 0,		ï
		¶ α0
					i=1	j=0					ï
		¶ s			n	m					ï
				= 2 å( yi -		åα j xij )xi = 0,					ï
	¶ α1
					i=1	j=0					ý
				.....................................							ï
											ï
	¶ s				n	m
							j	m			ï
				= 2 å( yi -		åα j xi		)xi		= 0.ï
¶ αm
					i=1	j=0					þ

После преобразований получим систему m+1 линейных уравнений:

		n	m	n	ü
	å å a j xij			= å yi ,	ï
	i=1		j =0	i=1	ï
	n	m		n	ï
	å å a j xij +1			= å yi xi , ï
	i=1	j =0		i=1	ý
					ï
		.............................
	n				ï
		m	j + m	n
					m ï
	å å a j xi			= å yi xi .
	i=1	j =0		i=1	ï
					þ

Система уравнений (15) в сокращенной записи имеет вид

n m	n
å å a j xij + k = å yi xik ; k = 0, 1, ..., m
i=1 j =0	i =1

(15)

(16)

и решается относительно неизвестных коэффициентов aj (j = 0,1,...,m), например, методом Гаусса.

3.3.1 Оценка точности аппроксимации

Рассмотрим следующее тождество:

yi – y = ( y€i – y )+( yi – y€i ),

где yi – значение функции в i-ом узле аппроксимации; y€i – значение функции регресии в i-ом узле аппроксимации;

y =

– среднее значение функции.

Возведя обе части этого выражения получим

n

å yi

i=1

n

в квадрат и просуммировав от i = 1 до n,

n	n	n
å(yi - y)2	= å(y€i - y)2	+ å(yi - y€i )2 .	(17)
i=1	i=1	i=1

53

Формула (17) играет фундаментальную роль в оценке точности аппроксимации функции. Согласно ей сумма квадратов относительно среднего равна сумме квадратов, обусловленной регрессией плюс сумма квадратов относительно регрессии. Отсюда следует, что разброс y-ов относительно их среднего можно приписать в некоторой степени тому факту, что не все значения функции лежат на линии регрессии. А если бы это было не так, то сумма квадратов относительно регрессии была бы равна нулю, т. е.

n

å(yi − y€i )2 = 0.

i=1

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

	n
	å(y€i − y)2
R2 =	i=1	(18)
	n
	å(yi − y)2

i=1

будет не слишком сильно отличать от единицы. Величину R2 называют множест-

венным коэффициентом детерминации или квадратом множественного коэф- фициента корреляции.

3.3.2Алгоритм и программа аппроксимации функции

Вкачестве аппроксимирующей функции выбираем полином степени m:

m

Pm(x, a) = åa j x j .

j =0

Алгоритм аппроксимации функции может быть представлен следующими шагами.

1.Ввести значения: n – число узлов аппроксимации, xi, yi (i = 1, 2,…, n) – значения аргумента и функции в i-ом узле аппроксимации, m – степень аппроксимирующего полинома Pm(x, a).

2.Вычислить коэффициенты матрицы системы линейных уравнений:

n

bkj = åxij +k , k = 0, 1, …, m, j = 0, 1, …, m.

i=1

3.Вычислить свободные члены системы линейных уравнений:

n

bk,m+1 = å yi xik , k = 0, 1, …, m.

i=1

4.Решить систему линейных уравнений методом Гаусса относительно неиз-

вестных a = (a0, a1, …, am).

5.Вычислить величину R2 по формуле (18).

54

6.Вывод результатов: a = (a0, a1, …, am) – коэффициенты аппроксимирующего полинома; R2 – коэффициент множественной детерминации.

7.Конец вычислений.

Пример 3. Функцию f(x), определенную таблицей, аппроксимировать многочленом второй степени, оценить погрешности аппроксимации.

X	xi	1	2	3	4	5
Y	yi	2	2,69	3,1	3,39	3,61

3.3.2.1 Решение в системе Turbo Pascal

program Approksimazia; uses crt;

const

m=2; {наибольшая степень полинома в нашей программе} n=5;{максимальное число узлов аппроксимации}

type

TB=array[0..m+1,0..m+1] of real; {тип расширенной матрицы системы линейных уравнений}

TA=array[0..m] of real;{тип коэффициэнтов аппроксимирующего полинома} TXY=array[1..n] of real;{тип элементов улов аппроксимации}

{процедура решения системы линейных уравнений методом Гаусса} procedure pGauss(A:TB; n: byte; Var x:TA);

var

k,i,j:byte; d:real; begin

{прямой ход в методе Гаусса}

for k:=0 to n-1 do {цикл удаления k-й неизвестной} for i:=k+1 to n do{из i-го уравнения}

begin

d:=A[i,k]; for j:=k to n+1 do A[i,j]:=A[i,j]-A[k,j]/A[k,k]*d; end;

{обратный ход в методе Гаусса} for k:=n downto 0 do

begin

x[k]:=A[k,n+1]; for j:=k+1 to n do x[k]:=x[k]-A[k,j]*x[j]; x[k]:=x[k]/A[k,k]; end;

end;

{функция вычисления k-ой степени числа x<>0} function sq(x:real; k:integer):real;

var i:integer; p:real; begin

p:=1; for i:=1 to k do p:=x*p; sq:=p; end;

55

{основная программа} const

x:TXY = (1, 2, 3, 4, 5); {значения аргумента в узлах аппроксимации}

y:TXY = (2, 2.69, 3.1, 3.39, 3.61); {значения функции в узлах аппроксимации} var

b:TB; {расширенная матрица системы линейных уравнений} a:TA; {коэффициэты полинома}

i,j,k:integer; {для индексации} R:real; {для оценки аппроксимации}

Ycp,p,s1,s2:real; {рабочие переменные} begin

{занулить матрицу b}

for k:=0 to m+1 do for j:=0 to m+1 do b[k,j]:=0; {формирование расширенной матрицы}

for k := 0 to m do begin

for i:=1 to n do

b[k,m+1]:=b[k,m+1]+y[i]*sq(x[i],k); {формирование свободных членов уравнений системы}

for j := 0 to m do {формирование коэффициентов системы уравнений} for i:=1 to n do {цикл по числу узлов аппроксимации}

begin	p:=sq(x[i], k+j);	b[k,j] := b[k,j] + p;	end;
end;

pGauss(b,m,a);{решить систему уравнений метомом Гаусса} {вывод на экран коэффициэнтов аппроксимирующего полинома} for j:=0 to m do writeln('a[',j,']=',a[j]:10:5);

{вычислить среднее значение функции}

Ycp:=0; for i:=1 to n do Ycp:=Ycp+y[i]; Ycp:=Ycp/n; {формирование оценки аппроксимации}

s1:=0; s2:=0; for i:=1 to n do begin

{вычисление суммы квадратов относительно среднего} s1:=s1+sqr(y[i]-Ycp);

{вычисление значения аппроксимирующего полинома в узле x[i] по схеме Горнера}

p:=a[m]; for j:=m-1 downto 0 do p:=p*x[i]+a[j]; {вычисление суммы квадратов, обусловленной регрессией} s2:=s2+sqr(p-Ycp);

end;

R:=s2/s1;

writeln('Множественный коэффициент детерминации R^2 ',R:10:5); readkey;

end.

56

3.3.2.2 Алгоритм решения в электронной таблице Excel

1.Поместить координаты вектора X в диапазон ячеек C2:G2.

2.Поместить координаты вектора Y в диапазон ячеек C3:G3.

3.Выделить диапазон ячеек C2:G3.

4.Вызвать МАСТЕР ДИАГРАММ.

4.1.В диалоговом окне ТИП ДИАГРАММЫ выбрать тип – ТОЧЕЧНАЯ, вид – ТОЧЕЧНАЯ.

4.2.Нажимаем кнопку ДАЛЕЕ, …, а затем ГОТОВО.

4.3.Диаграмму размещаем на отдельном листе.

5.Выделить диаграмму. В меню ДИАГРАММА выбрать команду ДОБАВИТЬ ЛИНИЮ ТРЕНДА. В диалоговом окне ЛИНИЯ ТРЕНДА на вкладке ТИП выбрать тип линии ПОЛИНОМИАЛЬНАЯ степень 2. В диалоговом окне ЛИНИЯ ТРЕНДА на вкладке ПАРАМЕТРЫ установить флажки ПОКАЗЫВАТЬ УРАВНЕНИЕ и ПОМЕСТИТЬ ВЕЛИЧИНУ ДОСТОВЕРНОСТИ АППРОКСИМАЦИИ и нажать кнопку ОК. Для обеспечения наглядности выполнить форматирование диаграммы.

6.Конец вычислений.

Задания к данной теме приведены в приложении В.

57

4 ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ

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

4.1 Вычисление производной по ее определению

Пусть функция y = f(x) определена в некоторой окрестности точки x0 и имеет производную в этой точке, т. е. существует предел отношения приращения функции y к приращению аргумента x при стремлении x к нулю

y ′(x0)=f ′(x0)= lim	y	, x = x - x0, y = f(x0 + x) - f(x).	(1)
x→0	x

Значение производной можно получить, переходя к пределу (1) по последовательности целых чисел n, полагая, например,

x = ( x)n= (10xn)0 ,

где ( x)0 – некоторое начальное приращение аргумента, n = 1, 2, ... . Тогда значение производной функции f(x) в точке x0 запишется так:

y ′(x0)=f ′(x0)= lim	(	y )n	, ( y)n = f(x0 + ( x)n) - f(x0).
		x )n
n→∞ (

Отсюда получаем приближенное равенство

y ′(x0)≈	(	y )n	.	(2)
	(
		x )
		n

Для достижения заданной точности ε приближения производной при определенном числе вычислений можно воспользоваться неравенством

	(	y )n	−	(	y )n−1		< ε .
	(	x )		(	x )
		n			n−1

Тогда при вычислении производной по ее определению находят отношение

			58
	y	=	f ( x + x ) − f ( x )
	Dx		Dx

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

4.2 Конечно-разностные аппроксимации производных

Пусть отрезок [a, b] разбит на n (n ³ 2) равных частей точками

{xi}: a = x0 < x1 < ... < xi-1 < xi < xi+1 < ... < xn = b.

Разность между соседними значениями аргумента постоянна, т. е. шаг h = xi – xi-1 = const (i = 1, 2, ..., n). Далее пусть на отрезке [a, b] определена функция y = f(x), значения которой в точках xi равны

yi = f(xi) (i = 0, 1, ..., n) .

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

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

y ¢(xi)»	yi , xi = xi+1 – xi = h, yi = yi+1 – yi,,
	x
	i
y ¢(xi)»		yi+1 − yi	, i = 0, 1, ..., n-1;	(3)
		h

б) аппроксимация с помощью левых разностей:

y ¢(xi)»		yi	, xi = xi-1 – xi = -h, yi = yi-1 – yi ,
		x
		i
y ¢(xi)»	yi −1	− yi		=	yi − yi −1	, i = 1, 2, ..., n;	(4)
					h
	− h

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

y ¢(xi)»	yi , xi = xi+1 – xi-1 = 2h, yi = yi+1 – yi-1 ,,
	x
	i
	y ¢(xi)»	yi+1 − yi−1	, i = 1, 2, ..., n-1.	(5)
		2h

Отметим, что соотношение (3) и (5) не позволяют вычислить производную в точке xn = b, а (4) и (5) – в точке x0 = a.

Погрешность аппроксимации производных с помощью левых и правых разностей имеет один и тот же порядок 0(h), а погрешность аппроксимации центральными разностями имеет порядок 0(h2).

59

Приближенное значение производной второго порядка в точке xi выражается через значения функции yi-1, yi, yi+1. Для этого представим вторую производную с помощью правой разности:

y ′′(xi)≈

y'i

,

xi = xi+1 – xi = h,

y′i = y′i+1 – y′i ,

x

i

а приближенные

производные первого порядка y′i+1, y′i

– с помощью левых раз-

ностей:

'

'

yi+1 − yi

'

'

yi − yi−1

yi

+1 = y

( xi+1 ) =

, yi = y

( xi ) ≈

.

h

h

Окончательно получим:

y'i+1 − y'i

yi+1 − yi

− yi − yi−1

y'' ( xi ) =

=

h

h

=

yi+1 − 2 yi + yi−1

,

h

h

h2

где i = 1, 2, ..., n–1.

Погрешность аппроксимации производной второго порядка имеет порядок 0(h2) для функции f(x), имеющей непрерывную производную до четвертого порядка включительно на отрезке [a, b]. Для повышения точности вычислений строят интерполяционный полином, дифференцируя который получают соответствующие формулы вычисления производных.

4.3 Другие методы вычисления производных первого и второго порядков

Пусть отрезок [a, b] разбит на n равных частей с шагом h= b −n a , и в точках

xi (i = 0,1,...,n; x0 = a, xn = b) некоторая функция принимает значения yi. Для переменной x, принадлежащей отрезку [a, b], требуется вычислить значения первой и второй производной, имеющих порядок аппроксимации 0(h2) в зависимости от шага. Полагаем, что значения производных y′ и y′′ в точках x, близких к точкам xi, равны соответствующим значениям y′i и y′′i. Будем считать точку x близкой к xi, если она принадлежит промежутку

xi - 2h ≤ x≤ xi + 2h .	xi-	h	xi	xi +	h	xi+1	x
	2			2

Точки x, близкие к xi, имеют одно и то же значение параметра:

60

i = целая часть( x −h a + h2 ).

В зависимости от i и n ≥ 3 имеем три типа формул. Так, при i = 0:

y′0 ≈		1		(–3y0 + 4y1 – y2),	y′′0 ≈		1		(2y0 – 5y1 + 4y2 – y3);		(6)
		2h					h2
при i = 1, 2, ..., n–1:					1
y′i ≈	1		(-yi-1 + yi+1), y′′i ≈			(yi-1 – 2yi + yi+1);					(7)
	2h				h2
при i = n:								1
y′n ≈		1		(yn-2 + 4yn-1 + 3yn),	y′′n ≈					(–yn-3 + 4yn-2 - 5yn-1 + 2yn).	(8)
		2h						h2

4.3.1 Алгоритм вычисления первой и второй производной функции

1.Задать начальные условия: границы отрезка [a, b], число разбиений отрезка.

2.Задать значение x, при котором необходимо вычислить производную функции.

3.Задать значения у в узлах интерполяции для вычисления функции.

4.Вычислить значения h, i.

5.Вычислить значения первой и второй производной у1 и у2.

6.Вывод результатов: x, y1,y2.

7.Конец вычислений.

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

x	2,4	2,7	3,1	3,4	3,7	4,1	4,4	4,7	5,1	для
y	3,53	3,78	3,95	4,04	4,12	4,16	4,25	4,37	4,51	x=4,13

Вычисления выполнить по формулам (6)–(8).