| Автор cайта: Владимир Потемкин fortran-90@yandex.ru |
 |
Интерполяция, сглаживание и аппроксимация
PL10. Вычисление значения функции с помощью интерполяционного процесса Эйткена-Лагранжа.
PL12. Вычисление значения функции с помощью интерполяционного процесса Эйткена-Эрмита.
SPL. Вычисление коэффициентов кубического интерполяционного сплайна.
SPS. Вычисление коэффициентов кубического сглаживающего сплайна.
Hermite. Построение аппроксимирующей кирвой Эрмита.
ORT. Аппроксимация функциональной зависимости полиномом заданной степени.
PL10
Программа PL10 вычисляет значение функции y(x) в заданной точке x с помощью интерполяционного процесса Эйткена-Лагранжа.
Функция y(x) задана таблично массивами Xi и Yi, i=1, 2, ..., n. Массив Xi должен быть упорядочен в возрастающем порядке, т.е. для всех его элементов xi<xi+1. Максимальное число узлов m, используемых для интерполяции, задается при вызове программы.
Процесс интерполяции построен по итерационной схеме. В качестве первого приближения используется значение функции в ближайшем узле интерполяции. Для получения более точных приближений используются последовательно вычисляемые значения многочлена Лагранжа L(k)(x), степень которого возрастает на каждом шаге итерации:
где ω(k)(x) = (x-x1)(x-x2) ... (x-xk+1), k=1, 2, ..., n-1. Итерационный процесс продолжается до получения заданной точности интерполяции.
Вычисления проводятся по треугольной интерполяционной схеме Эйткена:
x1 y1
x2 y2 L12
x3 y3 L13 L123
x4 y4 L14 L124 L1234
. . . . . . . . . . . . . . . . . .
xm ym L1m L12m L123m ... L123...m
Каждый многочлен Лагранжа более высокой степени получается из двух предыдущих многочленов по рекуррентным формулам:
Элементы L12...k нисходящей диагонали этой схемы представляют собой интерполяционные многочлены Лагранжа L(k)(x) последовательно возрастающих степеней (k=2, ..., m).
Для определения значений многочлена L(k)(x) при росте k в вычислительный процесс вовлекаются все более удаленные от x узлы интерполяции. Поэтому значения xk и yk (k=1, 2, ..., m) предварительно упорядочиваются таким образом, чтобы модули разностей |xk - x| возрастали при росте индекса k:
|xk - x| ≤ |xk+1 - x|.
Вызов программы
CALL PL10(X, Y, XX, M, Eps, YY, Error)
Параметры программы
X, Y, XX, M, Eps - входные параметры;
YY, Error - выходные параметры;
Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке;
Real XX - значение аргумента, в котором требуется вычислить значение функции;
Integer M - число узлов, которое разрешается использовать при интерполяции, 1<=M<=N;
Real Eps - погрешность, с которой требуется вычислить значение функции, X[1]<=XX<=X[N];
Real YY - значение функции, вычисляемое программой;
Integer Error - индикатор ошибки, принимающий при выполнении программы одно из следующих значений:
Error=0 - заданная точность интерполяции достигнута, ошибок нет;
Error=1 - заданная точность не достигается или достигнутая точность не может быть проверена, т.к. мало значение M;
Error=2 - заданная точность не достигается из-за роста вычислительной погрешности или плохого поведения старших производных интерполируемой функции.
Пример
! Интерполяция функции многочленом Лагранжа
! по интерполяционной схеме Эйткена.
program TestPL10
use NML
implicit none
integer, parameter:: n=40
real:: X(n), Y(n), XX, YY, YE, Eps
integer:: M, Er
integer i
!begin
M=6
Eps=1.E-6
do i=1, n
X(i)=0.5*float(i-1)
Y(i)=exp(-X(i))*sin(X(i))
end do
XX=4.2
! YE - точное значение по формуле
YE=exp(-XX)*sin(XX)
! YY - значение, полученное интерполяцией
call PL10(X, Y, XX, M, Eps, YY, Er)
print 10, Er
print 11, XX, YE, YY
10 format(/' Error =',I2/)
11 format(' X = ',F6.2,' YE =',F11.6,' YY =',F11.6/)
end program TestPL10
────────────
Error = 1
X = 4.20 YE = -0.013070 YY = -0.013074
Вернуться к оглавлению Скачать PL10
PL12
Программа PL12 вычисляет значение функции y(x) в заданной точке x с помощью интерполяционного процесса Эйткена-Эрмита.
Функция y(x) задана таблично массивами Xi и Yi, i=1, 2, ..., n. Массив Xi должен быть упорядочен в возрастающем порядке, т.е. для всех его элементов xi<xi+1. Максимальное число узлов m, используемых для интерполяции, задается при вызове программы.
Процесс интерполяции построен по итерационной схеме. В качестве первого приближения используется интерполяционный многочлен Эрмита первой степени. На каждой новой итерации степень многочлена увеличивается на 1. Итерационный процесс продолжается до получения заданной точности.
Для построения многочленов Эрмита нечетной степени 2k-1 используются значения функции и ее производной во всех k узлах, для построения многочленов четной степени 2k используется значение функции в k узлах и значение производной в k-1 узле.
Вычисления производятся по треугольной итерационной схеме Эйткена:
x1 y1 V1 V12 V13 V14 V15 ... V1,2m-1
y1' V2 V23 V24 V25 ...
x2 y2 V3 V34 V35 ...
y2' V4 V45 ...
x3 y3 V5 ...
y3' ...
. . . . . . . . . . . . . . . . . .
xm ym V2m-1
ym'
Элементы V1, V12, ..., V2m-1 первой строки этой схемы представляют собой значения многочленов Эрмита последовательно возрастающих степеней:
V1 = H(1)(x), V12 = H(2)(x), ..., V1,2n-1 = H(2m-1)(x).
Элементы первого столбца V1, V2, ..., V2m-1 вычисляются по формулам:
Vi = yj + yj'(x - xj)
Vi+1 = yj + (yj+1 - yj)[(x - xj)/(xj+1 - xj)],
где i=1, 3, 5, ..., 2m-1, и j=(i+1)/2. Остальные элементы этой схемы вычисляются (вверх по диагонали) по следующей рекуррентной формуле:
где j=(i+1)/2, l=(k+3)/2.
Для определения значений многчлена H(k)(x) при росте k в вычислительный процесс вовлекаются все более удаленные от x узлы интерполяции. Поэтому значения xk, yk и yk' (k=1, 2, ..., m) предварительно упорядочиваются таким образом, чтобы модули разностей |xk - x| возрастали при росте индекса k:
|xk - x| ≤ |xk+1 - x|.
Вызов программы
CALL PL12(X, Y, DY, XX, M, Eps, YY, Error)
Параметры программы
X, Y, DY, XX, M, Eps - входные параметры;
YY, Error - выходные параметры;
Real X(1:N), Y(1:N), DY(1:N) - массивы, содержащие значения аргумента, функции и производной. Массив X должен быть упорядочен в возрастающем порядке;
Real XX - значение аргумента, в котором требуется вычислить значение функции;
Integer M - число узлов, которое разрешается использовать при интерполяции, 1<=M<=N;
Real Eps - погрешность, с которой требуется вычислить значение функции, X[1]<=XX<=X[N];
Real YY - значение функции, вычисляемое программой;
Integer Error - индикатор ошибки, принимающий при выполнении программы одно из следующих значений:
Error=0 - заданная точность интерполяции достигнута, ошибок нет;
Error=1 - заданная точность не достигается или достигнутая точность не может быть проверена, т.к. мало значение M;
Error=2 - заданная точность не достигается из-за роста вычислительной погрешности или плохого поведения старших производных интерполируемой функции.
Пример
! Интерполяция функции многочленом Эрмита
! по интерполяционной схеме Эйткена.
program TestPL12
use NML_PL12
implicit none
integer, parameter:: n=40
real:: X(n), Y(n), DY(n), XX, YY, YE, Eps
integer:: M, Er
integer i
!begin
M=6
Eps=1.E-8
do i=1, n
X(i)=0.5*float(i-1)
Y(i)=exp(-X(i))*sin(X(i))
DY(i)=exp(-X(i))*(-sin(X(i))+cos(X(i)))
end do
XX=5.2
YE=exp(-XX)*sin(XX)
call PL12(X, Y, DY, XX, M, Eps, YY, Er)
print 10, Er
print 11, XX, YE, YY
10 format(/' Error =',I2/)
11 format(' X = ',F6.2,' YE =',F11.6,' YY =',F11.6/)
end program TestPL12
────────────
Error = 0
X = 5.20 YE = -0.004874 YY = -0.004874
Вернуться к оглавлению Скачать PL12
SPL
Программа SPL вычисляет коэффициенты кубического интерполяционного сплайна для функциональной зависимости yi = y(xi), i=1, 2, ..., n. Кубический сплайн представляет собой многочлен третьей степени на каждом из отрезков [xi, xi+1], i=1, 2, ..., n-1:
│s1(x) = y1 + b1(x - x1) + c1(x - x1)2 + d1(x - x1)3
│ . . . . . . . . . . . . . . . . . . . .(1)
│sn-1(x) = yn-1 + bn-1(x - xn-1) + cn-1(x - xn-1)2 + dn-1(x - xn-1)3
Для нахождения коэффициентов сплайна bi, ci, di используют условия сшивания для самих многочленов и для их первой и второй производной в точке xi:
│s1(x2) = y1 + b1h1 + c1h12 + d1h13 = y2
│ . . . . . . . . . . . . . . . . . . . .(2)
│sn-1(xn) = yn-1 + bn-1hn-1 + cn-1hn-12 + dn-1hn-13 = yn
│s1'(x2) = s2'(x2)
│ . . . . . . . . . . .(3)
│sn-2'(xn-1) = sn-1'(xn-1)
│s1"(x2) = s2"(x2)
│ . . . . . . . . . . .(4)
│sn-2"(xn-1) = sn-1"(xn-1)
где h1 = x2 - x1, ..., hn-1 = xn - xn-1. Из (4) получаем:
│c1 + 3d1h1 = c2
│ . . . . . . .
│cn-2 + 3dn-2hn-2 = cn-1
и дополняем уравнением cn-1 + 3dn-1hn-1 = cn. Вводим обозначения
│c1 = 3σ1
│ . . .
│cn = 3σn
│∆1 = (y2 - y1)/h1
│ . . . . . . .
│∆n-1 = (yn - yn-1)/hn-1
Имеем с помощью (2):
│bi = ∆i - hi(σi+1 + 2σi)
│ci = 3σi i=1, ..., n-1
│di = (σi+1 - σi)/hi
Подставляем в (3) с учетом (2):
│h1σ1 + 2(h1 + h2)σ2 + h2σ3 = ∆2 - ∆1
│h2σ2 + 2(h2 + h3)σ3 + h3σ4 = ∆3 - ∆2
│ . . . . . . . . . . . . . . . . . . . .(5)
│hn-3σn-3 + 2(hn-3 + hn-2)σn-2 + hn-2σn-1 = ∆n-2 - ∆n-3
│hn-2σn-2 + 2(hn-2 + hn-1)σn-1 + hn-1σn = ∆n-1 - ∆n-2
Мы получили систему из n-2 уравнений и n неизвестных.
Дополним полученную систему (5) уравнениями s1"(x1) = 0 и sn-1"(xn) = 0, или, что то же самое, σ1=σn=0. Тогда получаем полную систему для определения неизвестных bi, ci, di коэффициентов сплайна (1). Матрица полученной системы является симметричной и с диагональным преобладанием.
Построенный сплайн минимизирует функционал
b
J(s) = ∫(s"(x))2dx,
a
где a=x1, b=xn.
Применение краевых условий s1"(x1)=0 и sn-1"(xn)=0 не является оптимальным, поскольку в этом случае значение сплайна вблизи граничных точек может сильно отличаться от значения интерполируемой функции. Поэтому в программе применяются следующие типы краевых условий:
- краевые условия 1-го типа. s'(a) = y'(a), s'(b) = y'(b) - на концах промежутка [a,b] задаются значения 1-ой производной интерполируемой функции;
- краевые условия 2-го типа. s"(a) = y"(a), s"(b) = y"(b) - на концах промежутка [a,b] задаются значения 2-ой производной интерполируемой функции;
- краевые условия 3-го типа. s'(a) = s'(b), s"(a) = s"(b) - периодические краевые условия, которые применяются в случае, если интерполируемая функция является периодической с периодом T=b-a;
- краевые условия 4-го типа. Применяются, когда неизвестны ни первые, ни вторые производные в граничных точках и функция не является периодической. В этом случае недостающие два уравнения получаются из аппроксимации третьих производных на концах промежутка по четырем первым и четырем последним из заданных точек.
Недостающие уравнения выглядят следующим образом:
для краевых условий 1 типа
2h1σ1 + h1σ2 = ∆1 - y'(x1),
hn-1σn-1 + 2hn-1σn = y'(xn) - ∆n-1;
для краевых условий 2 типа
6σ1 = y"(x1), 6σn = y"(xn);
для краевых условий 3 типа
2(h1 + hn-1)σ1 + h1σ2 + hn-1σn-1 = ∆1 - ∆n-1,
σ1 = σn;
для краевых условий 4 типа
-h1σ1 + h1σ2 = h12∆1(3),
hn-1σn-1 - hn-1σn = hn-12∆n-3(3),
где ∆1(3) = [(∆3 - ∆2)/(x4 - x2) - (∆2 - ∆1)/(x3 - x1)]/(x4 - x1), ∆n-3(3) = [(∆n-1 - ∆n-2)/(xn - xn-2) - (∆n-2 - ∆n-3)/(xn-2 - xn-3)]/(xn - xn-3) - аппроксимации третьих производных.
Вызов программы
CALL SPL(X, Y, ib, D1, DN, B, C, D)
Параметры программы
X, Y, ib, D1, DN - входные параметры;
B, C, D - выходные параметры.
Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке;
Integer Ib - параметр, который задает тип граничных условий, 1<=Ib<=4;
Real D1, DN - значения этих параметров используются в граничных условиях 1 и 2 типа;
Real B(1:N), C(1:N), D(1:N) - массивы, содержащие коэффициенты сплайна;
При работе программы SPL используется программа LA30.
Вернуться к оглавлению Скачать SPL
SPF
Программа SPF вычисляет значение интерполяционного (SPL) или сглаживающего (SPS) сплайна в заданной точке X.
Для поиска прмежутка [xi..xi+1], в который попадает точка X, применяется двоичный поиск. Программа запоминает во внутренней переменной L номер промежутка, в котором вычислялся сплайн при предыдущем вызове. Если при последующем вызове точка X попадает в тот же промежуток, то новый промежуток можно не искать и сплайн вычисляется быстрее.
Вызов программы
CALL SPF(X, Y, B, C, D, XX, Id)
Параметры программы
X, Y, B, C, D, XX, Id - входные параметры;
Real Spf - возвращаемое значение.
Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции;
Real B(1:N), C(1:N), D(1:N) - массивы, содержащие значения коэффициентов сплайна;
Real XX - значение аргумента, в котором требуется вычислить значение функции;
Integer Id - параметр, задающий режим вычисления:
id=0 - программа вычисляет значение функции;
id=1 - программа вычисляет значение первой производной;
id=2 - программа вычисляет значение второй производной.
Пример
!Сплайн - интерполяция
program TestSPL
use NML
implicit none
integer, parameter:: n=60
real:: X(n), Y(n), B(n), C(n), D(n), XX, YY, YD, SY, SD
real, parameter:: pi=3.14159265359
integer:: ib=4
real:: H, D1, DN
integer i
!begin
H=2.0*pi/float(n-1)
do i=1, n
X(i)=H*float(i-1)
Y(i)=exp(-X(i))*sin(X(i))
end do
XX=1.2
YY=exp(-XX)*sin(XX)
YD=exp(-XX)*(cos(XX)-sin(XX))
if (ib==1) then
D1=1.0 ! D1=exp(0)*(cos(0)-sin(0))
DN=exp(-2.*pi) ! DN=exp(-2pi)*(cos(2pi)-sin(2pi))
else if (ib==2) then
D1=-2.0 ! D1=-2*exp(0)*cos(0);
DN=-2.*exp(-2.*pi) ! DN=-2*exp(-2pi)*cos(2pi)
else
D1=0.0; DN=0.0
end if
call SPL(X, Y, ib, D1, DN, B, C, D)
SY=SPF(X, Y, B, C, D, XX, 0)
SD=SPF(X, Y, B, C, D, XX, 1)
print 10, XX, YY, YD
print 11, SY, SD
10 format(/' X = ',F6.2,' YY =',F12.7,' YD =',F12.7)
11 format(/16x,'SY =',F12.7,' SD =',F12.7)
end program TestSPL
X = 1.20 YY = 0.2807248 YD = -0.1715847
SY = 0.2807250 SD = -0.1715749
Вернуться к оглавлению Скачать SPF
SPS
Программа SPS вычисляет коэффициенты кубического сглаживающего сплайна для функциональной зависимости yi = y(xi), i=1, 2, ..., n. Предполагается, что точные значения yi неизвестны, а заданы с некоторой погрешностью ρi, 0 ≤ ρi ≤ 1. В случае, если все ρi=0, то коэффициенты ai совпадают с yi и сглаживающий сплайн превращается в интерполяционный. Это означает, что чем точнее заданы величины yi, тем меньше должны быть соответствующие значения ρi. Если же требуется, чтобы сплайн прошел в точности через точку (xi,yi), то соответствующую ей погрешность ρi следует положить равной нулю.
Кубический сплайн представляет собой многочлен третьей степени на каждом из отрезков [xi, xi+1], i=1, 2, ..., n-1:
si(x) = ai + bi(x - xi) + ci(x - xi)2 + di(x - xi)3.
Для нахождения коэффициентов сплайна ai, bi, ci, di используются условия сшивания для самих многочленов и для их первой и второй производной в точках xi:
│s1(x2) = a1 + b1h1 + c1h12 + d1h13 = s2(x2) = a2
│ . . . . . . . . . . . . . . . . . . . .
│sn-1(xn) = an-1 + bn-1hn-1 + cn-1hn-12 + dn-1hn-13 = sn(xn) = an
│s1'(x2) = b1 + 2c1h1 + 3d1h12 = s2'(x2) = b2
│ . . . . . . . . . . . . . . . . . . . .
│sn-2'(xn-1) = bn-2 + 2cn-2hn-2 + 3dn-2hn-22 = sn-1'(xn-1) = bn-1
│s1"(x2) = 2c1 + 6d1h1 = s2"(x2) = 2c2
│ . . . . . . . . . . . . . . . . . . . .
│sn-2"(xn-1) = 2cn-2 + 6dn-2hn-2 = sn-1"(xn-1) = 2cn-1
где h1 = x2 - x1, ..., hn-1 = xn - xn-1. Дополнительно n уравнений получают из условия минимизации функционала
bn
J(s) = ∫(s"(x))2dx + ∑(1/ρi)(si - yi)2, x1=a, xn=b,
ai=1
которые имеют вид [Б5]:
│s1(x1) + ρ1s1'''(x1+) = y1
│s2(x2) + ρ2(s2'''(x2+) - s2'''(x2-)) = y2
│ . . . . . . . . . . . . . . . . . . . .
│sn-1(xn-1) + ρn-1(sn-1'''(xn-1+) - sn-1'''(xn-1-)) = yn-1
│sn(xn) - ρnsn'''(xn-) = yn
где si'''(xi+) и si'''(xi-) означают третью производную в точке xi справа и слева.
Введя обозначения
│c1 = 3σ1
│ . . .
│cn = 3σn
│∆1 = (y2 - y1)/h1
│ . . . . . . .
│∆n-1 = (yn - yn-1)/hn-1
и выполнив элементарные преобразования (см. прогр. SPL), получаем следующую систему уравнений:
│D2σ1 + A2σ2 + B2σ3 + C2σ4 = G2
│E3σ1 + D3σ2 + A3σ3 + B3σ4 + C3σ5 = G3
│ . . . . . . . . . . . .(*)
│En-2σn-4 + Dn-2σn-3 + An-2σn-2 + Bn-2σn-1 + Cn-2σn = Gn-2
│En-1σn-3 + Dn-1σn-2 + An-1σn-1 + Bn-1σn = Gn-2
где
Ai = 2(hi-1 + hi) + pi-1/(hi-1)2 + pi(1/hi-1 + 1/hi)2 + pi+1/(hi)2, i=2, ..., n-1
Bi = hi - (1/hi)[pi(1/hi-1 + 1/hi) + pi+1(1/hi + 1/hi+1)], i=2, ..., n-2
Bn-1 = hn-1 - (1/hn-1)[pn-1(1/hn-2 + 1/hn-1) + pn/hn-1]
Ci = pi+1/(hihi+1), i=2, ..., n-2
D2 = h1 - (1/h1)[p1/h1 + p2(1/h1 + 1/h2)]
Di = hi-1 - (1/hi-1)[pi-1(1/hi-2 + 1/hi-1) + pi(1/hi-1 + 1/hi)], i=3, ..., n-1
Ei = pi-1/(hi-2hi-1), i=3, ..., n-1
Gi = ∆i - ∆i-1, i=2, ..., n-1
Здесь pi = 6ρi.
Мы получили систему n-2 уравнений относительно n неизвестных. Недостающие два уравнения получаются из того или иного типа краевых условий.
Краевые условия 1-го типа. s'(a)=y'(a), s'(b)=y'(b) - на концах промежутка [a,b] задаются значения 1-ой производной сглаживаемой функции.
Краевые условия 2-го типа. s"(a)=y"(a), s"(b)=y"(b) - на концах промежутка [a,b] задаются значения 2-ой производной сглаживаемой функции.
Краевые условия 3-го типа. s(a)=s(b), s'(a)=s'(b), s"(a)=s"(b) - периодические краевые условия, которые применяются в случае, если сглаживаемая функция является периодической с периодом T=b-a.
Если значения ни первых, ни вторых производных на концах промежутка не известны и функция не является периодической, то следует положить s"(a)=s"(b)=0, в этом случае достигается эффект наилучшего сглаживания исходной таблично-заданной зависимости.
Недостающие уравнения выглядят следующим образом.
Для краевых условий 1-го и 2-го типа:
│A1σ1 + B1σ2 + + C1σ3 = G1
│Enσn-2 + Dnσn-1 + Anσn = Gn
Значения коэффициентов для краевых условий первого типа:
A1 = 2h1 + (p1 + p2)/(h1)2
B1 = h1 - p1/(h1)2 - p2(1/h1)(1/h1) + 1/h2)
C1 = p2/(h1h2)
En = pn-1/(hn-2hn-1)
Dn = hn-1 - pn-1(1/hn-1)(1/hn-2) + 1/hn-1) - pn/(hn-1)2
An = 2hn-1 + (pn-1 + pn)/(hn-1)2
G1 = ∆1 - y'(a), Gn = y'(b) - ∆n-1.
Значения коэффициентов для краевых условий второго типа:
A1=An=6, B1=C1=En=Dn=0, G1=y"(a), Gn=y"(b).
Краевые условия 3-го типа образуют сисему n-1 уравнений относительно n-1 неизвестного, поскольку в этом случае σ1=σn.
│A1σ1 + B1σ2 + C1σ3 + E1σn-1 + D1σn = G1
│D2σ1 + A2σ2 + B2σ3 + C2σ4 + E2σn = G2
│E3σ1 + D3σ2 + A3σ3 + B3σ4 + C3σ5 = G3
│ . . . . . . . . . . . . . . . .
│En-3σn-5 + Dn-3σn-4 + An-3σn-3 + Bn-3σn-2 + Cn-3σn-1 = Gn-3
│Cn-2σ1 + En-2σn-4 + Dn-2σn-3 + An-2σn-2 + Bn-2σn-1 + = Gn-2
│Bn-1σ1 + Cn-1σ2 + En-1σn-3 + Dn-1σn-2 + An-1σn-1 = Gn-1,
где
A1 = 2(hn-1 + h1) + pn-1/(hn-1)2 + p1(1/hn-1 + 1/h1)2 + p2/(h1)2,
An-1 = 2(hn-2 + hn-1) + pn-2/(hn-2)2 + pn-1(1/hn-2 + 1/hn-1)2 + p1/(hn-1)2,
B1 = h1 - (1/h1)[p1(1/hn-1 + 1/h1) + p2(1/h1 + 1/h2)],
Bn-1 = hn-1 - (1/hn-1)[pn-1(1/hn-2 + 1/hn-1) + p1(1/hn-1 + 1/h1)],
C1 = p2/(h1h2),
Cn-1 = p1/(hn-1h1),
Ai, Bi, Ci, i=2, ..., n-2 совпадает с (*),
D1 = Bn-1, D2 = B1, E1 = Cn-2, E2 = Cn-1,
Di+1 = Bi, Ei+1 = Ci-1, i=2, ..., n-2.
Коэффициенты сплайна находятся из уравнений:
│ai = yi - ρiHi
│bi = ∆i - hi(σi+1 + 2σi) - (1/hi)(pi+1Hi+1 - piHi)
│ci = 3σi
│di = (σi+1 - σi)/hi, i=1, ..., n-1,
an = yn - ρnHn,
где H1 = (σ2 - σ1)/h1, Hi = (σi+1 - σi)/hi - (σi - σi-1)/hi-1 i=2, ..., n-1, Hn = (σn - σn-1)/hn.
Вызов программы
CALL SPS(X, Y, P, ib, D1, DN, A, B, C, D)
Параметры программы
X, Y, P, ib, D1, DN - входные параметры;
A, B, C, D - выходные параметры.
Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке;
Real P(1:N) - массив, содержащий значения весовых коэффициентов, 0<=P(i)<=6. Значение P(i) должно быть равно 0, если значение ординаты Y(i) известно точно и должно равняться 6, если значение Y(i) нужно игнорировать;
Integer Ib - параметр, который задает тип граничных условий, 1<=Ib<=4;
Real D1, DN - значения этих параметров используются в граничных условиях 1 и 2 типа;
Real A(1:N), B(1:N), C(1:N), D(1:N) - массивы, содержащие коэффициенты сплайна;
При работе программы SPS используется программа LA31.
Пример
! Построение сглаживающего сплайна
! и графика сглаженной функции.
module SplPar
integer, parameter:: n=20
real:: X(n), Y(n), P(n), A(n), B(n), C(n), D(n)
integer:: ib
real:: D1, DN
end module SplPar
program TestSPS
$objcomment lib: "D:/Fortran/Bib_GR/Debug/Graphica"
use SplPar
use NML
implicit none
real, parameter:: pi=3.14159265359
integer:: seed(1)
real:: Q(n)
real:: H, AX=0.0, BX=2.0*pi
real, external:: FS
integer i, j
!begin
ib=2
P(:)=0.4
seed(1)=6
call random_seed(put=seed)
call random_number(Q)
H=(BX-AX)/float(n-1)
do i=1, n
X(i)=H*float(i-1)
Y(i)=sin(X(i))+0.2*(Q(i)-0.5)
end do
if (ib==1) then
D1=1.0
DN=1.0
else if (ib==2) then
D1=0.0
DN=0.0
else
D1=0.0; DN=0.0
end if
call SPS(X, Y, P, ib, D1, DN, A, B, C, D)
call Graphic(X(1), X(n), FS)
!график выведен в файл 'gf.txt'
end program TestSPS
real function FS(XX)
use SplPar
use NML
real, intent(in):: XX
FS=SPF(X, A, B, C, D, XX, 0)
return
end
Вернуться к оглавлению Скачать SPS
Hermite
Задано некоторое множество точек (вершин) P1 ... Pn на плоскости или в пространстве, n>1. Программа Hermite строит составную аппроксимирующую кривую Эрмита γ(1)ᑌ ... ᑌγ(n), которая является объединением элементарных кубических кривых γ(i), i = 1, ..., n-1. Элементарная кубическая кривая Эрмита определяется двумя вершинами Pi и Pi+1 и ненулевыми касательными векторами в этих точках Qi и Qi+1 (Рис 1).
Координаты точек кривой между вершинами Pi и Pi+1 можно найти из векторного уравнения
R(t) = (1 - 3t2 + 2t3) Pi + t2(3 - 2t) Pi+1 + t(1 - 2t + t2) Qi - t2(1 - t) Qi+1, 0 ≤ t ≤ 1
Это уравнение в координатной (матричной) форме может быть записано как
R(t) = GMT, 0 ≤ t ≤ 1
где (для точек P1 и P2)
Матрица M называется базисной матрицей кубической кривой Эрмита, а матрица G - её геометрической матрицей. Здесь x, y и z - пространственные координаты вершин P1 и P2, u, v, w - координаты касательных векторов Q1 и Q2.
Составная кривая Эрмита (Рис. 2) является объединением элементарных кривых Эрмита и определяется массивом вершин P1 ... Pn, n>1. Каждая i-ая кривая описывается уравнением
R(i)(t) = (Pi Pi+1 Qi Qi+1)MT, 0 ≤ t ≤ 1
где M - базисная матрица, а промежуточные векторы Q2 ... Qn-1 можно найти из условий равенства вторых производных в базисных точках P2 ... Pn-1, которые сводятся к матричному уравнению
Свойства составной кривой Эрмита описаны в [Б5].
Программная реализация
Для начала вычислений нужно задать вершины, через которые должна пройти кривая, и векторы касательных в конечных вершинах. Программа HermiteTangents вычисляет значения касательных в промежуточных вершинах. Программа Hermite вычисляет координаты заданной точки на построенной кривой. Точка на кривой задаётся параметром t, который должен лежать в пределах от 0 до n-1. Целая часть параметра определяет номер i элементарной кривой Эрмита, а дробная часть задаёт положение точки на этой кривой.
[Б5]
Вызов программы HermiteTangents
call HermiteTangents(XP, xq1, xqN, XQ)
Параметры программы
XP(1:n), xq1, xqN - входные параметры;
XQ(1:n) - выходной параметр.
Real XP(1:n) - массив, содержащий x-, y- или z-координаты точек Pi, i=1, ..., n;
Real xq1, xqN - x-, y- или z-координаты касательных векторов в точках P1 и Pn:
Real XQ(1:n) - массив, содержащий x-, y- или z-координаты касательных векторов Qi в точках Pi, i=1, ..., n;
integer n - задаёт количество точек P и размерность массивов XP и XQ;
Вызов программы Hermite
Hermite(XP, XQ, t)
Программа Hermite должна быть вызвана после вызова программы HermiteTangents
Параметры программы
XP(1:n), XQ(1:n), t - входные параметры;
Real Hermite - возвращаемое значение;
Real XP(1:n) - массив, содержащий x-, y- или z-координаты точек Pi, i=1, ..., n;
Real XQ(1:n) - массив, содержащий x-, y- или z-координаты касательных векторов Qi в точках Pi, i=1, ..., n (вычислены с помощью программы HermiteTangents);
integer n - задаёт количество точек P и размерность массивов XP и XQ;
Real t - значение параметра, в котором нужно найти x-, y- или z-координату составной кривой Эрмита, 0 ≤ t ≤ n-1
Примеры
program TestHermit
! Построение аппроксимирующей кривой Эрмита
use GRAPHICA
use NML
implicit none
integer, parameter:: n=8, div=32
real XP(n), YP(n), XQ(n), YQ(n), xq1, yq1, xqN, yqN
real T(div), XT(div), YT(div), dt
integer i
data XP/0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.80, 0.95/
data YP/0.90, 0.50, 0.40, 0.50, 0.25, 0.10, 0.50, 0.90/
data xq1/0.1/, yq1/-0.02/, xqN/0.1/, yqN/0.01/
! Вычисление X и Y координат касательных во внутренних точках
call HermiteTangents(XP, xq1, xqN, XQ)
call HermiteTangents(YP, yq1, yqN, YQ)
! Заполнение массива с параметром
dt=float(n-1)/float(div)
do i=1,div
T(i)=(i-1)*dt
enddo
! Построение аппроксимирующей кривой
do i=1,div
XT(i)=Hermite(XP, XQ, T(i))
YT(i)=Hermite(YP, YQ, T(i))
enddo
! Вывод результатов вычислений
print 10
print 11, (T(i), XT(i), YT(i), i=1, div)
call Graphar(XT, YT)
10 format(8X,'t',10X,'x',14X,'y')
11 format(2X,F10.7,2E15.7)
end program TestHermit
График построен программой Advanced Grapher
program TestHermit2
! Построение кривой Эрмита c самопересечением
use GRAPHICA
use NML
implicit none
integer, parameter:: n=8, div=64
real XP(n), YP(n), XQ(n), YQ(n), xq1, yq1, xqN, yqN
real T(div), XT(div), YT(div), dt
integer i
data XP/0.10, 0.30, 0.50, 0.70, 0.95, 0.70, 0.40, 0.10/
data YP/0.30, 0.50, 0.65, 0.55, 0.45, 0.30, 0.30, 0.45/
data xq1/0.05/, yq1/0.02/, xqN/-0.05/, yqN/-0.04/
call HermiteTangents(XP, xq1, xqN, XQ)
call HermiteTangents(YP, yq1, yqN, YQ)
dt=float(n-1)/float(div)
do i=1,div
T(i)=(i-1)*dt
XT(i)=Hermite(XP, XQ, T(i))
YT(i)=Hermite(YP, YQ, T(i))
enddo
! Вывод результатов вычислений
print 10
print 11, (T(i), XT(i), YT(i), i=1, div)
call Graphar(XT, YT)
10 format(8X,'t',10X,'x',14X,'y')
11 format(2X,F10.7,2E15.7)
end program TestHermit2
Вернуться к оглавлению Скачать Hermite
ORT
Программа ORT осуществляет разложение таблично-заданной зависимости yi=y(xi), i=1, ..., n в ряд по системе многочленов, ортогональных на множестве точек {xi}:
m
y(x) = ∑ αkuk(x),
k=0
где αk - коэффициенты этого разложения, а uk(x) - ортогональные многочлены последовательно возрастающих степеней:
k
uk(x) = ∑ βj(k)xj k=0, ..., m.
j=0
Предполагается, что значения yi заданы с некоторой неопределённостью (с весовыми коэффициентами) pi, 0 < pi ≤ 1. Значение pi=1 означает, что соответствующее значение yi задано точно, и наоборот, чем ближе pi к нулю, тем менее учитывается данная точка при аппроксимации. Весовые коэффициенты легко вычислить, если известна погрешность ∆i величин yi, т.е. yi ± ∆i. Тогда весовые коэффициенты будут равны
pi = 1 - (∆i/∆max),
где ∆max - максимальное значение погрешности.
Система ортогональных многочленов строится по рекуррентному соотношению
A*uk+1(x) = (x-A)uk(x) - Buk-1(x), u0=0, u1=1/SQRT(∑pi),
где A=∑pixiuk2(xi), B=∑pixiuk(xi)uk-1(xi). A* находится из условия нормировки
n
∑ piuk+12(xi) = 1.
i=1
Коэффициент разложения αk определяется из формулы
n
αk = ∑ piyiuk(xi).
i=1
Программа последовательно увеличивает количество ортогональных многочленов, участвующих в аппроксимации заданной дискретной зависимости y(xi), пока степень старшего многочлена не достигнет заданного значения m (m≥0) или пока из-за накопления вычислительной погрешности среднеквадратичное отклонение
не перестанет уменьшаться. Достигнутая программой степень старшего многочлена возвращается в переменной Order.
Для упрощения дальнейших расчетов коэффициенты всех многочленов при одинаковых степенях x суммируются и результатом выполнения программы является массив коэффициентов c0, c1, ..., COrder. Значение функции в точке x (x1 ≤ x ≤ xn) можно вычислить по формуле
Order
y(x) = ∑ ckxk.
k=0
Краткое описание алгоритма построения многочленов, ортогональных на множестве точек, можно найти в [А4]. Данный метод более устойчив к накоплению ошибок вычислений, чем метод наименьших квадратов.
Вызов программы
CALL ORT(X, Y, P, C, M, Order, Sigma)
Параметры программы
X, Y, P, M - входные параметры;
C, Order, Sigma - выходные параметры.
Real X(1:N), Y(1:N) - массивы, содержащие значения аргумента и функции. Массив X должен быть упорядочен в возрастающем порядке;
Real P(1:N) - массив, содержащий значения весовых коэффициентов, 1>=P(i)>0. Значение P(i) должно быть равно 1, если значение ординаты Y(i) известно точно и стремиться к 0, если значение Y(i) можно игнорировать;
Real C(0:M) - массив, содержащий коэффициенты многочлена, аппроксимирующего исходную зависимость;
Integer M - наибольшая степень многочлена, которую разрешается достигнуть программе;
Integer Order - наибольшая степень многочлена, которую удалось достигнуть программе, Order ≤ M;
Real Sigma - среднеквадратичное отклонение, достигнутое при аппроксимации.
Пример
! Аппроксимация функции многочленами,
! ортогональными на множестве точек.
module OrtPar
integer, parameter:: M=20
real:: C(0:M)
integer:: Order
end module OrtPar
program TestORT
$objcomment lib: "D:/Fortran/Bib_GR/Debug/Graphica"
use OrtPar
use NML
implicit none
integer, parameter:: n=40
real:: X(n), Y(n), P(n)
data P/n*1.0/
real, external:: OrtPol
real, parameter:: pi=3.1415926
real:: H, XX, YY, YYY, Sigma
integer:: i
!begin
H=2.0*pi/float(n-1)
do i=1, n
X(i)=H*float(i-1)
Y(i)=sin(X(i))
end do
XX=0.5
YY=sin(XX)
call ORT(X, Y, P, C, M, Order, Sigma)
YYY=OrtPol(XX)
print 10, XX, YY, YYY, Order, Sigma
10 format(/' X = ',F6.2,' YY =',F12.7,' YYY =',F12.7 &
/' Order =',I3,' Sigma =',E15.7/)
call Graphic(X(1), X(n), OrtPol)
print 20, (i, C(i), i=0, Order)
20 format(8x,I2,4x,E14.6)
end program TestORT
real function OrtPol(X)
use OrtPar
real, intent(in):: X
real:: S
integer i
S=0.0
do i=Order, 1, -1
S=(S+C(i))*X
end do
OrtPol=S+C(0)
return
end function OrtPol
X = 0.50 YY = 0.4794255 YYY = 0.4794317
Order = 9 Sigma = 0.2819624E-04
0 0.361258E-05
1 0.999756E+00
2 0.152848E-02
3 -0.170313E+00
4 0.445399E-02
5 0.517216E-02
6 0.138082E-02
7 -0.574159E-03
8 0.610403E-04
9 -0.215886E-05
Вернуться к оглавлению Скачать ORT
