330 LETX=1>K!LETV=1/V!60T0 460

340 LETX=1X:LETV=L0G<V):60TO 406

350 LETX=1/X!LETV=LGT<V>!&0T0 400

360 LETX=X-H:60T0 400

390 REMПИНЕйИАЯ РЕГРЕССИЯ

400 LETXt=Xl+X:LETVl=Vl+V;LETX-2=X-2+X*X

410 LErV2=V2+Vi«V!LETP=P+Xi«V!NEXT I

420 LETU= С X1 *V 1 - J«P > У к X1 *X 1 -jj«X2 )

430 LETW-<V1-ижХl)yj

44Й LETR=cP--Xl*VlJ)SGR<:c:X2-Xli«Kt/J>*i;V£-Vl*VlJ-> 450 GOTO 60

500 REMСВЕДЕНИЕ ЛИНЕЙНОЙ РЕГРЕССИИ К НЕЛИНЕЙНОЙ 510 GOTO 90 520 GOTO 90 530 GOTO 90 540 GOTO 90

550 LETA=10-bl5LETB=10"U!GOTO 90 ..

560 LETft=EXPaO!G0T0 98

570 LEtA=10"WsGOTO 98 • .

580 GOTO 98 59© GOTO 570 600 GOTO 90 618 GOTO 90

628 LETA=1.U!LETB=WA)EG0T0 98 638 LETA=t-W!LETB=U/W!GOTO 98 640 GOTO 568 658 GOTO 578 660 GOTO 98

780 REMВЫЧИСЛЕНИЕ V-IX) ПО ЗАДАННОМУ X 718 LETV=A+B-*Xi!GOTO 138

728 LETVI-CA+BSKXJSGOTO 138 • . ; • s

730 LETV=A+B/-X3G0T0 138 • r; :

748 LETV=X-(A+B3«X)!&0T0 138

758 LETV=A*B"X!GOTO 138

768 LETV=A*EXP<E*X)!60T0 138

770 .LETV=A»ie4U*X>s60T0 138

780 LETV=1/<:A+B*EXP<-X>>!60T0 130

790 LETV=A»X"B!60T0 138

888 LETV=A+B»LGT<X>!80Tu 138

818 LETV=A+B*L06«>!G0T0 138

820 LETV=A/:B+X)sGOTO 138

838 LETV=A!iiX/<B+XV.60T0 138

840 LETV=A*EXP<:B-X)!G0T0 138

858 LETV=A«10"":e/X>!&OTO 130

868 LETV=A+B*X"H!60T0 130

900 PRINTТАБЛИЦА КОДОВ И ЭМПИРИЧЕСКИХ ФУНКЦИЙ <г<Ю 910 PRINT1. V=A+b»X 2. V=1/<A+B*X>

920 PRINTS. V=A+B/X 4. У=Х-<А+ВжХ)

-938 PRINTS. V=A*B"X 6. V=A*EXP<b*X)

948 PRINT7. VASKie-iBiitX) 8. ¥=1/<А+ВжЕХР<-Х>) 958 PRINT9. V=A!«Xb 18.V=A+B*LGT<:x)

960 PRINT11.V=A+BsifLOG<X> 12.V=A-<:b+X)

978 PRINT13.V=A*X/<B+X) 14.V=A-*EXP<B/-X>

980 PRINT15.V=ft*10":B/X) 16.V=A+B*<;XN>

998 INPUTЗАДАЙТЕ КОД РЕГРЕССИИ К=К 1000 GOTO 65:END

Контрольные примеры (везде xi == 1, 5. i/i=2, у2 = 4. j;3==ll, У4 = 31 (/4 = 1,0117,

л:2 = 2. л:з = 3,5, л:4 = 5, номер совпадает с кодом к, В = 1,98157, i? = 0,99999, у (2) =3,97262, у (4) =

а результат дан в скобках). =!5-,59902).

1. 1/1=0,4,2 = 0,9,1/3 = 2,5,1/4 = 4,2 ( = 1,13537, 6. у,=7,4, № = 53, 1/з = 1040, у4 = 20300 {А = В=1,0558, R = 0,99908, у (2) =0,97619, у (4)= =1,0168, В=1,9800, « = 0,999999, у {2) = =-3,0878). =53,34312, у (4) ==-2798,46803).

