Функции, заданные своими значениями, сведенными в таблицу, часто встречаются на практике. Математики называют их решетчатыми. С точки зрения программирования эти объекты существенно отличаются от зависимостей, заданных математическими выражениями.
Функции в форме таблиц - это результаты численного решения дифференциальных или интегральных уравнений, либо эмпирические данные, результаты эксперимента. Их включение во многие вычислительные алгоритмы невозможно без предварительного преобразования, изменяющего тип области определения. (Таблица определена на поле натуральных чисел, тогда как методы математического анализа требуют области определения функции на поле действительных чисел.) Такое преобразование может быть осуществлено двумя операциями: интерполяцией или аппроксимацией.
При интерполяции кривая проходит точно через указанные точки таблицы, при аппроксимации - вблизи точек, но так, что суммарное отклонение в некотором смысле минимально. Аппроксимацию применяют тогда, когда необходимо получить общее выражение для большого множества точек. В библиотеке представлены: интерполяция степенными многочленами, сплайнами и рядами Фурье; алгебраическая и гармоническая аппроксимация по методу наименьших квадратов.
Сплайн - интерполяционная формула, обеспечивающая непрерывность как самой функции, так и нескольких ее производных. При последующем дифференцировании сплайн дает меньшую погрешность, нежели интерполяционный многочлен того же порядка.
Интерполяция рядами Фурье применяется к периодическим функциям, для которых обеспечивает минимальную погрешность на всей оси X.
Если необходимо заменить сложное выражение функциональной зависимости более простым, но менее точным, применяют аппроксимацию. Так как узлы в этом случае можно выбирать по своему усмотрению, используют метод Чебышева.
Методы численного интегрирования сгруппированы в разделе "Функции в форме выражений". Ниже даны лишь рекомендации по их адаптации к интегрированию таблиц.
1. В общем случае табличная функция предварительно интерполируется степенным многочленом, после чего к ней может быть применена любая формула интегрирования. Целесообразно выбирать порядок интерполяции, совпадающий с высшей степенью многочлена, для которого интегрирование проходит без погрешности.
2. Если границы интервала интегрирования совпадают с узлами таблицы, можно подобрать параметры интегратора так, чтобы для вычисления требовались только точки-узлы таблицы. На эту операцию ориентированы формулы Ньютона - Котеса (программы 86, 87), отличающиеся постоянным шагом интегрирования.
3. Часто требуется получить таблицу интеграла с переменным верхним пределом, в которой узлы табуляции совпадали бы с узлами исходной функции. В этом случае рационально интегрировать функцию по варианту 2, затем интерполировать первообразную, вычисляя недостающие значения в промежуточных узлах.
4. Если в отрезок интегрирования попадает не более 5:7 узлов таблицы, можно вычислить коэффициенты интерполяционного многочлена, проходящего через все узлы (программа 97), а интеграл от него взять аналитически (программа 29). При большем числе узлов наблюдается резкий рост вычислительной погрешности.
Численному дифференцированию с равной легкостью поддаются как функции-таблицы, так и функции-выражения. Однако не следует дифференцировать таблицу, содержащую шум, так как при дифференцировании он усиливается. Зашумленные таблицы предварительно аппроксимируют по методу наименьших квадратов.
Программа 92. Вычисление трех производных по семи точкам. В узлах регулярной решетки xi=x0+is, где s - шаг, x0 - центральная точка, i=-3;3, программа вычисляет значения трех первых производных. Формулы точны для многочлена 6-го порядка.
y'=45{(y1-y-1)-9(y2-y-2)+(y3-y-3)}/60/s;
y"=490y0+270(y1+y-1)-27(y2+y-2)+2(y3+y-3)/180/s2;
y'''=-13(y1-y-1)+8(y2-y-2)-(y3-y-3)/8/s3.
После расчета производной в одной точке центральный узел сдвигается на s вправо и процесс повторяется.
10 INPUT "s",A,"x0",B:B=B-4*A:FOR
C=1 TO 7:B=B+A:GOSUB 60
20 NEXT C:PRINT "x";B-3*A,"y'";(45*(H-F)-9*(I-E)+J-D)/60/A
30 PRINT "y2";(270*(H+F)-490*G-27*(I+E)+2*(J+D))/180/A/A
40 PRINT "y3";(13*(F-H)-8*(E-I)+D-J)/8/A-3
50 B=B+A:GOSUB 60:FOR C=1 TO 7:C(C)=D(C):GOTO 20
60 PRINT "y(";B;")";:INPUT C(C):RETURN
Размер: 234, K
Пример: y=x6; s=1; x0=0.
Модификация программы:
10 INPUT "s",A,"x0",B:B=B-4*A:FOR
C=1 TO 7:B=B+A:GOSUB 60
20 NEXT C:PRINT "y";G,"y1";(45*(H-F)-9*(I-E)+J-D)/60/A
30 PRINT "y2";(270*(H+F)-490*G-27*(I+E)+2*(J+D))/180/A/A
40 PRINT "y3";(13*(F-H)-8*(E-I)+D-J)/8/A-3,"X";B-A-A
50 B=B+A:GOSUB 60:FOR C=1 TO 7:C(C)=D(C):GOTO 20
60 C(C)=B-6:RETURN
Ответ: y=0; y'=0; y"=0; y'''=0; x=1; y=1; y'=6; y"=30; y'''=120;...
Программа 93. Вычисление первой производной с интерполяцией по пяти точкам. Необходимо задать шаг решетки s и значения функции в узлах интерполяции yi=f(x0+is), i=-2;2. Теперь для любого p=(x-x0)/s производная будет вычисляться по формуле:
f'(x0+ps)={(y0(2p2-5)/2 - (y1+y-1)(2p2-4)/3+(y2+y-2)(2p2-1)/12)p+(y1-y-1)(3p2-4)/6-(y2-y-2)(2p2-1)/12}/s.
Формула точна для многочлена 4-го порядка.
10 INPUT "s",A:FOR B=-2 TO 2:PRINT
B;:INPUT G(B):NEXT B
20 INPUT "p",B:C=B*B
30 D=(C-.5)/6*(E+I)-(C-2)/1.5*(F+H)+(C-2.5)*G
40 D=D*B+(3*C-1)/12*(I-E)-(3*C-4)/6*(H-F):PRINT D/A:GOTO 20
Размер: 140, I
Пример: s=0,5; y-2=1; y-1=1/16; y0=0; y1=1/16; y2=1; x=0,5.
Ответ: f'(0,25)=0,0625.
Интерполяционная формула Лагранжа:
f(x)=е(i)yiХ(j)(x-xj)/(xi-xj); i=0...n, j=0...n.
где {xi, yi} - узлы интерполяции.
10INPUT A:FOR B=0 TO A:PRINT B;:INPUT
G(B+B),"y",H(B+B):NEXT B
20 INPUT C:D=0:FOR B=0 TO A:E=1:FOR F=0 TO A
30 IF B¦F;E=E*(C-G(F+F))/(G(B+B)-G(F+F))
40 NEXT F:D=D+E*H(B+B):NEXT B:PRINT D:GOTO 20
Размер: 131, H(2*n)
Пример: n=3; x0=-1; y0=-1; x1=0; y1=0; x2=1; y2=1; x3=2; y3=8; x=1,5.
Ответ: f(1,5)=3,375.
Программа 95. Случай регулярной решетки. Справедлива замена:
p=(x-x0)/s, где s - шаг решетки.
10 INPUT A,"x",D,"s",G:FOR B=0 TO
A:PRINT B;:INPUT H(B):NEXT B
20 INPUT C:C=(C-D)/G:D=0:FOR B=0 TO A:E=1
30 FOR F=0 TO A:IF B¦F;E=E*(C-F)/(B-F)
40 NEXT F:D=D+E*H(B):NEXT B:PRINT D:GOTO 20
Размер: 123, H(n)
Пример: n=3; x0=-1; s=1; y0=-1; y1=0; y2=1; y3=8; x=1,5.
Ответ: f(1,5)=3,375.
Интерполяционная формула Ньютона:
f(x)=r0+(x-x0){r1+(x-x1){r2+ ... (x-xn-1) rn} ... };
где ri - разделенная разность порядка i.
f(x) = е(i)aixi, i=0...n
где ai - неизвестные коэффициенты.
10 INPUT A:FOR B=0 TO A:PRINT B;:INPUT
E(B+A),"y",D(B):NEXT B
20 FOR B=1 TO A:FOR C=A TO BSTEP -1
30 D(C)=(D(C)-C(C))/(E(C+A)-E(C+A-B)):NEXT C:NEXT B
40 B=B-1:IF B=0;FOR B=0 TO A:PRINT B;D(B):NEXT B:GOTO 10
50 FOR C=B TO A:C(C)=C(C)-D(C)*D(B+A):NEXT C:GOTO 40
Размер: 171, E(2*n)
Пример: n=4; x0=0,2; y0=5,9376; x1=0,3; y1=6,5321; x2=0,4; y2=7,2336; x3=0,5; y3=8,0625; x4=0,6; y4=9,0416.
Ответ: a0=5; a1=4; a2=3; a3=2; a4=1.
Программа 97. Случай регулярной решетки. Справедлива замена:
p=(x-x0)/s, где s - шаг решетки.
10 INPUT "n",A,"x0",D,"s",E:FOR
B=0 TO A:PRINT B*E+D;
20 INPUT F(B):NEXT B:FOR B=1 TO A:FOR C=A TO BSTEP -1
30 F(C)=(F(C)-E(C))/B/E:NEXT C:NEXT B
40 B=B-1:IF B=0;FOR B=0 TO A:PRINT B;F(B):NEXT B:GOTO 10
50 FOR C=B TO A:E(C)=E(C)-F(C)*(D+B*E-E):NEXT C:GOTO 40
Размер: 170, F(n)
Пример: n=4; x0=0,2; s=0,1; y0=5,9376; y1=6,5321; y2=7,2336; y3=8,0625; y4=9,0416.
Ответ: f(x)=5+4x+3x2+2x3+x4.
f(t)= е(i)ciTi(t), i=0...n
где Ti(t) - многочлен Чебышева порядка i. Узлы интерполяции в относительном масштабе:
ti=cos(p(2i+1)/(n+1)/2), i=0...n.
Относительный масштаб соответствует интервалу интерполяции [-1;1]. При произвольном интервале [a,b] узлы определяются по формуле:
xi=(-a+b+ti(b-a))/2.
Вычислив значения функции в узлах yi=f(xi), можно определить коэффициенты Чебышева
f(t)= 2/(n+1)е(k)yicos(ip(2k+1)/(n+1)/2), k=0...n
Программа 98. Аппроксимация методом Чебышева. Значение интерполяционного многочлена для произвольной точки x внутри интервала [a,b] определяется алгоритмом:
t=(2x-b-a)/(b-a);
bk=ck-bk+2+2tbk+1,
k=n...1;
f(x)=c0/2-b2-tb1.
10 VAC:MODE 5:INPUT "n",C,"a",A,"b",B:FOR
D=0 TO C
20 F=(A+B)/2+(B-A)/2*COS (p*(D+D+1)/(C+C+2)):GOSUB
70
30 FOR F=0 TO C:I(F)=I(F)+2/(C+1)*E*COS (p*(D+.5)*F/(C+1))
40 NEXT F:NEXT D:I=I/2
50 INPUT "x",E:E=(E+E-B-A)/(B-A):F=0:G=0:FOR D=1 TO C
60 H=J(C-D)-F+2*G*E:F=G:G=H:NEXT D:PRINT I-F+G*E:GOTO 50
70 E= ...f{F}...
:RETURN
Размер: 234+, I(n)
Пример: f(x)=5-6x+37x2-12x3+2x4; n=4; a=-1; b=1; x=0,5.
В программе: 70 E=(((2*F-12)*F+37)*F-6)*F+5:RETURN
Ответ: f(0,5)=9,875.
Программа 99. Вычисление коэффициентов аппроксимирующего многочлена. Программа нормализует интервал аппроксимации, вычисляет коэффициенты Чебышева, затем пересчитывает их в коэффициенты обычного степенного многочлена, для которого проводит обратную замену переменных.
10 VAC:MODE 5:INPUT "n",C,"a",A,"b",B:FOR
D=0 TO C
20 F=(A+B)/2+(B-A)/2*COS (p*(D+D+1)/(C+C+2)):GOSUB
80
30 FOR F=0 TO C:I(F)=I(F)+2/(C+1)*E*COS (p*(D+.5)*F/(C+1))
40 NEXT F:NEXT D:I=I/2:FOR E=2 TO C:FOR F=C TO E STEP -1
50 G(F)=G(F)-I(F):I(F)=2*I(F):NEXT F:NEXT E:D=2/(B-A)
60 B=(B+A)/(A-B):FOR E=0 TO C:G=0:H=1:FOR F=E TO C:G=G+H*I(F)
70 H=H*B*(F+1)/(F+1-E):NEXT F:PRINT E;G*D-E:NEXT
E:END
80 E= ...f{F}...
:RETURN
Размер: 298+, I(C)
Пример: f(x)=5-6x+37x2-12x3+2x4; n=4; a=-1; b=1.
В программе: 80 E=(((2*F-12)*F+37)*F-6)*F+5:RETURN
Ответ: a0=5; a1=-6; a2=37; a3=-12; a4=2.
f(x)=(xi+1-x)2(2(x-xi)+s)yi/s3+(x-xi)2(2(xi+1-x)+s)yi+1/s3+(xi+1-x)2(x-xi)y'i/s2+(x-xi)2(xi+1-x)y'i+1/s2;
i=0..n-1; xi≤x<xi+1; s=xi+1-xi.
Программа 100. Интерполяция локальным сплайном. Локальный кубический сплайн требует задания значений функции и ее производной в узлах регулярной решетки и сохраняет первую производную непрерывной.
Если значения производной неизвестны, они могут быть получены численным дифференцированием:
y'i=(yi+1-yi-1)/2/s.
Но тогда интервал интерполяции сужается до [x0+s, x0+(n-1)s], так как производная на краях не определена.
10 INPUT "n",A,"x0",B,"s",C
20 FOR D=0 TO A:PRINT B+D*C;:INPUT G(D):NEXT D
30 FOR D=1 TO A-1:G(A+D)=(H(D)-F(D))/2/C:NEXT D
40 INPUT "x",E:E=(E-B)/C:D=INT E:E=FRAC E
50 F=(1-E)*(1-E)*((E+E+1)*G(D)+E*C*G(A+D))
60 F=F+E*E*((3-E-E)*H(D)-(1-E)*C*H(A+D)):PRINT F:GOTO 40
Размер: 202,F(2*n)
Пример: Значения производной в узлах известны.
n=1; x0=1; s=0,1; y0=10; y'0=10; y1=11,051; y'1=11,03; x=1,05.
Модификация программы:
10 INPUT "n",A,"x0",B,"s",C:FOR
D=0 TO A:PRINT B+D*C;
20 INPUT G(D),"y'",H(A+D):NEXT D
30 INPUT "x",E:E=(E-B)/C:D=INT E:E=FRAC E
40 F=(1-E)-2*((E+E+1)*G(D)+E*C*H(A+D))
50 F=F+E*E*((3-E-E)*H(D)-(1-E)*C*I(A+D)):PRINT F:GOTO 30
Ответ: f(1,05)=10,512625.
f(x)= ((1-p)3y"i+p3y"i+1)s2/6+(yi-y"is2/6)(1-p)+(yi+1-y"i+1s2/6)p;
где s - шаг решетки; yi - значения интерполируемой функции в точках xi=x0+is; i=Int[(x-x0)/s]; p=Frac[(x-x0)/s]; y"i - неизвестные значения второй производной сплайн-функции; определяются из системы уравнений:
y"i+0,25(y"i-1+y"i+1)=1,5(yi-1-2yi+yi+1)/s2; i=1...n;
y"0=y"n, y"n+1=y"1- граничные условия для периодического сплайна.
Система решается методом Зейделя.
Глобальный кубический сплайн обеспечивает непрерывность второй производной.
Программа 101. Вычисление вторых производных сплайна
10 INPUT "n",A,"s",C:FOR B=1 TO
A:PRINT B;:INPUT C(B):NEXT B
20 FOR B=0 TO A-1
30 D(A+B)=D(A*FRAC ((A+B-1)/A))-2*D(B)+D(A*FRAC ((B+1)/A))
40 D(A+B)=D(A+B)*6/C/C:NEXT B
50 $=CHR 0:FOR B=0 TO A-1
60 C=(D(A+B)-D(A*FRAC ((A+B-1)/A))-D(A*FRAC ((B+1)/A)))/4
70 IF D(B);D(B)=C:$=CHR 1
80 NEXT B:IF $=CHR 1 THEN 50
90 FOR B=1 TO A:PRINT B;C(B):NEXT B
Размер: 242, C(2*n)
Пример: n=4; s=p/2; y1=1; y2=0; y3=-1; y4=0.
Ответ: y"1=-1,21585419; y"2=0; y"3=1,21585419; y"4=0.
Программа 102. Интерполяция глобальным сплайном
10 INPUT "n",A,"x0",D,"s",C:FOR
B=0 TO A-1:PRINT B;
20 INPUT H(A+B):NEXT B:FOR B=0 TO A-1:GOSUB 110
30 H(A+A+B)=H(A+F)-2*H(A+B)+H(A+G):NEXT B
40 $=CHR 0:FOR B=0 TO A-1:GOSUB 110:E=(H(A+A+B)-H(F)-H(G))/4
50 IF H(B);H(B)=E:$=CHR 1
60 NEXT B:IF $=CHR 1 THEN 40
70 INPUT "x",E:E=A*FRAC ((E-D)/C/A):IF E<0;E=E+A
80 B=INT E:E=FRAC E:GOSUB 120:F=1-E
90 F=H(B)*F-3+H(G)*E-3+(H(A+B)-H(B))*F+(H(A+G)-H(G))*E
100 PRINT F:GOTO 70
110 F=A*FRAC ((A+B-1)/A)
120 G=A*FRAC ((B+1)/A):RETURN
Размер: 342,G(3*n)
Пример: n=4; x0=0, s=90; y0=1; y1=0; y2=-1; y3=0; x=60.
Ответ: f(30)=0,481481482.