2. у=0.1, S2 = 0,05, 1/3 = 0,027, i/, = 0,02 {А= 7. i/,=3, i/2 = 5, уз=\0. i/4 = 20 ( = 1,90498, = 0,11086, В = 10,13857, « = 0,99845. i/ (3) = В = 0,20501, « = 0,99979. i/ (1) =3,05418, i/ (3) = = 0,032758, у (1,5)=0.06528). =7.85066).

3.4,1 = 10.3, 1/2 = 5, 4/3 = 3, 4,4 = 2,4 (A= 8. у,=2,1, 1/2 = 4,2, 1/з = 7,7, j/4 = 9,3 (A =

-=0,20763. B= 10,0062, « = 0,99880, у (I) = =0,099793, В =1,02288, « = 0,099999. J/(2) =

= 10,21384, у (3) =3,543034). =4,19771, у (4) =8,43686).

4. i/i = 0,17, 1/2 = 0,18, 1/3 = 0,19, (/4 = 0,192 9. i/, = 1, 1/2 = 4, 1/3= 12, у4 = 25 ( =0,99935,

(=0,93817. B = 5,01769, « = 0,99993, J/(1)= fi= 1.99493, « = 0,99998. 4, (2) =3,98338, 1/ (4) =

= 0,16790,/у (3)=0,18760). =15,87764).

Ш. j;, = l, jrj=!,6, jr3 = 2,i, 6Г4 = 2,4 (/4 = =0,99924, В = 2,00965, «=0,99996, y(2)=-= 1,60421. у (4) =2,20917).

И. y, = l, 1/2=2,4, бГз=3,5, i/,=4,2 (/4 = = i .00831, fi = i ,98764, R = 0,99997, у (2)= 2,38604. у (4) =3,76377).

12. {,,=3,3, {/e==2.5, № = i,8, ;У4 = 1,4 (/4 = =9,70274. B= 1,91056, R = 0,99986. у (2) = = 2,48116, у (4) = l,64i59).

13. 5г, = 3,33, j;j = 5, 1,3=6,36, y, = 7,li. {A = = 10,00373, B = 2,00375, i?==0,999999, y{2) = = 4,99719, у (4) =6,66499).

14. 1,1 = 74, s,! = 27, t,3=17, {,,= 15 (Л = =9,80145, B = 2.01995, R = 0,99963, j; (2) = = 26,91020. у (4) = 16,24066).

15. /у, = 500. у2=50. 1,3=19, ;У4=13 (/4 = = 5,16152. В.= 1,98393, R=i0,99998, у{2) = = 50,66928, I, (4) = 16,1719).

16. j;,=2,!, 1,2 = 2,28, г/з = 2,65, 1,4 = 3.12 {А = = 1,99742, B = 0,i0022, /? = 0,99998. jr (2) = =2,28001, у (4) =2,79915).

В этих примерах результаты указаны с точностью до о-6 знаков после запитой. Время выполнения каждого контрольного примера составляет несколько секунд (при четырех парах xi, yi). С ростом числа пар xi, у,- оно пропорционально растет.

§ П5.11. Сплайн-аппроксимация, интерполяция и экстраполяция

Полиномиальная интерполяция и аппроксимация не обеспечивает непрерывность производных функции у {х) -л может давать значительные по-грешао.;тй в промежутках Ь5ежду узлами. Кроме того, она плохо приспособлена дли экстраполяции и, как правй.!!о, не обеспечивает правильное асимптотическое поведение у {х) яри изменении аргумента х за пределами интервала интерполяции. Нередко с увеличением числа узлов погрешность такой интерполяции не только не уменьшается, но и начинает расти.

От зтйх недостатков свободна аппроксимация и интерполяции с помощью сплайн-функций. Сплайн (spiain) в переводе с английского языка означает «гибкая линейка». Сплайн-функцию *«ожно наглядно трактовать как линию, которую образует гибкая линейка, будучи закрепленной в ряде точек - узлах иитерполяции. Математически сплайн - специальный многочлен, принимающий в узлах значения у (х) =у, = у (xi) и обеспечивающий непрерывность в них производных. Обычно достаточно обеспечить непрерывность первой и второй производных, дли чего

достаточно использовать сплайн-многочлены третьего порядка (кубические сплайны) [18, 24, 42].

Для каждого отрезка {jr,-, jc+ij изменения х кубическая сплайн-функция записывается в виде

и (X) --

6fti

[т,- х)-)-т,+ , (x-x,)J-f

Г/ m,ft? \

{XiJl-X) -\г

) (х-х.)

где /г, = лг;+,-х„ \,(,х)у(х), т, = Г («/) и /=1, 2,..., п (п - число узлов). При известных х„ у, и mi эта формула задает сплайн-аппроксимацию.

Если потребовать выполнения условия f, {х) •= то приведенное выше выражение для кубических полиномов-сплайнов приводит к системе линейных уравнений, из которых находятся т,:

hitni+2 (if-f-fti+i) /п,-4.im,-+2 =

\ h,+ { ft, J

Однако эта система не полностью определяет пи. Чтобы полностью определить все т., нужно задать дополнительные граничные условия. Если они заданы в виде mi=0 и гп„~0, получаем нормальные сплайн-функции; при m„=mi и mn+i~nti имеем периодические сплайн-функции и т. д.

Приведенная ниже программа для n = N значений и у, вычисляет по данному .алгоритму коэффициенты mi==i" (xs) нормального кубического сплайна и при klxuxn] обеспечивает интерполяцию с помощью описанной выше сплайн-функции. Если x<xt, то выполняется линейная экстраполяция яо формуле

f (х)=У,-{(х,.-х,) т2/6+{у2-у,)/(х2~

-x,)) {х~а,),

а при х> х„ - экстраполяция по формуле

/ {x)yn+i(x„~x..,-i) га„ ,/6 +

+ {Уп-уг.-.,)/{х„-х„.,)) {х-х„).

При экстраполяции предполагается, что наклон линейного участка /(х) равен первой производной сплайн-функции в точке {х\, у,) при экстраполяции для x<Xi и в точке {х,„ у„) при х> х„.

Программа П5.И.

£118 PRINTСПЯййН-АППРОКСИМйЦИЯ,ИНТЕРПОЛЯЦИЯ И ЭКСТРАПОЛЯЦИЯ-•

еге inputзадайте число узлов n=w 838 в I к кн)V ю, L < N) ? и >; W), R < н), S н >

046 FCR 1=1 то NiPRINT!2=0!fl.nH 1=I 858 1НРиТВВЕДИТЕ Kj.VX<I)/Sa):NEKT I

ese LETB=Mi;£>--:;;.;i>3LETE=(:v!:2)-va)).-Ti

878 FOR K=2 TO H-lsLETH=B:LETB=X<K+n-X.;K)

em LETF=E5LETE=<VCK+l)-V>;K))B:LETLtK>=B<B+H5

ЗЭ8 L,ETRtK)=l-L<;K>:LETS(K)=6ж<E-F)/tH-Б): NEXT К

168 FOR K=2 TO H-i5LETP=l<;R(K)!tiL(K-l)+2)!LETL<;K)=-L(K>*P

US lets;k)=(S<k)-R(K)!«S(K-i))*p: next к

128 leth<N)=0sletL;h-i)=S;h-ij.-LETMCH-i)=L;h-i>

138 FOR K=N-2 TO 1 STEP-1

140 LETL<K>=L<K)!t4.<K+l>+S<K)sLETM<K)=L<K)!NEXT К

158 FOR К=1 TO N:PRINT!2.0!MK!F1.9! = M(;K>:NEKT R

166 IHPUT-ВВЕДИТЕ АРГУМЕНТ X=X:l.ETI=8! IFK>X<H> THEH 236

170 IF X<=XC1> THEN £50

180 LETI=I+1!IF Х>ка> THEN 188

130 LETJ=I-1!LETD=X СI>-X С J):LETH=X-X С J):LETK-X iI>-X

200 LETP=Ii*.D.-6: LETS=;M(: J)*R-3+Mf I жН-З.-б-Б

210 LETV=V+«; CVC J)-M< J)-*P>*R+< П-Нс I .!*Р>*И>-Б

220 PRINTЗНАЧЕНИЕ VOO=VjGOTO 168

230 letii=X(:n)-?u:h-i)sLETV=rii*M(:N-i>у&-и\><.н>-\":Н1?у.-ь

240 LETV=Vi«-;X-X<H)>+VO-O!e0T0 220

250 LETD=X с 2 ) -X 1): LETV=-Ii»:M С 2 ) "6+ С V С 2 > -V (1 л ) В

260 LETV=V*«:X-K<i:!>+Va>:GOTO 22e:END

Контрольный пример Пусть надо про- инвертирование матрицы 3X3 (если индексация

вести сплайн-аппроксимацию для интеграла Фре- массивов идет с 1) или 4X4 (если индексация

неля С (л;), заданного таблицей из Л = 7 пар идет с 0). xt и" уг.

х, = 0 0,3 0,5 0,8 ! 1.3 1,5

у, = 0 0.299401 0,4923442 0,7228442 0.7798934 0,6385505 0,4452612

Введя Л, Xi и у„ получим: mi=0, ms=-= - 0,08877862525, тз = - 0,5547268737. т, = = -2,019378004, «15=-3,563639669, тв = = - 1,902729099, im = 0. Далее, задавая х, можно провести интерполяцию и экстраполяцию:!/ (0,2)=0,2000938813 (0,1999211),у (0,48) = = 0,4737408625 (0,473751). у (1.6) =0,3422741197 и т. д. (в скобках после результата даны точные значения).

В приведенных ниже программах индексация массивов идет с О, операторы DIM задают точные размеры матриц (с учетом яринято-й индексации), в операторах МАТ матрицы задаются лишь именами массивов. Даны ссылки иа программы, обеспечивающие те же вычисления без использования матричных операторов и сравнение программ.

Программа П5.12/1.

10 РЕМТРнНСПОНИРОВйНИЕ МАТРИЦЫ А

20 INPUTЗАДАЙТЕ H==Ns .Ч=Н~1!ГаИ HHHJi-BOuН:>

38 DATA 1» 2» 3,4» 5» 6*7» 8/9

40 МнТ REhD asмат B=1RH<A)

•58 PRIHTIPI.SМАССИВ ТРйНСПОИИРОВАИИОй МАТРИЦЫ ПОСТРОЧНО

68 МНТ PRINT B=EHli

Хотя сплайн-интерполяция осуществима при переменном шаге hi--x,j.,~x„ значительное изменение hi (ft,„пкс/А,м„„> 4) нецелесообразно, поскольку в этом случае преимущества сплакн-интерполяций теряются и погрешность ее заметно возрастает. С применением сплайн-интерполяции возможно численное интегрирование, спектральный анализ и другие математические операции.

§ П5.12. Пакет программ с матричными операторами

Применение матричных операторов (см. § 2.2) сокращает текст программ, содержащих операции с матрицами, и позволяет в 5-20 раз (в зависимости от типа ПЭВМ) уменьшить время вычислений. Сокращение времени вычислений обусловлено тем, что основные матричные операции при этом выполняются на машинно-ориентированном языке ЭВМ.

Перед использованием матричных операторов с помощью операторов DIM задаются размеры массивов (основных и служебных). Далее с по-•мощью операторов DATA. МАТ READ, INPUT или МАТ INPUT вводятся элементы исходных матриц. В некоторых версиях бейсика размеры матриц до 10X10 можно не задавать, а размеры больших матриц задавать с запасом. При этом точные размеры матриц указываются в операторе, например, МАТ 1NV А (3,3) означает

Пример. Эта программа обеспечивает транспонирование матрицы, взятой из примера к программе 3.48. После задания размерности N квадратной матрицы-Л, уменьшения /V на 1 и задания размерности массивов А {N. Л) и В (Л, Л) в строке 20 с помощью оператора DATA (строка 30) построчно вводятся элементы исходной матрицы А. В строке 40 оператор МАТ READ А обеспечивает формирование матрицы А, а оператор М.4Т В =TRN (.4) -транспонирование матрицы а. После комментария в строке 50 массив транспонировавной матрицы В выводится построчно на индикацию с помощью оператора МАТ PRINT В (строка 60). Отметим, что в состав последнего в данной версии бейсика нельзя включать знаки изменения формата. Поэтому эти знаки (!Р1.5Г; даны после оператора PRINT в предшествующей строке 50. Результат получается в виде чисел 1, 4, 7, 2, 5, 8, 3, 6 и 9 (соответствует матрице А для контрольного примера к программе 3.48).

Программа П5.12/2.

Эта программа обеспечивает инвертирование матрицы А. Она примерно в 5 раз короче про--граммы 3.56 и обеспечивает уменьшение времени обращения примерно в 20 раз. Значения элементов исходной матрицы и результаты вычислений соответствуют результатам, приведенным для контрольного примера к программе 3.55.