Программы для МК-85. Предварительные замечания

0. Предварительные замечания

0.1. О пользователе

0.1.1 Кому адресована библиотека

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

0.1.2 Что нужно иметь и уметь

Для практической работы с библиотекой нужно иметь:
- карманный компьютер;
- Руководство пользователя к нему;
- конкретную программу действий.

Чтобы провести расчет с помощью библиотечной программы, нужно уметь:
- свести задачу к математической постановке;
- общаться с компьютером на удовлетворительном уровне взаимопонимания.

Чтобы составлять свои программы, нужно уметь:
- читать бейсик-программу;
- изменять программу;
- находить ошибки.

0.1.3 Где искать совместимые программы

Образцы хорошего стиля программирования в приложении к простым задачам можно найти в "Библиотеках программ" карманных компьютеров фирм Casio и Sharp.

Кроме того, можно сравнительно легко адаптировать вычислительные программы, написанные на стандартном бейсике, если только они не содержат матриц. Самый доступный сборник таких программ - Справочник по алгоритмам и программам на языке бейсик для персональных ЭВМ (Дьяконов В.П. - М., Наука, 1987).

Те, кто имеет опыт программирования калькуляторов, могут смело переводить алгоритмы с мнемокода на бейсик. Существует множество доступных источников, например Микрокалькуляторы в физике (Шелест А.Е. - М., Наука, 1988).

0.2. О библиотеке

0.2.1 Классификация программ

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

Библиотека разбита на разделы по типу обрабатываемых математических объектов. Тип объекта определяет структуру данных в программе. Грамотное построение структуры данных - наиболее сложный аспект программирования "Электроники МК 85".

Каждому типу математических объектов соответствует свой набор операций обработки. В выделяемых по этому признаку подразделах материал связан общей постановкой задачи и формой результатов.

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

0.2.2 Стиль программирования

Стиль программирования определен следующими требованиями:
- вычислительные ресурсы "МК 85" ограничены и должны расходоваться экономно;
- текст вводится только с клавиатуры, поэтому при оптимизации программ в спорных случаях следует жертвовать быстродействием в пользу компактности;
- комментарии при вводе и выводе данных должны исключать двусмысленность их толкования;
- программы должны поддаваться модификации, надстраиваться, сопрягаться друг с другом.

0.2.3 Стиль описания программ

При записи текстов программ используются некоторые вспомогательные обозначения. "Y=f{X,U,V}" обозначает определяемую пользователем процедуру-функцию с входными параметрами в ячейках с адресами X,U,V и результатом в ячейке Y. За текстом программы обязательно следует конкретный пример записи.

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

Символы сообщений в операторах INPUT и PRINT набраны обычным шрифтом и в точности соответствуют символам, использованным при описании программы.

Справа перед текстом программы указывается объем памяти, необходимый для ее хранения, а также старший адрес, занимаемый переменными. Эти данные понадобятся, когда Вы столкнетесь с проблемой распределения памяти.

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

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

Для конкурирующих вычислительных методов указаны предпочтительные области применения.

Для программ с динамическими массивами указаны их допустимые размеры (при объеме памяти 1221 шаг).

В конце описательной части приводится пример использования программы. Следование чисел в примере строго соответствует запрограммированной процедуре обмена данными. Приведенные ответы могут в младших разрядах отличаться от значений, получаемых на "МК 85".

0.3. О творческом использовании библиотеки

0.3.1 Внесение изменений в программу

Бейсик "МК 85" не позволяет создать библиотеку объектных модулей, ориентированную на технологию модульного программирования. Правильнее всего рассматривать программы библиотеки как заготовки, отталкиваясь от которых и комбинируя которыми можно без особого труда создать необходимую прикладную программу.

Существуют приемы, позволяющие сопрягать готовые блоки доступными средствами.

При совместном использовании нескольких крупных блоков удобно расположить их в разных программных областях и осуществлять передачу управления операторами GOTO #P, GOSUB #P ... RETURN, IF ... THEN #P, где P - номер программной области.

При внесении в готовый текст дополнительных операторов не следует забывать, что между двумя заполненными строками есть еще девять свободных. Любая многооператорная строка, не содержащая "IF ... ; ...", может быть разбита на несколько последовательных строк без риска нарушить правильность передачи управления, а конструкция вида
20 IF A>B; ОПЕРАТОРЫ
30 ...
заменяется на конструкцию
20 IF A≤B THEN 30
21 - 29 ОПЕРАТОРЫ
30 ...

В "МК 85" имя переменной однозначно связано с ее адресом в памяти, что является досадным отступлением от стандарта бейсика. Осваивая такую версию языка, удобно пользоваться понятиями адресного оператора A(n)= и адресной функции A(n). Первое означает "записать число в ячейку памяти с номером n", второе - "прочесть число из ячейки n".

Если ограничиться обозначениями A(n), программа будет вполне работоспособна, но текст ее можно сделать нагляднее и короче. Бейсик "МК 85" использует систему обозначений, согласно которой

...A(25)=B(24)= ... =Y(1)=Z(0)=Z;
A(24)=B(23)= ... =Y(0)=Y=Z(-1); ...

Пользуясь этим тождеством, из выражения для индекса исключают постоянное слагаемое, например, I(A-1) заменяют на H(A). Если адрес не имеет переменной составляющей, скобки не нужны (например, I(1) заменяют на H). Следовательно, все скаляры удобно хранить в младших адресах, а вектора переменной размерности - в старших.

В библиотечных программах переменные распределяются по адресам подряд, начиная с A. Так легче следить за расходом памяти.

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

строка | A | B | C | D | E | F
10 | a | b | c | x | yn | n
20 |  |  |  |  | yj | j
30 |  |  |  |  | y0 | -1

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

0.3.2 Средства обнаружения ошибок

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

Набирая программу, хорошо иметь перед глазами ее текст на бумаге.

Если программа состоит из нескольких крупных блоков, надежнее будет сначала проверить блоки по отдельности, а уж затем - программу целиком.

Многие описки обнаруживаются автоматически встроенной диагностикой бейсика. Много хлопот вызывают другие ошибки - проявляющиеся лишь в неверности результатов расчета. Для обнаружения таких ошибок применяют тест - набор контрольных примеров с известными ответами. Тест должен быть таким, чтобы срабатывала каждая ветвь алгоритма. Исходные данные должны быть нетривиальными. (Не следует, в частности, задавать x=1, если в программе x возводится в степень.) Кстати, примеры, приведенные в библиотеке, не являются тестами.

Поиск ошибки при ее обнаружении сводится к отслеживанию передачи управления и просмотру промежуточных результатов.

Просмотр ячеек с переменными производится в калькуляторном режиме после приостановки счета программы. Например, для просмотра содержимого ячейки A в момент, когда программа выполняется, следует нажать [STOP]. На дисплее появится номер текущей строки программы. Теперь наберем a [EXE]. Дисплей высветит значение переменной A. Повторное нажатие [EXE] запустит программу с точки останова.

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

Оператор RUN n. Позволяет запускать программу не сначала, а с n-ой строки.

Отладочный вывод промежуточного значения. Пример: 33 CSR 0:PRINT A;. CSR отменяет режим бегущей строки, ";" отменяет прерывание после вывода. Используется обычно для контроля сходимости итеративного процесса.

Клавиша останова [STOP]. Ее пускают в ход, когда есть подозрение, что алгоритм "зациклился".

Оператор STOP. Вызывает прерывание счета программы. Счет продолжится, если нажать [EXE].

Режим пооператорного исполнения программы TRACE. Наиболее мощное отладочное средство, которое применяется преимущественно для контроля передачи управления.

Серьезное разочарование вызывают огрехи в микросхемах и встроенном математическом обеспечении. Проведите такой опыт:

[MODE] [5] COS 0

- Образцы 1990г. выпуска возвращают -1. Похожие проблемы возникают и при попытке вычислить ATN 0 в любом из режимов. Более подробное исследование ошибок опубликовано в журнале "Наука и жизнь" N 6 1991г., стр. 133. Разумеется, программы библиотеки не учитывают столь экзотические свойства встроенной математики. Остается только порекомендовать Пользователю проявлять бдительность при вызове тригонометрических функций.

Еще один полезный опыт для всех:

10 A=1.00000000001:B=1:IF A¦B;PRINT A;" ¦";B

- Проверки отношения распространяются и на скрытые разряды.

0.3.3 Средства диалога

Удобство программы при ее неоднократном использовании определяется в первую очередь тем, как организован ввод-вывод данных. Программы библиотеки не изобилуют разнообразием форм диалога. Вы можете "перекроить" ввод-вывод по своему вкусу.

Вот некоторые варианты:

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

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

Вычисляемый комментарий. Пример: PRINT A;:INPUT B. Применяется, когда значение B зависит от значения A.

Вывод по укороченному формату. Может устанавливаться оператором SET n, который управляет форматом вывода (n - число знаков мантиссы). Альтернативный путь - использование в операторе PRINT функции RND (A,n), которая округлит число A до n-го разряда. Укороченный формат позволяет ускорить вывод и улучшить наглядность результатов при низкоточных вычислениях.

Вывод без прерывания. Избавляет от необходимости каждый раз нажимать [EXE] для продолжения счета. Пример: PRINT A;.

Отмена "бегущей строки". Осуществляется оператором управления курсором CSR и позволяет резко снизить затраты процессорного времени на вывод.

Прерывание перед выводом. Исключает возможность потери информации при выводе строки, содержащей более 12 символов. Пример: 10 PRINT "E",-EXP 3000.

Ввод выражения в ответ на запрос оператора INPUT. Используется, когда входную переменную программы нужно предварительно вычислить либо она зависит от переменных, хранящихся в памяти.

Например, вместо длинного числа 1,57079633 можно ввести выражение p/2.


Программы для МК-85. Числа

1. Числа

1.0. Введение

Число - простой математический объект, является базовым при построении более сложных объектов. Как правило, для хранения числа в программе нужно выделить одну переменную. Многие языки программирования различают переменные по типам и по точности.

1.0.1 Действительные числа

В бейсике "МК 85" всего один тип числовых переменных - действительные с плавающей точкой (есть также строки символов). На одну переменную расходуется 8 байт памяти, которые позволяют закодировать 12 десятичных разрядов мантиссы и порядок в диапазоне +4096. 2 младших разряда мантиссы на дисплей не выводятся и называются скрытыми. Благодаря им выводимое число всегда корректно округлено, и ситуация, когда 1 отображается как 0.9999999999, практически исключена.

Арифметические операции над действительными числами, поддерживаемые бейсиком:
A+B - сложение; A-B - вычитание;
A*B - умножение; A/B - деление.

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

1.0.2 Логические переменные

Бейсик "Электроники МК 85" дает возможность проводить следующие операции сравнения чисел:
A¦B - A не равно B; A=B - A равно B;
A<B - A меньше B; A≤B - A не больше B;
A>B - A больше B; A≥B - A не меньше B.

Эти операции разрешены только в операторе IF, который является арифметическим, а не логическим, как в стандартном бейсике.

Для организации логической переменной приходится использовать либо числовую (A), либо символьную ($, A$) переменную бейсика. Переменную, играющую роль логической, называют флагом. Операцию A=1 называют выставить флаг, операцию A=0 - сбросить флаг, IF A=1; - проверка.

1.0.3 Бинарные числа

Принимают всего два значения: 0 и 1. В одной переменной программы может храниться в упакованном виде до 50 таких чисел. Однако в "Библиотеке" нет ни одной программы, работающей с бинарными матрицами, так как их упаковка и распаковка - слишком трудоемкие операции.

1.0.4 Комплексные числа

Для хранения комплексного числа используется пара действительных переменных. Справедливы две формы определения комплексного числа: алгебраическая Z=A+jB и тригонометрическая Z=rejq, где j=Ц-1, A=Re Z - действительная часть, B=Im Z - мнимая часть, r- модуль, q - фаза.

Процедуры перевода числа из тригонометрической формы в алгебраическую и обратно:

RAD :A=C*COS D:B=C*SIN D
RAD :C=SQR (A*A+B*B):D=ASN (B/C)

Здесь C - r¦0, D - q.

Арифметические операции над комплексными числами в алгебраической форме:

A=C+E:B=D+F - сложение;
A=C*E-D*F:B=C*F+E*D -
  умножение;
B=E*E+F*F:A=(C*E+D*F)/B:B=(D*E-C*F)/B
  - деление;
где (A,B) - результат,
  (C,D),(E,F)
  - операнды.

1.1. Числа с произвольным основанием

Число состоит из нескольких разрядов, между которыми определено старшинство. Количество возможных комбинаций разряда называют основанием числа. Одно и то же число может иметь разные формы записи в зависимости от используемого основания. Во многих диалектах бейсика предусмотрена возможность вводавывода чисел в форме двоичного, восьмеричного и шестнадцатиричного. В бейсике "МК 85" допустим только десятеричный формат.

Преобразование десятеричного числа в число с другим основанием проводится по формуле:

d = е(i)eimi, i=l...h

где d - десятеричное число;
h - номер старшего ненулевого разряда;
l - номер младшего значимого разряда;
m - основание;
ei - значение i-го разряда числа с основанием m (0≤ei<m).

Программа 1. Перевод произвольного числа в десятеричное.

10 D=0:INPUT "m",A,"h",B,"l",C:FOR B=B TO C STEP -1:PRINT B;
20 INPUT E:D=D+E*A-B:NEXT B:PRINT D:GOTO 10

Размер: 64, E

Пример: m=100; h=2; l=-1; e2=45; e1=28; e0=81; e-1=32.

Ответ: d=452881,32.

Программа 2. Перевод десятеричного числа в число с требуемым основанием.

10 INPUT "d",D,"m",A:B=INT (LOG D/LOG A+1):E=B-10/LOG A
20 C=INT (D/A-B):IF C>0;PRINT B;C:D=D-C*A-B
30 B=B-1:IF D"`0;IF B>E THEN 20

Размер: 92, E

Пример: d=452881,32; m=100.

Ответ: e2=45; e1=28; e0=81; e-1=32.

Программа 3. Перевод десятеричного числа в число с основанием от 2 до 16. Значения от 10 до 15 выводятся на дисплей буквами от A до F.

10 INPUT "d",D,"m",A:B=INT (LOG D/LOG A):E=B-10/LOG A
20 C=INT (D/A-B):D=D-C*A-B:IF C>9;C=C+7
30 PRINT CHR (48+C);:B=B-1:IF B≥0 THEN 20
40 IF D¦0;IF B>E THEN 60
50 GOTO 10
60 IF B=-1;PRINT ".";
70 GOTO 20

Размер: 137, E

Пример: d=63; m=16.

Ответ: Hex d=3F.

1.2. Делимость натуральных чисел

В результате деления A на B получаются частное D и остаток R такие, что R<B и B*D+R=A. Если R=0, то говорят, что A кратно B. Натуральное число, которое делится только на само себя и на 1, называют простым.

Программа 4. Наибольший общий делитель двух чисел.

10 INPUT A,B
20 C=A:D=INT (B/A):A=B-D*A:IF A¦0;B=C:GOTO 20
30 PRINT C:GOTO 10

Размер: 52, D

Пример: 625; 1225.

Ответ: НОД(625;1225)=25.

Программа 5. Разложение числа на простые сомножители.

10 B=2:INPUT A
20 IF A<B*B;PRINT A:GOTO 10
30 C=A/B:IF INT C=C;A=C:PRINT B:GOTO 30
40 B=B+1:IF B>3;B=B+1
50 GOTO 20

Размер: 75, C

Пример: 1234567890.

Ответ: 2; 3; 3; 5; 3607; 3803.

Программа 6. Наименьшее общее кратное двух чисел.

10 INPUT A,B:E=1:D=2
20 IF A=1;IF B=1;PRINT E:GOTO 10
30 F=0:C=A/D:IF INT C=C;A=C:F=1
40 C=B/D:IF INT C=C;B=C:F=1
50 IF F=1;E=E*D:GOTO 30
60 D=D+1:IF D¦3;D=D+1
70 GOTO 20

Размер: 124, F

Пример: 45; 125.

Ответ: НОК(45;125)=1125.

Программа 7. Наибольший общий делитель нескольких чисел.

10 INPUT "n",A:C=1:D=1:FOR B=1 TO A:PRINT B;:INPUT D(B):NEXT B
20 FOR B=1 TO A:IF D(B)≤D;PRINT C:GOTO 10
30 NEXT B:D=D+1:IF D>3;D=D+1
40 FOR B=1 TO A:IF FRAC (D(B)/D)¦0 THEN 20
50 NEXT B:C=C*D:FOR B=1 TO A:D(B)=D(B)/D:NEXT B:GOTO 40

Размер: 145, D(n)

Пример: n=3; m1=45; m2=125; m3=225.

Ответ: НОД(45;125;225)=5.

Программа 8. Наименьшее общее кратное нескольких чисел.

10 INPUT "n",A:C=1:D=1
20 FOR B=1 TO A:PRINT B;:INPUT F(B):NEXT B
30 FOR B=1 TO A:IF F(B)>1 THEN 50
40 NEXT B:PRINT C:GOTO 10
50 D=D+1:IF D>3;D=D+1
60 F=0:FOR B=1 TO A:E=F(B)/D:IF INT E=E;F(B)=E:F=1
70 NEXT B:IF F=0 THEN 30
80 C=C*D:GOTO 60

Размер: 154, F(n)

Пример: n=3; m1=45; m2=125; m3=225.

Ответ: НОК(45;125;225)=1125.

1.3. Комбинаторика

Программа 9. Число перестановок. Расчетная формула:

n! = Х i, i=1...n

10 INPUT A:B=1:FOR C=1 TO A:B=B*C:NEXT C:PRINT B:GOTO 10

Размер: 32, C

Пример: n=5.

Ответ: 5!=120.

Примечание: Программа дает правильный ответ и при n=0.

Программа 10. Число размещений.

Amn=n! / (n-m)!

10 INPUT "n",A,"m",B:C=1:D=A-B
20 D=D+1:IF D≤A;C=C*D:GOTO 20
30 PRINT C:GOTO 10

Размер: 57, D

Пример: n=8; m=5.

Ответ: A58=6720.

Примечание: Алгоритм не повторяет приведенную выше формулу. Прямое ее копирование привело бы к увеличению затрат памяти, дополнительной погрешности вычислений, большей трудоемкости и более узкому диапазону допустимых n.

Программа 11. Число сочетаний (биномиальные коэффициенты).

Cmn=n! / m! / (n-m)!

10 INPUT "n",A,"m",B:C=1:D=B
20 D=D+1:IF D≤A;C=C*D/(A-D+1):GOTO 20
30 PRINT C:GOTO 10

Размер: 63, D

Пример: n=8; m=5.

Ответ: C58=56.

Примечание: Использование вместо оператора D=B последовательности D=A-B:IF D<B;D=B улучшает быстродействие алгоритма при m < n/2.

1.4. Числовые последовательности

В основе любого итерационного (многошагового) алгоритма лежит представление о некоторой последовательности

{xi}, i=0,1,2,3,...

где xi - член последовательности,
i - его порядковый номер.

Программа 12. Последовательность простых чисел.

Алгоритм относится к классу переборных: искомая последовательность получается путем поэлементного отсеивания членов базовой последовательности. Попытки как-то сузить область перебора получили название метода решета. В программе использовано самое примитивное решето, исключающее из рассмотрения все четные числа после двойки.

10 PRINT 2:A=1
20 A=A+2:B=3
30 IF A<B*B;PRINT A:GOTO 20
40 IF FRAC (A/B)=0 THEN 20
50 B=B+2:GOTO 30

Размер: 64, B

Ответ: 2; 3; 5; 7; ...

Программа 13. Последовательность чисел Бернулли. Расчетная формула:

B0=1; Bk=1/(k+1) е(i)Bi (k+1)! / (k+1-i) / i!, i=1 ... k-1

Этот алгоритм относится к нерекуррентным, так как для вычисления нового члена последовательности используются все предыдущие.

10 E=1:A=1
20 A=A+1:D(A)=0:C=1/A:B=0
30 B=B+1:D(A)=D(A)-C*D(B):IF B<A-1;C=C*(A-B+1)/B:GOTO 30
40 PRINT A-1;D(A):GOTO 20

Размер: 101, E(k)

Ответ: B1=-1/2; B2=1/6; B3=0; ...

Программа 14. Последовательность s-чисел Фибоначчи.

Mi=1, i=0...s;
Mi+1=Mi+Mi-s, i>s.

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

Как правило, мы имеем дело с рекуррентными алгоритмами.

10 INPUT A:FOR B=0 TO A:D(B)=1:NEXT B
20 E(A)=D(A)+D:PRINT D;:FOR B=0 TO A:D(B)=E(B):NEXT B:GOTO 20

Размер: 64, E(s)

Пример: s=1.

Ответ: M0=1; M1=1; M2=2; M3=3; M4=5; ...

1.5. Сортировка чисел

Сортировка сводится к упорядочиванию некоторого перечислимого числового множества по возрастанию (убыванию) его членов:

{xi}, xi+1 ≥ xi, i=1 ... n-1

Программа 15. Метод "пузырька". Программа переставляет индексы r исходного массива чисел er так, чтобы числа er(i) расположились в порядке возрастания по i.

10 INPUT "n",A:FOR B=1 TO A:C=B*2:PRINT
  B;
20 INPUT E(C):D(C)=C:NEXT B:FOR B=2 TO A:C=B*2
30 C=C-2:IF E(D(C))≤E(F(C)) THEN 50
40 D=D(C):D(C)=F(C):F(C)=D:IF C>2 THEN 30
50 NEXT B:FOR B=1 TO A:PRINT B;D(B+B)/2:NEXT B

Размер: 146, F(2*n)

Пример: n=4; e1=103; e2=17; e3=67; e4=25.

Ответ: r1=2; r2=4; r3=3; r4=1.

Примечание: Для упорядочивания по убыванию измените знак отношения в строке 30 на ≥.


Программы для МК-85. Матрицы

2. Матрицы

2.0. Введение

С точки зрения программирования матрица - это двумерный массив чисел. Так как бейсик "МК 85" не поддерживает массивы, организация матрицы сопряжена с некоторыми трудностями: приходится следить за тем, чтобы в участок адресов, отведенный под матрицу, не попадали другие переменные.

Матрица обрабатывается поэлементно. Обход элементов чаще всего организован как двойной вложенный цикл:

FOR I=0 TO A-1
FOR J=0 TO B-1
PRINT "M(";I;",";J;")=";M(I*B+J)
NEXT J
NEXT I

Здесь внутренний FOR изменяет номер столбца J, а внешний - номер строки I.

Транспонирование матрицы сводится к изменению порядка обхода:

FOR J=0 TO B-1
FOR I=0 TO A-1
PRINT "M(";J;",";I;")=";M(I*B+J)
NEXT I
NEXT J

Симметричные матрицы могут храниться в уплотненной форме - в "треугольном" массиве:

FOR I=0 TO A-1
FOR J=I TO A-1
PRINT "M(";I;",";J;")=";M(I*A+J)
NEXT J
NEXT I

Записаны только элементы, находящиеся на главной диагонали и справа от нее. При обращении к матрице должно выполняться отношение I≤J, что можно обеспечить введением новых индексов:

FOR I=0 TO A-1
FOR J=0 TO A-1
K=I:L=J:IF I>J;L=I:K=J
PRINT "M(";I;",";J;")=";M(K*A+L)
NEXT J
NEXT I

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

FOR I=0 TO A
K=INT N(I+I):L=100*FRAC N(I+I)
PRINT "M(";K;",";L;")=";M(I+I)
NEXT I

Здесь A - число ненулевых элементов матрицы, M(I+I) - значение I-го ненулевого элемента, N(I+I) - его индексы, упакованные в одну ячейку. Номер строки составляет целую часть, а номер столбца - дробную: N(I+I)=K+L/100. Формула записана для случая, когда число столбцов меньше ста. Уплотнение оправдано лишь тогда, когда число нулевых элементов больше, чем ненулевых.

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

Представленные в разделе программы записаны для матриц произвольного размера. Работая с матрицей размера 2 или 3, нет смысла использовать такие длинные программы. Применяйте явное выражение детерминанта. Примеры:

1. Вычисление детерминанта матрицы 3-го порядка.

10 FOR A=1 TO 3:FOR B=1 TO 3:PRINT A;B;:INPUT A(B+3*A):NEXT B
20 NEXT A:D=E*I*M-E*J*L+F*J*K-F*H*M+G*H*L-G*I*K:PRINT D

Адресация элементов матрицы: при вводе - индексная, при вычислении детерминанта - буквенная. Размещение матрицы в памяти:

ж E F G ц
з H I J ч
и K L M ш

2. Решение системы из двух линейных уравнений методом Крамера.

10 INPUT "A11",A,"A12",B,"A21",C,"A22",D,"B1",E,"B2",F
20 G=A*D-B*C:PRINT "X1";(E*D-B*F)/G,"X2";(A*F-E*C)/G

2.1. Матричное преобразование

r'=Ar

Программа 16. Матрица вращения в трехмерном пространстве - типичный пример матричного преобразования. Здесь r - исходный радиус-вектор точки; r' - радиус-вектор точки после поворота вокруг начала координат на углы a, b, g; A - матрица вращения:

cosa cosb | -cosa sinb cosg + sina sing | cosa sinb sing + sina cosg
sinb | cosb cosg | -cosb sing
-sina cosb | sina sinb cosg + cosa sing | -sina sinb sing + cosa cosg

10 MODE 4:INPUT "a",A,"b",B,"g",C:F=COS A:I=SIN B:J=COS C
20 K=SIN C:L=SIN A:G=L*K-F*I*J:H=F*I*K+L*J:M=L*I*J+F*K
30 N=F*J-L*I*K:D=COS B:F=F*D:J=J*D:K=-K*D:L=-L*D
40 INPUT "x",A,"y",B,"z",C:FOR D=0 TO 2:O=0:FOR E=0 TO 2
50 O=O+A(E)*F(E+D*3):NEXT E:PRINT CHR (88+D);O:NEXT D:GOTO 40

Размер: 212, O

Пример: a=30; b=45; g=60; x=1; y=0; z=0.

Ответ: x'=0,612372439; y'=0,70710678; z'=-0,353553389.

2.2. Произведение матриц

C=AB

Программа 17. Расчетная формула:

cij=е(k)aikbkj , k=1...r

10 INPUT "m",A,"r",B,"n",C:DEFM A*B+B*C-18:FOR G=0 TO 1
20 FOR D=0 TO A(G)-1:FOR E=1 TO B(G):PRINT D+1;E;
30 INPUT G(E+D*B(G)+G*A*B):NEXT E:NEXT D:NEXT G:FOR D=1 TO A
40 FOR E=1 TO C:F=0:FOR G=0 TO B-1:F=F+H(G+D*B-B)*G(E+G*C+A*B)
50 NEXT G:PRINT D;E;F:NEXT E:NEXT D

Размер: 177, H(r*(m+n)). Ограничение размерности (n+m)r≤128, в частности, для квадратных матриц: n≤8.

Пример: m=2; r=3; n=1; A=((2, 7, 4), (8, 1, 3)); B=(6, 2, 4) T

Порядок ввода данных: 2; 3; 1; 2; 7; 4; 8; 1; 3; 6; 2; 4.

Ответ: C=(42, 64) T

2.3. Детерминант. Обратная матрица

Детерминант определяется формулой:

det A=е(k) (-1)k Х(l) ai(l,k),j(l,k); l=1...n

где k - номер перестановки,
n - размерность матрицы A,
aij - элементы матрицы.

Обратная матрица A-1 - это матрица, удовлетворяющая соотношению A-1A=E, где E - единичная матрица.

Программа 18. Модифицированный метод Гаусса - Жордана. Выбор ненулевого диагонального элемента обеспечивается перестановкой строк. Результат возвращается в ячейки исходного массива.

10 INPUT E:DEFM E*E+E-20:FOR A=1 TO E:FOR B=1 TO E
20 PRINT A;B;:INPUT F(A*E+B):NEXT B:NEXT A
30 D=1:FOR A=1 TO E:IF F(A*E+A)¦0;F(A)=A:GOTO 80
40 FOR B=1 TO E:IF F(B*E+A)¦0 THEN 60
50 NEXT B:PRINT 0:END
60 D=-D:FOR C=1 TO E:F=F(A*E+C)
70 F(A*E+C)=F(B*E+C):F(B*E+C)=F:NEXT C:F(A)=B
80 F(A*E+A)=1/F(A*E+A):FOR C=1 TO E:IF C=A THEN 100
90 F(A*E+C)=F(A*E+A)*F(A*E+C)
100 NEXT C:FOR B=1 TO E:IF A=B THEN 130
110 F=F(B*E+A):F(B*E+A)=0:FOR C=1 TO E
120 F(B*E+C)=F(B*E+C)-F*F(A*E+C):NEXT C
130 NEXT B:D=D/F(A*E+A):NEXT A
140 PRINT D:FOR A=1 TO E:FOR B=1 TO E:F=F(A*E+F(B))
150 IF F(B)¦B;F(A*E+F(B))=F(A*E+B):F(A*E+B)=F
160 PRINT A;B;F:NEXT B:NEXT A

Размер: 453, F(n*(n+1)). Допустимая размерность n≤10.

Пример: n=3; A=

ж1 0 1ц
з1 1 0ч
и0 1 1ш

Порядок ввода данных: 3; 1; 0; 1; 1; 1; 0; 0; 1; 1.

Ответ: det A=2; A-1=

ж0.5 0.5 -0.5ц
з-0.5 0.5 0.5ч
и0.5 -0.5 0.5ш

Проверка с использованием свойства обратимости: (A-1)-1 = A.

Повторный запуск программы командой RUN 30 дает исходную матрицу с удвоенной вычислительной погрешностью. Для плохо обусловленных A погрешность может быть значительной.

2.4. Характеристический многочлен

P(l) = det(A-lE)

Программа 19. Метод Данилевского. Для невырожденной матрицы программа возвращает коэффициенты характеристического многочлена вида

P(l) = ln + е(i)cili-1 , i=1...n.

Вырожденной матрице соответствует несколько независимых многочленов. Тогда программа выводит значения k, по которым можно определить степени многочленов. Степень первого многочлена равна k1, второго - (k2-k1), третьего - (k3-k2) и т.д. Первые k1 выводимых на дисплей коэффициентов относятся к первому многочлену, вторые (k2-k1) - ко второму и т.д. В этом случае коэффициенты выводятся с противоположным знаком.

10 INPUT A:DEFM A*A+A-20:FOR B=1 TO A:FOR C=1 TO A
20 PRINT B;C;:INPUT F(C+A*B):NEXT C:NEXT B:D=1
30 FOR B=D+1 TO A:E=0:FOR C=B TO A
40 IF E<ABS E(B+C*A);E=ABS E(B+C*A):F=C
50 NEXT C:IF E>0;FOR C=B-1 TO A:E=F(C+A*B):GOTO 70
60 FOR C=D TO B-1:F(C)=E(B+C*A):NEXT C:D=B:PRINT D-1;:GOTO 30
70 F(C+A*B)=F(C+F*A):F(C+F*A)=E:NEXT C:FOR C=D TO A:E=F(B+C*A)
80 F(B+C*A)=F(F+C*A):F(F+C*A)=E:NEXT C:FOR C=D TO A
90 F(C)=E(B+C*A):NEXT C:E=E(B+A*B):FOR C=B-1 TO A
100 F(C+B*A)=F(C+B*A)/E:FOR F=D TO A
110 IF F¦B;F(C+F*A)=F(C+F*A)-F(F)*F(C+B*A)
120 NEXT F:NEXT C:FOR C=D TO A:E=0:FOR F=B TO A
130 E=E+F(F+C*A)*F(F):NEXT F:IF C≤B;IF C¦D;E=E+E(C)
140 F(B+C*A)=E:NEXT C:NEXT B:FOR B=D TO A:F(B)=F(B*A+A):NEXT B
150 PRINT "":FOR B=1 TO A:PRINT B;F(B):NEXT B

Размер: 521, F(n*(n+1)). Ограничение размерности n≤9. Можно модифицировать программу для расчетов с n=10: убрать INPUT A:DEFM A*A+A-20 из строки 10 и строку вывода 150, а перед RUN дать команды A=10 и DEFM 90. После счета коэффициенты можно вызвать из ячеек G...P.

Пример: n=4; A=

ж1 1 4 3ц
з1 1 2 1ч
з0 0 4 3ч
и0 0 2 1ш

Порядок ввода данных:4; 1; 1; 4; 3; 1; 1; 2; 1; 0; 0; 4; 3; 0; 0; 2; 1.

Ответ: k1=2; c1=0; c2=2; c3=2; c4=5;
что соответствует многочленам: P1(l) = l2-2l; P2(l) = l2-5l-2.

2.5. Собственные значения

Собственные значения матрицы являются корнями характеристического многочлена:

P(li)=0, i=1...n

Для определения собственных значений несимметричных матриц пользуйтесь программой 19, затем 32. Симметричные матрицы приводятся к диагональным, что и реализовано в следующей программе.

Программа 20. Метод вращений с преградами. Матрица приводится к диагональной применением рекуррентной формулы

A0=A; Ak=UkTAk-1Uk;

где Uk - матрица плоского вращения, аннулирующая внедиагональный элемент, по модулю превышающий текущее значение барьера. От цикла к циклу величина барьера уменьшается.

В целях снижения вычислительной трудоемкости программа предусматривает задание погрешности - требуемого числа e верных знаков результата.

10 INPUT "e",B,"n",A:DEFM A*(A+1)/2-12:E=0:FOR C=1 TO A
20 FOR D=1 TO C:PRINT C;D;:INPUT F:IF C¦D;E=E+2*F*F
30 N(C*(C-1)/2+D)=F:NEXT D:NEXT C:E=SQR E/A:B=E*10--B
40 $=CHR 1:FOR C=2 TO A:FOR D=1 TO C-1:F=N(C*(C-1)/2+D)
50 IF ABS F<E THEN 150
60 $="":G=N(C*(C+1)/2):H=N(D*(D+1)/2):I=(H-G)/2
70 J=1:IF I¦0;J=SGN I*F/SQR (F*F+I*I)
80 J=J/SQR (2*(1+SQR (1-J*J))):I=SQR (1-J*J)
90 FOR K=1 TO A:FOR L=0 TO 1
100 M(L)=C(L)*(C(L)-1)/2+K:IF C(L)<K;M(L)=K*(K-1)/2+C(L)
110 NEXT L:L=N(M)*I-N(N)*J:IF C¦K;N(N)=N(N)*I+N(M)*J
120 N(M)=L:NEXT K:L=2*F*I*J
130 N(C*(C-1)/2+D)=(G-H)*J*I+F*(I*I-J*J)
140 N(D*(D+1)/2)=H*I*I+G*J*J+L:N(C*(C+1)/2)=H*J*J+G*I*I-L
150 NEXT D:NEXT C:IF $="" THEN 40
160 IF E>B;E=E/10:GOTO 40
170 FOR C=1 TO A:PRINT N(C*(C+1)/2):NEXT C

Размер: 570, N(n*(n+1)/2). Ограничение размерности n≤13.

Пример: e=6; n=4; A=

ж2======ц
з8 3====ч
з4 6 1==ч
и5 2 9 0ш

Ответ: l1 =-9,97435311; l2 =18,5917142; l3 =-4,26123764; l4 =1,64387682.

2.6. Системы линейных уравнений

Ax=b

К системам линейных уравнений с квадратной невырожденной матрицей применяют как точные методы, основанные на гауссовском приведении матрицы к треугольному виду, так и приближенные, использующие процесс пошаговой компенсации невязок. Последние хороши тем, что во-первых, не столь подвержены проклятью размерности, а во-вторых, позволяют учитывать разреженность матрицы коэффициентов. В этом разделе представлено по одному методу из каждого класса - метод Якоби и метод Зейделя, программы для которых оказались наиболее компактными. Системы с переопределенной матрицей и системы с вырожденной матрицей предварительно приводятся к нормальной системе (Используйте программу 17). Для переопределенной матрицы нормальная система получается по следующим формулам:

A'x=b'; где A'=A T A; b'=A T b.

Решение x соответствует минимуму скалярного произведения

(Ax-b)(Ax-b).

В случае вырожденной матрицы решается возмущенная система

(A T A+e E)x=A T b.

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

Нормализованная система решается с помощью программы 21.

Программа 21. Метод вращений.

10 INPUT "n",A:DEFM A*A+A-19:FOR B=1 TO A:FOR C=0 TO A
20 PRINT B;C;:INPUT G(C*A+B):NEXT C:NEXT B
30 FOR B=1 TO A-1:FOR C=B+1 TO A
40 F=G(B*A+B):G=G(B*A+C):IF F=0;IF G=0 THEN 80
50 E=SQR (F*F+G*G):D=G/E:E=F/E
60 FOR F=0 TO A:G=E*G(F*A+B)+D*G(F*A+C)
70 G(F*A+C)=E*G(F*A+C)-D*G(F*A+B):G(F*A+B)=G:NEXT F
80 NEXT C:NEXT B:FOR B=A TO 1 STEP -1
90 D=0:IF B<A;FOR C=B+1 TO A:D=D+G(C)*G(C*A+B):NEXT C
100 G(B)=(G(B)-D)/G(B*A+B):NEXT B
110 FOR B=1 TO A:PRINT B;G(B):NEXT B

Размер: 348,G(n*n+n). Допустимая размерность n≤10.

Пример: n=3; b=(31, 32, 33)T; A=

ж=0=30===1ц
з45 =0 -13ч
и12 13 ==8ш

Порядок ввода данных: 3; 31; 0; 30; 1; 32; 45; 0; -13; 33; 12; 13; 8.

Ответ: x=(1, 1, 1)T

Метод вращений применен также в программах: 81, 146.

Программа 151. Метод Гаусса с выбором максимального элемента.

10 INPUT A:D=A*(A+2):DEFM D+A-15
20 FOR B=1 TO A:FOR C=1 TO A+1:PRINT B;C;
30 INPUT K(C*A+B):NEXT C:NEXT B
40 FOR B=1 TO A:K(D+B)=B:NEXT B
50 FOR B=1 TO A-1:E=B*A:F=0:FOR H=B+1 TO A
60 I=K(E+K(D+H)):IF ABS I>F;F=ABS I:G=H
70 NEXT H:GOSUB 200
80 J=K(D+G):K(D+G)=K(D+B):K(D+B)=J
90 FOR H=B+1 TO A:G=K(D+H):F=K(E+G)/I:K(E+G)=F
100 FOR K=B+1 TO A+1:C=K*A:K(C+G)=K(C+G)-F*K(C+J)
110 NEXT K:NEXT H:NEXT B
120 FOR B=A TO 1 STEP -1:J=K(D+B):I=K(D-A+J)
130 IF B<A;FOR H=B+1 TO A:I=I-K(H*A+J)*K(H):NEXT H
140 IF I=0 THEN 160
150 F=K(B*A+J):GOSUB 200:I=I/F
160 K(B)=I:NEXT B
170 FOR B=1 TO A:PRINT "x";B,K(B):NEXT B:END
200 IF F¦0;RETURN
210 PRINT "DET=0":END

Размер: 448, L(n(n+3)). Предельная размерность n≤9.

Примечание: Метод Гаусса представлен еще в программе 105.

Программа 22. Метод Зейделя для разреженных матриц. Для обеспечения сходимости следует провести замену переменных так, чтобы для всех строк i=1...n выполнялось:

|aii|>|aij|, j=1...n, j¦i.

Ввод данных. Сначала задается размерность системы n, затем допустимая абсолютная ошибка e. Система считывается построчно. Для i=1...n вводится диагональный элемент матрицы коэффициентов aii, затем свободный член bi, затем ненулевые внедиагональные элементы aij текущей строки. Для ввода aij необходимо сначала указать столбец j, а затем - само значение aij. Для перехода к следующей строке вводите j=0. Начальное приближение x0 задается командой VAC. При необходимости увеличитьточность пользуйтесь командами: B=1E-10, RUN P1, что вызовет выполнение добавочных итераций.

DEFM 3*n+2*r-20
VAC

В программной области P0:

10 INPUT "n",A,"e",B:D=A*3:FOR C=1 TO A:PRINT C
20 INPUT "a",C(C*3),"b",D(C*3)
30 INPUT "j",E:IF E=0;NEXT C:END
40 D=D+2:INPUT "a",E(D):D(D)=C+E/100:GOTO 30

В программной области P1:

10 $=CHR 1:D=A*3:FOR C=3 TO A*3 STEP 3:E=D(C)
20 IF INT F(D)=C/3;D=D+2:E=E-E(D)*E(300*FRAC D(D)):GOTO 20
30 E=E/C(C):IF ABS (E(C)-E)≥B;$=""
40 E(C)=E:NEXT C:IF $="" THEN 10
50 FOR C=1 TO A:PRINT C;E(C*3):NEXT C

Размер: 110+148, F(3*n+2*r). Ограничение размерности при поочередной загрузке областей P0 и P1: 3n=2r≤154, где n - размерность системы, r - количество ненулевых внедиагональных элементов матрицы A.

Пример: e =10-6

м10x1+5x4=15
н10x2+2x3=12
п10x3-3x1= 7==
о10x4+6x2=16

Порядок ввода данных: 4; 1E-6;
10; 15; 4; 5; 0;
10; 12; 3; 2; 0;
10; 7; 1; -3; 0;
10; 16; 2; 6; 0.

Ответ: x=(1, 1, 1, 1) T

Метод Зейделя применен в программах: 80, 101, 102, 145.


Программы для МК-85. Степенные многочлены

3. Степенные многочлены

3.0. Введение

Степенным многочленом (полиномом) называют функцию вида

P(x)=е(i)aixi, i=0...n

где ai - действительный коэффициент,
x - действительная или комплексная переменная,
n - степень многочлена.

В бейсик-программе многочлен может быть представлен в виде вектора коэффициентов {ai}, i=0...n.

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

Раздел содержит широкий набор алгоритмов обработки многочленов.

3.1. Вычисление значений

Быстрый способ вычисления значений многочлена - схема Горнера:

P(x)=(...(anx+an-1)x+...)x+a0.

Программа 23. Случай действительного аргумента.

10 INPUT "n",A:FOR B=0 TO A:PRINT A-B;:INPUT E(B):NEXT B
20 INPUT C:D=0:FOR B=0 TO A:D=D*C+E(B):NEXT B:PRINT D:GOTO 20

Размер: 68, E(n)

Пример: P(x)=3x3+2x2+x; x=0,1.

Порядок ввода данных: 3; 3; 2; 1; 0; 0,1.

Ответ: 0,123.

Программа 24. Случай комплексного аргумента.

10 INPUT "n",A:FOR B=0 TO A:PRINT A-B;:INPUT H(B):NEXT B
20 INPUT "Re",C,"Im",D:E=0:F=0:FOR B=0 TO A:G=C*E-D*F+H(B)
30 F=C*F+D*E:E=G:NEXT B:PRINT "Re";E,"Im";F:GOTO 20

Размер: 116, I(n)

Пример: P(x)=3x3+2x2+x; x=0,1+j0,2.

Порядок ввода данных: 3; 3; 2; 1; 0; 0,1; 0,2.

Ответ: 0,007+j0,274.

3.2. Оценка устойчивости

Все корни многочлена степени n

Pn(x)=е(i)aixi-1, i=1...n

отрицательны только в том случае, если то же выполняется для многочлена степени n-1

Pn-1(x)= a1xn-1+(a2 - a3a0/a1)xn-2+ a3xn-3+(a4 - a5a0/a1)xn-4+ ...

На этом правиле основан рекуррентный алгоритм, который обычно применяется для оценки устойчивости моделей систем управления.

Программа 25. Критерий Рауса - Гурвица. В случае подтверждения устойчивости программа выдает сообщение "У.", неустойчивости - "Не". Следует выделить случай, когда на какой-либо итерации a1=0. Программа выдает номер шага i, после чего прерывается, либо (если i≥n-1) выдает результат. Если результат "Не", неустойчивость доказана. При любом другом исходе необходимо ввести в член ai малое возмущение ai Ь ai + d и провести счет заново.

10 INPUT "n",A:FOR B=0 TO A:PRINT B;:INPUT E(B):NEXT B
20 FOR B=0 TO A:IF E(B)<0;PRINT "Не":GOTO 10
30 IF E(B)=0;PRINT B;
40 IF B>A-3 THEN 70
50 IF F(B)=0;PRINT B+1:GOTO 10
60 D=E(B)/F(B):FOR C=B TO A-3 STEP 2:G(C)=G(C)-D*H(C):NEXT C
70 NEXT B:PRINT "У.":GOTO 10

Размер: 161, E(n)

Пример 1: P(x)=2x4+3x3+x2+5x+1.

Порядок ввода данных: 4; 2; 3; 1; 5; 1.

Ответ: Неустойчив.

Пример 2: P(x)=12x5+12x4+16x3+10x2+5x+1.

Порядок ввода данных: 5; 12; 12; 16; 10; 5; 1.

Ответ: Устойчив.

3.3. Арифметические операции над многочленами

Программа 26. Перемножение многочленов.

P(x)= P1(x)P2(x)=е(i)aixi е(j)bjxj = е(k)(е(i+j=k)aibj)xk, i=0...n, j=0...m.

10 VAC :INPUT "n",A:FOR B=0 TO A:PRINT B;:INPUT F(B):NEXT B
20 INPUT "m",C:FOR D=0 TO C:PRINT D;:INPUT E:FOR B=0 TO A
30 G(A+B+D)=G(A+B+D)+E*F(B):NEXT B:NEXT D:FOR B=0 TO A+C
40 F(B)=G(A+B):G(A+B)=0:PRINT B;F(B):NEXT B:A=A+C:GOTO 20

Размер: 147, G(2*n+m)

Пример: P(x)=(4+3x+2x2+x3)(x+2x2)(x+2x2).

Порядок ввода данных: 3; 4; 3; 2; 1; 2; 0; 1; 2.

Промежуточный результат: 4x+11x2+8x3+5x4+2x5.

Продолжение ввода данных: 2; 0; 1; 2.

Ответ: P(x)=4x2+19x3+30x4+21x5+12x6+4x7.

Программа 27. Деление многочленов.

Pn(x) = Tn-m(x) Qm(x) +Rm-1(x);

где делимое Pn(x) и делитель Qm(x) - исходные данные,
а частное Tn-m(x) и остаток Rm-1(x) - результаты.

10 INPUT "m",A:FOR B=0 TO A:PRINT
  A-B;:INPUT E(B):NEXT B
20 INPUT "n",C:FOR B=0 TO C:PRINT C-B;:INPUT F(A+B):NEXT B
30 FOR B=A TO C:F(B)=F(B)/E:PRINT C-B;F(B)
40 FOR D=1 TO A:F(B+D)=F(B+D)-E(D)*F(B):NEXT D:NEXT B
50 FOR B=C-A+1 TO C:PRINT "R";C-B;F(A+B):NEXT B

Размер: 166, F(m+n)

Пример: P5(x)=2x5+5x4+8x3+11x2+4x; Q2(x)=2x2+x.

Порядок ввода данных: 5; 2; 5; 8; 11; 4; 0; 2; 2; 1; 0.

Ответ: T(x)=x3+2x2+3x+4; R(x)=0.

Программа 28. Общий множитель двух многочленов. Два многочлена P1(x) и P2(x) могут быть представлены в виде

P1(x)=R(x)Q1(x);
P2(x)=R(x)Q2(x);

где R(x) - общий множитель;
Q1(x), Q2(x) - взаимно простые многочлены.
R(x) определяется по алгоритму Евклида.

Для приближенных вычислений следует изменить условие в строке 80. (Например, при использовании условия ABS H(A+C-E)<1E-5 все коэффициенты, по модулю меньшие 10-5, будут считаться нулевыми.)

10 INPUT "m",A:FOR B=0 TO A:PRINT A-B;:INPUT G(B):NEXT B
20 INPUT "n",C:FOR B=0 TO C:PRINT C-B;:INPUT H(A+B):NEXT B
30 IF A>C;E=C:GOTO 80
40 FOR B=A TO C:H(B)=H(B)/G
50 IF A>0;FOR D=1 TO A:H(B+D)=H(B+D)-G(D)*H(B):NEXT D
60 NEXT B:E=A
70 E=E-1:IF E<0;FOR B=0 TO A:PRINT A-B;G(B)/G:NEXT B:END
80 IF H(A+C-E)=0 THEN 70
90 FOR B=0 TO E:F=H(A+C):FOR D=0 TO E-B:H(A+C-D)=G(A+C-D)
100 NEXT D:FOR D=A+B+1 TO 0 STEP -1:G(D)=F(D)
110 NEXT D:NEXT B:C=A:A=E:GOTO 30

Размер: 298, H(m+n)

Пример: P1(x)=2x3-3x+45; P2(x)=3x3+10x2-8x-33.

Порядок ввода данных: 3; 2; 0; -3; 45; 3; 3; 10; -8; -33.

Ответ: R(x)=x+3.

3.4. Интегрирование и дифференцирование

Программа 29.

Первообразная многочлена:

те(i)aixidx=е(i)ai-1/i xidx +C;

где новые коэффициенты ai=ai-1/i, i=1...n+1
C - произвольный свободный член.

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

Производная от многочлена:

(е(i)aixidx)' = е(i)(i+1)ai+1xidx +C;

где новые коэффициенты ai= iai+1, i=0...n-1

Управляющие ключи программы:
I - интегрировать;
D - дифференцировать;
P - показать результат;
N - новый многочлен;
T - закончить.

10 INPUT "n",A:FOR B=A TO 0 STEP -1:PRINT B;:INPUT C(B):NEXT B
20 INPUT "Операция",$:IF $="T";END
30 IF $="D";A=A-1:FOR B=0 TO A:C(B)=(B+1)*D(B):NEXT B
40 IF $="P";FOR B=A TO 0 STEP -1:PRINT B;" ";C(B):NEXT B
50 IF $="N" THEN 10
60 IF $¦"I" THEN 20
70 A=A+1:FOR B=A TO 1 STEP -1:C(B)=B(B)/B:NEXT B
80 INPUT "C",C:GOTO 20

Размер: 202, C(n+)

Пример: P(x)=9x2+8x+7;

Операция = I; C= 6; Операция = P.

Результат: 3x3+4x2+7x+6.

Операция = D; Операция = P.

Результат: 9x2+8x+7.

3.5. Замена переменной

Программа 30. Линейная замена переменной.

Исходный многочлен: Qn(y)=е(i)qiyi, i=0...n

После замены y=ax+b:

P(x)=Q(ax+b)=е(k)ak(е(i)i!/k!/(i-k)!bi-kqi)xk, k=0...n, i=k...n.

Коэффициенты нового многочлена:

pk=akе(i)i!/k!/(i-k)!bi-kqi

10 INPUT "a",A,"b",B,"n",C:FOR D=0 TO C:PRINT D;:INPUT H(D)
20 NEXT D:FOR D=0 TO C:F=0:G=1:FOR E=D TO C:F=F+G*H(E)
30 G=G*B*(E+1)/(E+1-D):NEXT E:PRINT D;F*A-D:NEXT D

Размер: 113, H(n)

Пример: y=2x+3; Q(y)=5+6y+7y2+8y3+9y4.

Порядок ввода данных: 2; 3; 4; 5; 6; 7; 8; 9.

Ответ: P(x)=1031+2472x+2260x2+928x3+144x4.

3.6. Разложение на сомножители

Квадратное уравнение:

a2x2+a1x+a0=0.

Корни:

x1,2=-a1/2a2 ¦Ц((-a1/2a2)2-a0/a2)

Вычисления в режиме калькулятора:

Пример 1: 25x2+12x+31=0.

A=25 [EXE]
B=-12/2/A [EXE]
B*B-31/A [EXE]
Вывод: -1.1824
SQR -[ANS][EXE]
Вывод: 1.08738217
B [EXE]
Вывод:-0.24

Ответ: -0,24¦j1,08738217.

Пример 2: 25x2+112x+31=0.

A=25 [EXE]
B=-112/2/A [EXE]
B*B-31/A [EXE]
3.7776 [EXE]
A=SQR [ANS][EXE]
B-A [EXE]
Вывод: -4.18360489
B+A [EXE]
Вывод:-0.29639511

Ответ: -4,18360489; -0,29639511.

Программа 31. Кубическое уравнение.

a3x3+a2x2+a1x+a0=0

Действительный корень выделяется методом дихотомии, после чего аналитически решается оставшееся квадратное уравнение. Область поиска корней определяется из неравенства: |xi|<max(i){ai/a3}+1

10 F=1:FOR E=0 TO 3:PRINT 3-E;:INPUT
  A(E)
20 IF E>0;A(E)=A(E)/A:IF ABS A(E)+1>F;F=ABS A(E)+1
30 NEXT E:A=0
40 E=A+B:G=E*A+C:IF G*A+D<0;A=A+F
50 F=F/2:A=A-F:IF A+F¦A THEN 40
60 E=-E/2:G=E*E-G:PRINT A:IF G>0 THEN 80
70 PRINT "Re";E,"Im";SQR -G:GOTO 10
80 PRINT SQR G-E,SQR G+E:GOTO 10

Размер: 197, G

Пример: 100x3+99,99x2+999999,99x-100=0.

Порядок ввода данных: 100; 99,99; 999999,99; -100.

Ответ: x1=10-4; x2,3=-0,5¦j9,99875.

Программа 32. Метод Хичкока. Алгебраическое уравнение степени n:

е(i)aixi=0, i=0...n, an¦0, a0¦0

Уравнение решается численно, путем последовательного выделения квадратных трехчленов. Коэффициенты b1 и b0 квадратного трехчлена x2+b1x+b0 хранятся в ячейках C и D. Время сходимости процедуры зависит от того, насколько удачны начальные значения b1 и b0. В отдельных случаях сходимость медленная.

10 INPUT "n",A:FOR B=0 TO A:PRINT
  A-B;:INPUT J(B):NEXT B
20 D=1:E=1:C=K/J:IF A=1;PRINT -C:END
30 F=C*C/4-L/J:IF A=2 THEN 90
40 B=D:C=E:D=0:G=J:H=J:I=K-C*J:FOR E=3 TO A:F=I:I=I(E)-C*F-B*G
50 G=F:F=H:H=G-C*F-B*D:D=F:NEXT E:G=J(A)-B*G:E=H+C*F
60 D=F*F*B+E*H:E=C+(I*E-G*F)/D:D=B+(I*F*B+G*H)/D
70 IF ABS (B-D)+ABS (C-E)>(ABS D+ABS E)*1E-10
  THEN 40
80 FOR I=0 TO A-3:K(I)=K(I)-E*J(I)-D*I(I):NEXT I:F=E*E/4-D
90 G=-C/2:H=SQR ABS F:IF F>0;PRINT G-H,G+H:GOTO 110
100 PRINT "Re";G,"Im";H
110 A=A-2:IF A>0 THEN 20

Размер: 380, J(n).

Пример: x8-1=0.

Порядок ввода данных: 8; 1; 0; 0; 0; 0; 0; 0; 0; -1.

Ответ:
x1,2= ¦1; x3,4= -1/Ц2 ¦j/Ц2;
x5,6= ¦j; x7,8=1/Ц2 ¦j/Ц2.


Программы для МК-85. Степенные ряды

4. Степенные ряды

4.0. Введение

Степенным рядом называют конструкцию вида:

f(x)=е(i)aixi, i=0...г

где ai - коэффициенты, для которых задано правило вычисления.

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

Программирование рядов проводится по рекуррентной схеме: вычисляется частная сумма ряда, представляющая собой степенной многочлен

е(i)aixi, i=0...n

где n зависит, с одной стороны, от требуемой точности, с другой - от быстроты сходимости частных сумм к истинному значению ряда.

В программе переменные резервируются под новое и старое значения частной суммы, под текущий номер i и под ту часть выражения для члена ряда, которая содержит факториал и переменную степень. В цикле, организованном по возрастанию i, две рекуррентные формулы: одна мультипликативная, для нового члена ряда; другая аддитивная, для частной суммы. Условие выхода из цикла - равенство новой и старой частных сумм. Оно выполняется тогда, когда новый член на 12 порядков меньше суммы. Обычно для этого бывает достаточно просуммировать всего 6-10 членов ряда. Если сходимость плохая и алгоритм выполняет большее число итераций, вычислительная погрешность будет существенной. В этом случае для высокоточных расчетов следует применять специальные алгоритмы (см. программу 65).

Можно также предварительно вычислить значимые коэффициенты степенного ряда (пример - программа 40), сохранить их в явном виде и затем пользоваться процедурой вычисления многочлена, реализованной в программе 23.

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

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

4.1. Вычисление систем ортогональных многочленов

Программа 33. Многочлены Чебышева. Рекуррентная формула:

Ti+1(x)=2xTi(x)-Ti-1(x), где T0(x)=1, T1(x)=x.

10 INPUT "n",A:D=1:E=0:F=1:FOR B=2 TO A:C(B+B)=0:FOR C=0 TO B
20 IF C>0;C(B+B+C)=2*B(B+C)
30 IF C<B-1;C(B+B+C)=C(B+B+C)-D(C)
40 NEXT C:FOR C=0 TO B+B:D(C)=C(C+B):NEXT C:NEXT B
50 FOR C=0 TO A:PRINT D(C+A):NEXT C:GOTO 10

Размер: 157, C(3*n)

Пример: n=4.

Ответ: ... T4(x)=1-8x2+8x4.

Программа 34. Многочлены Эрмита. Рекуррентная формула:

Hi+1(x)=2xHi(x)-2iHi-1(x), где H0(x)=1, H1(x)=2x.

10 INPUT "n",A:D=1:E=0:F=2:FOR B=2 TO A:C(B+B)=0:FOR C=0 TO B
20 IF C>0;C(B+B+C)=2*B(B+C)
30 IF C<B-1;C(B+B+C)=C(B+B+C)-2*(B-1)*D(C)
40 NEXT C:FOR C=0 TO B+B:D(C)=C(C+B):NEXT C:NEXT B
50 FOR C=0 TO A:PRINT D(C+A):NEXT C:GOTO 10

Размер: 165,C(3*n)

Пример: n=4.

Ответ: ... H4(x)=12-48x2+16x4.

Программа 35. Многочлены Лежандра. Рекуррентная формула:

Pi+1(x)=((2i+1)xPi(x) - iPi-1(x))/(i+1), где P0(x)=1, P1(x)=x.

10 INPUT "n",A:D=1:E=0:F=1:FOR B=2 TO A:C(B+B)=0:FOR C=0 TO B
20 IF C>0;C(B+B+C)=(B+B-1)*B(B+C)
30 IF C<B-1;C(B+B+C)=C(B+B+C)-(B-1)*D(C)
40 C(B+B+C)=C(B+B+C)/B:NEXT C:FOR C=0 TO B+B:D(C)=C(C+B)
50 NEXT C:NEXT B:FOR C=0 TO A:PRINT D(C+A):NEXT C:GOTO 10

Размер: 189, C(3*n)

Пример: n=4.

Ответ: ... P4(x)=0,375-3,75x2+4,375x4.

Программа 36. Многочлены Лагерра. Рекуррентная формула:

Li+1(x)=(2i+1-x)Li(x)-i2Li-1(x), где L0(x)=1, L1(x)=1-x.

10 INPUT "n",A:D=1:E=1:F=-1:FOR B=2 TO A:C(B+B)=0:FOR C=0 TO B
20 IF C>0;C(B+B+C)=-B(B+C)
30 IF C<B;C(B+B+C)=C(B+B+C)+(B+B-1)*C(B+C)
40 IF C<B-1;C(B+B+C)=C(B+B+C)-(B-1)*(B-1)*D(C)
50 NEXT C:FOR C=0 TO B+B:D(C)=C(C+B):NEXT C
60 NEXT B:FOR C=0 TO A:PRINT D(C+A):NEXT C:GOTO 10

Размер: 209, C(3*n)

Пример: n=4.

Ответ: ... L4(x)=24-96x+72x2-16x3+x4.

4.2. Вычисление значений ортогональных многочленов

Многочлены Чебышева имеют явное выражение через тригонометрические функции:

Tn(x)=cos n arccos x.

Вычисления в режиме калькулятора: A=n, B=x, C=COS A ACS B, C - Tn(x). Следует отметить, что эта формула более трудоемка и менее точна нежели рекуррентная.

Программа 37. Многочлены Лежандра.

10 INPUT "n",A,"x",B:C=1:D=B:FOR E=2 TO A
20 F=((E+E-1)*B*D-E*C+C)/E:C=D:D=F:NEXT E:PRINT F:GOTO 10

Размер: 74, F

Пример: n=3; x=0,1.

Ответ: P3(0,1)=-0,1475.

Программа 38. Многочлены Лагерра.

10 INPUT "n",A,"x",B:C=1:D=1-B:FOR E=2 TO A
20 F=(E+E-1-B)*D-(E-1)*(E-1)*C:C=D:D=F:NEXT E:PRINT F:GOTO 10

Размер: 80, F

Пример: n=3; x=0,1.

Ответ: L3(0,1)=4,289.

Программа 39. Многочлены Эрмита.

10 INPUT "n",A,"x",B:B=B+B:C=1:D=B:FOR E=1 TO A
20 F=B*D-2*E*C:C=D:D=F:NEXT E:PRINT C:GOTO 10

Размер: 68, F

Пример: n=3; x=0,1.

Ответ: H3(0,1)=-1,192.

4.3. Вычисление функций с помощью рядов

Программа 40. Интеграл вероятностей. Вычисление коэффициентов ряда

erf x = 2/Цpт0x
      exp(-t2)dt = 2/Цpе(n)(-1)nx2n+1/n!/(2n+1),
      n=0...г

Задача - получить численные значения коэффициентов ряда ai. Количество коэффициентов должно быть таким, чтобы обеспечить вычисление функции erf x в интервале [-xm; xm] с абсолютной погрешностью не выше e.

Переменные:

A | B | C | D | E | F | G
xm | a2n+1 | 2/Цp (-1)n/n! | n | e | x2n+1 | 2n+1

Вычисления прекращаются, когда выполнится |a2n+1xm2n+1| ≤ e.

10 INPUT "e",E,"xm",A:F=A:C=2/SQR p:D=0
20 G=D+D+1:B=C/G:PRINT G,B:IF ABS (F*B)≤E;END
30 D=D+1:F=F*A*A:C=-C/D:GOTO 20

Размер: 91, E

Пример: e =10-8; xm=1/2.

Ответ:
a1=1.128379167;
a3=-0.376126389;
a5=0.1128379167;
a7=-0.0268661706;
a9=0.0052239776;
a11=-0.0008548327;
a13=0.0001205537;
a15=-0.00001492565.

Программа 41. Интеграл вероятностей. Вычисление значений

10 INPUT A:B=A:C=A:D=0
20 E=B:D=D+1:C=-C*A*A/D:B=B+C/(D+D+1):IF E¦B THEN 20
30 PRINT B*2/SQR p:GOTO 10

Размер: 76, E

Переменные:

A | B | C | D | E
x | еn | (-1)n x2n+1/n! | n | еn-1

Пример: x=2.

Ответ: erf 2=0,995322275.

Программа 42. Интегральный синус

Si(x)=т0x sin t / t dt = е(n)(-1)nx2n+1/(2n+1)!/(2n+1), n=0...г

10 INPUT A:B=A:C=A:D=1
20 E=B:D=D+2:C=-C*A*A/D/(D-1):B=B+C/D:IF E¦B THEN 20
30 PRINT B:GOTO 10

Размер: 71, E

Пример: x=0,1.

Ответ: Si(0,1)=0,0999444611.

Программа 43. Интегральный косинус

Ci(x)= g+ln x +т0x cos (t -1)/ t dt = g+ln x + е(n)(-1)nx2n/(2n)!/(2n), n=0...г

где g =0,5772156649 - постоянная Эйлера.

10 INPUT A:B=.5772156649+LN A:C=1:D=0
20 E=B:D=D+2:C=-C*A*A/D/(D-1):B=B+C/D:IF E¦B THEN 20
30 PRINT B:GOTO 10

Размер: 84, E

Пример: x=0,1.

Ответ: Ci(0,1)=-1,727868408.

Программа 44. Совместное вычисление Ci(x) и Si(x).

10 INPUT A:B=.5772156649+LN A:C=0:D=0:E=1
20 F=C:GOSUB 40:E=-E:GOSUB 40:IF F¦C THEN 20
30 PRINT "Ci";B,"Si";C:GOTO 10
40 D=D+1:E=E*A/D:G=2*FRAC (D/2):B(G)=B(G)+E/D:RETURN

Размер: 125, G

Пример: x=0,1.

Ответ: Ci(0,1)=-1,727868408; Si(0,1)=0,0999444611.

Интегралы Френеля:

C(x)=т0xcos (p/2 t2) dt, n=0...г
S(x)=т0xsin (p/2 t2) dt, n=0...г

Программа 45. Случай x≤3. Разложение в ряд:

C(x)= е(n)(-1)n(p/2)2nx4n+1/(2n)!/(4n+1),
      n=0...г
S(x)= е(n)(-1)n(p/2)2n+1x4n+3/(2n+1)!/(4n+3),
      n=0...г

10 INPUT A:B=0:C=A:D=C:E=0:G=1:IF A>3 THEN 50
20 F=D:GOSUB 40:E=E+C/G:C=-C:GOSUB 40:D=D+C/G:IF F¦D THEN 20
30 PRINT "C";D,"S";E:GOTO 10
40 B=B+1:G=G+2:C=C*p/2*A*A/B:RETURN

Размер: 125, G

Пример: x=2.

Ответ: C(2)=0,4882534075; S(2)=0,343415678364.

Программа 46. Случай x>3. Степенной ряд не сходится, и вычисления проводятся по асимптотическим формулам:

 C(x)= 1/2+ cos(p x2/2)/p x
      A(x) - sin(p x2/2)/p x
      B(x)
S(x)= 1/2+ cos(p x2/2)/p x
      A(x) + sin(p x2/2)/p x
      B(x)

A(x)=е(n)(-1)n+1(4n+1)!!!!/(p x2)2n+1, n=0...г
B(x)=е(n)(-1)n+1(4n-1)!!!!/(p x2)2n, n=0...г

50 B=p*A*A:C=-1:D=0:E=C:F=1:G=1:I=1:MODE
  5
60 H=D:GOSUB 100:F=-F*C:D=D+F/I:GOSUB 100
70 G=-G*C:E=E-G/I:IF H¦D THEN 60
80 B=B/2:C=p*A:F=COS B/C:G=SIN B/C
90 PRINT "C";.5+D*F-E*G,"S";.5+D*G+E*F:GOTO 10
100 C=C+2:I=I*B:RETURN

Размер: 175, I

Пример: x=4.

Ответ: C(4)=0,498419978; S(4)=0,420516754.

Программа 47. Функция Лапласа (интеграл нормального распределения вероятности).

P(x) = 1/Ц(2p) т-гx exp(-t2/2)dt = 1/2 + exp(-x2/2)/Ц(2p)е(n)x2n+1/(2n+1)!!, n=0...г

10 INPUT A:B=A:C=1:D=0
20 D=D+B:C=C+2:B=B*A*A/C:IF D+B¦D THEN 20
30 PRINT .5+D/EXP (A*A/2)/SQR (p+p):GOTO 10

Размер: 79, D

Пример: x=3.

Ответ: P(3)=0,998650065.

Программа 48. Дилогарифм. Разложение в ряд в интервале (0;2).

f(x) = т1x ln(t)/(1-t) dt =е(n)(-1)i(x-1)i/i2, n=0...г

10 INPUT A:B=0:C=1:D=0
20 E=B:D=D+1:C=-C*(A-1):B=B+C/D/D:IF E¦B THEN 20
30 PRINT B:GOTO 10

Размер: 67, E

Пример: x=0,1.

Ответ: f(0,1)=1,29971424.


Программы для МК-85. Функции, заданные аналитически

5. Функции в форме выражений

5.0. Введение

5.0.1 Понятие функции

В математике под функцией понимают отображение множества A на множество B, при котором каждому элементу из A однозначно соответствует некоторый элемент из B. A называют областью определения, B - областью значений.

Функция в программировании - это конструкция языка, представляющая собой алгоритмическую процедуру, описанную отдельно от тела основной программы. В бейсике "МК 85" нет ни функций, ни каких-либо других программных модулей. Это обстоятельство создало ряд проблем на пути унификации бейсик-программ, заставило включить в их тексты служебные обозначения. Так, определяемая пользователем функция обозначается C=f{A,...,B}, где A,...,B - параметры, C - результат.

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

5.0.2 Встроенные функции бейсика "Электроники МК 85"

Арифметические:

-X - инверсия знака,
      "минус";
SGN X
      - "знак": 1 при x>0, -1 при x<0 и 0 при x=0;
ABS X - абсолютное
      значение, "модуль";
INT X - целая часть;
FRAC X - дробная часть;
RND(X,N) - корректное
      округление, оставляет N знаков после десятичной точки.

Тригонометрические и обратные тригонометрические:

SIN X - синус,
ASN X - арксинус (|x|≤1);
COS X - косинус,
ACS X - арккосинус (|x|≤1);
TAN X - тангенс,
ATN X - арктангенс.

Логарифмические, показательные, степенные:

LOG X - десятичный логарифм x>0;
LN X - натуральный логарифм x>0;
EXP X - экспонента x;
X-P - x в степени p (если p нецелое, x≥0);
SQR X - квадратный корень из x≥0.

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

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

5.0.3 Запись функций-выражений на бейсике

В простейшем случае функция задается одним оператором присвоения. Например, гиперболический тангенс может быть явно выражен через экспоненты: Y=(EXP X-EXP -X)/(EXP X+EXP -X). Он будет вычисляться вдвое быстрее, если экспоненты X и -X вычислить заранее: A=EXP X:B=EXP -X:Y=(A-B)/(A+B). Если вычислять гиперболический тангенс через его собственный степенной ряд, быстродействие увеличится еще вдвое.

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

Часто встречаются функции, для которых задано несколько выражений, справедливых на отдельных участках области определения. Например, экспоненциальная плотность вероятности:

j(x)= 0 при x<0;
j(x)= e-l x при x≥0

Она может быть описана строкой F=0:IF X>0;F=EXP (-L*X). Первое присвоение выполняется всегда. В данном случае это оправдано, так как делает компактнее описание процедуры, не увеличивая ее трудоемкости. В общем случае следовало бы записать:

50 IF X≤0;F=0
51 IF X>0;F=EXP (-L*X)

Аналогично описываются особые точки. Например, для f(x)=sin(x)/x:

50 F=1:IF X¦0;F=SIN X/X

5.0.4 Использование циклов и подпрограмм.

Любой вычислительный алгоритм может быть описан на бейсике набором операторов присвоения и условных арифметических операторов ( IF ... THEN ...), однако бейсик располагает и более мощными средствами управления, облегчающими процесс программирования.

Цикл FOR - NEXT удобно применять при табулировании функции. Пример:

FOR X=A TO B STEP S:Y=SIN X+COS X:PRINT "Y(";X;")=";Y:NEXT X.

Здесь (A,B) - интервал табулирования, S - его шаг.

Если функция зависит от нескольких параметров, их изменение описывается вложенными циклами:

50 FOR X=A TO B STEP S
60 FOR V=C TO D STEP R
70 PRINT X,V,SQR (X-1/V)
80 NEXT V
90 NEXT X

Здесь для каждого значения X вычисляется ряд значений функции SQR (X-1/V) в интервале (C,D) с шагом R. Таблица результатов будет двумерной.

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

50 X=A
60 ...
...
90 X=X+S:IF SGN S*X¦SGN S*B THEN 60

Обычно SGN S определен в алгоритме жестко и может быть исключен из условия (X≤B при S>0 и X≥B при S<0).

Чтобы вольное использование цикла не приводило к ошибкам, детально исследуйте его действие, придумав соответствующие тесты.

Отрабатывая оператор FOR - TO - STEP - NEXT, система выполняет операции: установить или переустановить цикл, продолжить цикл, сбросить цикл. Цикл устанавливается, когда выполнение программы доходит до оператора FOR. Если управляющая переменная цикла уже есть в списке циклов, новый цикл не устанавливается, а старый восстанавливается в исходное положение. (Именно поэтому разрешен выход из цикла минуя оператор NEXT.) В системную область памяти записывается имя управляющей переменной, выполняются все операторы присвоения, встречающиеся в заголовке цикла, запоминаются конечное значение и шаг. Обработка заголовка программы завершается запоминанием адреса следующего оператора, чтобы NEXT I знал, куда передать управление. Выполнение системой оператора NEXT может закончиться либо сбросом цикла, либо его продолжением. По параметру (I) система определяет, к какому циклу относится NEXT. После этого изменяется значение управляющей переменной (I=I+S) и проверяется, не превысило ли оно заданное конечное значение. Если превысило, цикл сбрасывается и управление передается следующему за NEXT оператору. (Физически сброс выражается в освобождении соответствующей строки списка циклов в системной области памяти.) Если условие завершения не выполнилось, управление передается первому оператору тела цикла.

Подпрограмма GOSUB - RETURN часто используется для того, чтобы вынести функцию за пределы алгоритма ее обработки. Это делает программу нагляднее. Если значение функции должно вычисляться более, чем в одной точке алгоритма, выделение ее в подпрограмму экономит память.

Определенный интерес представляет рекурсивный вызов подпрограммы. Вот, например, один из способов вычисления факториала:

10 INPUT N:M=1:GOSUB 20:PRINT "N!=";M:END
20 IF N>1;M=M*N:N=N-1:GOSUB 20
30 RETURN

Бейсик "МК 85" выдаст ошибку при N>8, так как уровень вложения подпрограмм не может превышать 8. Как и цикл, подпрограмма хранит свои параметры в системной области памяти. Оператор GOSUB #P,L запоминает адрес следующего за ним оператора и передает управление в программную область P в начало строки L. Оператор RETURN передает управление оператору, адрес которого записан в последней строке адресов возврата из подпрограмм и "вычеркивает" эту строку. Следовательно, нельзя выходить из подпрограммы, минуя RETURN.

5.1. Вычисление значений

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

5.1.1 Интегральные показательные функции

Программа 49. Функция En(x).

En(x) = т1гe-xt/tn dt.

Рекуррентная формула и разложение в степенной ряд:

En+1(x)=(e-x-xEn(x))/n
E0(x)=e-x/x
E1(x)= -g-ln x - е(i)(-1)ixi/i!/i, i=1...г

10 INPUT "n",A,"x",B:G=EXP -B:IF A=0;C=G/B:GOTO 60
20 E=1:C=0:D=0
30 D=D+1:F=C:E=-E*B/D:C=C-E/D:IF C¦F THEN 30
40 C=C-.5772156649-LN B:D=1
50 IF D¦A;C=(G-B*C)/D:D=D+1:GOTO 50
60 PRINT C:GOTO 10

Размер: 150, G

Пример: n=2; x=0,1.

Ответ: E2(0,1)=0,722545036.

Программа 50. Функция Ei(x).

Ei(x)= тxгe-t/t dt, x>0

Разложение в ряд:

Ei(x)= g+ln x + е(i)xi/i!/i, i=1...г

10 INPUT A:B=0:C=.5772156649+LN A:D=1
20 E=C:B=B+1:D=D*A/B:C=C+D/B:IF C¦E THEN 20
30 PRINT C:GOTO 10

Размер: 75, F

Пример: x=0,1.

Ответ: Ei(0,1)=-1,62281282.

Программа 51. Функция an(x).

an(x) = т1гtne-xt dt.

Рекурсивная формула:

an(x)=(e-x+n an-1(x))/x
a0(x)=e-x/x

10 INPUT "n",A,"x",B:C=0:E=EXP -B:D=E/B
20 IF C=A;PRINT D:GOTO 10
30 C=C+1:D=(E+C*D)/B:GOTO 20

Размер: 69, E

Пример: n=3; x=2.

Ответ: a3(2)=0,321421302.

Программа 52. Функция bn(x).

bn(x) = т-11tne-xt dt.

Рекурсивная формула:

bn(x)=( (-1)nex-e-x+n bn-1(x) )/x
b0(x)= (ex-e-x)/x

10 INPUT "n",A,"x",B:C=0:F=1:E=EXP -B:D=(1/E-E)/B
20 IF C=A;PRINT D:GOTO 10
30 C=C+1:F=-F:D=(F/E-E+C*D)/B:GOTO 20

Размер: 88, F

Пример: n=2; x=0,1.

Ответ: b2(0,1)=0,6686087.

5.1.2 Факториальная функция и ее производные

x!=G(x+1)= т0гe-tt x dt

Программа 53. Гамма-функция. Аппроксимация в интервале (0;1) многочленом восьмой степени (максимальная погрешность 3Ч10-7).

Расчетная формула:

G(x+1)= е(i)cixi; i=0...8

где ci - числовые коэффициенты (см. текст программы).

При |x|>1 действует рекуррентная формула:

G(x+1)=x G(x).

При x<0 действует формула

G(x+1)= |x|p / G|x+1| / sin(p x+p)

100 C=ABS A:B=FRAC C:D=(35868343*B-193527818)*B+482199394
110 D=(((D*B-756704078)*B+918206857)*B-897056937)*B+988205891
120 D=(D*B-577191652)*B*1E-9+1
130 B=B+1:IF B≤C;D=D*B:GOTO 130
140 IF A≥0;RETURN
150 IF FRAC A¦0;MODE 5:D=p*A/D/SIN
  (p*A):RETURN
160 D=1E4000:RETURN

Размер: 211, D

Пример: x=-0,5.

Обращение к подпрограмме:

10 INPUT A:GOSUB 100:PRINT D:GOTO 10

Ответ: G(0,5)=1,77245369.

Программа 54. Гамма-функция. Аппроксимация формулой Стирлинга (дает 5 верных знаков в интервале (-20;50). Расчетная формула:

G(x+1) ¦ Ц{(2p)/(x+21)}
      exp{(x+21)(ln(x+21)-1) +1/12(x+21)}/Х(i)(x+i),
      i=0...20

10 INPUT A:B=1:FOR C=A+1 TO A+20:B=B*C:NEXT C
20 A=EXP (C*LN C-C+1/12/C)*SQR (2*p/C)/B:PRINT A:GOTO 10

Размер: 69, C

Пример: x=-0,5.

Ответ: G(0,5)=1,77245371.

Программа 55. Дигамма-функция.

Y(x) = d ln G(x)/dx

Усеченный ряд

Y(x)=ln x+1/2/x + 1/12/x2+1/120/x4

обеспечивает хорошую точность только при x>20. При x<20 используется рекуррентная формула

Y(x-1) = Y(x) - 1/x

10 INPUT A:B=20+FRAC A:IF B<A;B=A
20 C=LN B+1/2/B-1/12/B/B+1/120/B-4
30 IF A>B;PRINT C:GOTO 10
40 C=C-1/B:B=B-1:GOTO 30

Размер: 88, C

Пример: x=1.

Ответ: Y(1)=-0,57721562.

Программа 56. Тригамма-функция.

Y'(x)=dY(x) / dx

Усеченный ряд для x>20:

Y'(x)=1/x - 1/2/x2+1/6/x3 - 1/30/x5

Рекуррентная формула для x<20:

Y'(x+1)= Y'(x)+1/x2.

10 INPUT A:B=20-FRAC A:IF B<A;B=A
20 C=1/B-1/2/B/B+1/6/B-3-1/30/B-5
30 IF A>B;PRINT C:GOTO 10
40 C=C+1/B/B:B=B-1:GOTO 30

Размер: 91, C

Пример: x=0.

Ответ: Y'(0)=1,64493406.

Программа 57. Неполная Гамма-функция. Разложение в ряд:

g(a, x)=e-xx a/a е(i)xia!/(a+i)!, i=0...г

10 INPUT "a",A,"x",B:C=1:D=A:E=1
20 F=C:D=D+1:E=E*B/D:C=C+E:IF F¦C THEN 20
30 PRINT C/EXP B*B-A/A:GOTO 10

Размер: 79, F

Пример: a =0,1; x=10.

Ответ: g(0,1; 10)=9,51350196.

5.1.3 Функции Бесселя

Программа 58. Функции Бесселя дробного порядка. Расчетные формулы:

Jv(x) = (x/2)v е(i)(-x2/4)v/(v+i)!, i=0...г
Yv(x) = (Jv(x)cos(vp) - J-v(x))/sin(vp)

В окрестности целых значений n погрешность вычисления Jv(x) не ограничена. Если порядок отличается от целого менее, чем на 10-4, пользуйтесь программой 59.

10 INPUT "n",A,"x",E:E=E/2:GOSUB 100:GOSUB 40
20 PRINT "J";F:G=F:A=-A:GOSUB 100:GOSUB 40:B=-A*p
30 PRINT "Y";(G*COS B-F)/SIN B:GOTO 10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=-D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN

100 Программа 53

Размер: 363, G

Пример: n=2,1; x=3.

Ответ: J2.1(3)=0,4761627444; Y2.1(3)=-0,205851878.

Программа 59. Функции Бесселя натурального порядка.

10 INPUT "n",A,"x",E:E=E/2
20 C=E-A:G=0:D=2*LN E+1.1544313298:IF A=0
  THEN 50
30 F=1/A:B=0
40 B=B+1:D=D-1/B:C=C/B:G=G-F:IF B<A;F=F*E*E/B/(A-B):GOTO 40
50 B=0:G=G/C:F=0
60 F=F+C:G=G+D*C:IF F+C=F;PRINT "J";F,"Y";G/p:GOTO
  10
70 B=B+1:C=-C*E*E/B/(B+A):D=D-1/B-1/(B+A):GOTO 60

Размер: 229, G

Пример: n=2; x=3.

Ответ: J2(3)=0,4860912606; Y2(3)=-0,16040039.

Программа 60. Модифицированные функции Бесселя дробного порядка. Расчетные формулы:

В окрестности целых значений n погрешность вычисления Kv(x) не ограничена. Если порядок отличается от целого менее, чем на 10-4, пользуйтесь программой 61.

Iv(x) = (x/2)v е(i)(x2/4)v/(v+i)!/i!, i=0...г
Kv(x) = p/2 (I-v(x) - Iv(x))/sin(vp)

10 INPUT "n",A,"x",E:E=E/2:GOSUB 100:GOSUB 40
20 PRINT "I";F:G=F:A=-A:GOSUB 100:GOSUB 40
30 PRINT "K";p/2/SIN (A*p)*(G-F):GOTO 10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN
100 Программа 53

Размер: 360, G

Пример: n=2,1; x=3.

Ответ: I2.1(3)=2,086691962; K2.1(3)=0,06513768076.

Программа 61. Модифицированные функции Бесселя натурального порядка.

10 INPUT "n",A,"x",E:E=E/2
20 C=E-A:G=0:D=2*LN E+1.1544313298:IF A=0
  THEN 50
30 F=1/A:B=0
40 B=B+1:D=D-1/B:C=C/B:G=G+F:IF B<A;F=-F*E*E/B/(A-B):GOTO 40
50 B=0:G=G/C:F=0
60 F=F+C:G=G-D*C*-1-A:IF F+C=F;PRINT "I";F,"K";G/2:GOTO
  10
70 B=B+1:C=C*E*E/B/(B+A):D=D-1/B-1/(B+A):GOTO 60

Размер: 234, G

Пример: n=2; x=3.

Ответ: I2(3)=2,245212441; K2(3)=0,0615104577.

Программа 62. Сферические функции Бесселя.

jn(x)= Ц(p/2/x) Jn+0.5(x);
yn(x)= Ц(p/2/x) Jn+0.5(x);

10 INPUT "n",A,"x",E:E=E/2:A=A+.5:MODE 5:GOSUB 70:GOSUB 40
20 G=F:A=-A:GOSUB 70:GOSUB 40:C=SQR (p/E)/2
30 PRINT "j";G*C,"y";-1-(.5-A)*F*C:GOTO 10
40 D=E-A/D:C=0:F=D
50 B=F:C=C+1:D=-D*E*E/C/(A+C):F=F+D:IF B¦F THEN 50
60 RETURN
70 D=SQR p:C=ABS A:FOR B=.5 TO C:D=D*B:NEXT B
75 IF A<0;D=-1-(A+.5)*p*C/D
80 RETURN

Размер: 232, G

Пример: n=2; x=3.

Ответ: j2(3)=0,298637497; y2(3)=-0,267038335.

5.1.4 Функции Кельвина

Программа 63. Функции Кельвина bervx и beivx.

 bervx = е(i)cos(3pv/4+p i/2)/(v+i)!/i!(x/2)v+2i,
      i=0...г
beivx = е(i)sin(3pv/4+p i/2)/(v+i)!/i!(x/2)v+2i,
      i=0...г

10 INPUT "n",A,"x",E:E=E/2:MODE 5
20 GOSUB 100:D=E-A/D:C=0:F=0:G=0:H=D
30 IF ABS D≤1E-11*ABS H;PRINT "ber";G,"bei";F:GOTO 10
40 FOR B=0 TO 1:F(B)=F(B)-D*SIN (p*(A/4+(C+B)/2)):NEXT B
50 C=C+1:D=D*E*E/C/(C+A):GOTO 30
100 Программа 53

Размер: 372, H

Пример: n=2; x=3.

Ответ: ber23=0,808368452; bei23=-0,891022351.

Программа 64. Функции Кельвина bernx и beinx, kernx и keinx.

Область P0:

10 INPUT "n",A,"x",E:E=E/2:MODE 5:D=E-A:FOR B=1 TO A
20 D=D/B:NEXT B:C=0:F=0:G=0:H=D
30 IF ABS D≤1E-11*ABS H;PRINT "ber";G,"bei";F:GOTO #1
40 FOR B=0 TO 1:F(B)=F(B)-D*SIN (p*(A/4+(C+B)/2)):NEXT B
50 C=C+1:D=D*E*E/C/(C+A):GOTO 30

Область P1:

10 H=p/4*G+LN E*F:F=p/4*F-LN E*G:I=0:J=0:K=0:L=0
20 C=1:G=1:IF A>1;FOR D=2 TO A:C=C*D-C:G=G*D:NEXT D
30 D=1:M=.75*A*p:B=0:GOSUB 120
40 R=COS M:S=SIN M:T=(Q+P)/G*D:K=K+R*T:L=L+S*T:IF B=0;U=T
50 IF B<A;I=I+R*C*D:J=J+S*C*D:GOTO 70
60 IF U-T=U THEN 100
70 B=B+1:D=D*E*E:G=G*B*(A+B):GOSUB 120:M=M+p/2
80 IF B<A;C=C/B:IF B<A-1;C=C/(A-B-1)
90 GOTO 40
100 G=E-A:F=F+I/G/2+K*G/2:H=L*G/2-J/G/2-H
110 PRINT "ker";F,"kei";H:GOTO #0
120 N=B+1:GOSUB 130:Q=P:N=N+A:GOSUB 130:RETURN
130 P=.577215664:IF N=1;RETURN
140 FOR O=1 TO N-1:P=P+1/O:NEXT O:RETURN

Размер: 589, U

Пример: n=3; x=1

Ответ: ber31=0,0156287683;
bei31=-0,0137879837;
ker31=4,0858681054;
kei31=-6,8922228687.

5.1.5 Гипергеометрические функции

Программа 65. Гипергеометрическая функция.

F(a,b,c,x)=1+ е(i)Х(j)x(a+j)(b+j)/(1+j)/(c+j), |x|<1< i=0...г, j=0...i-1.

Ряд сходится медленно. При плохой сходимости становится значительной вычислительная погрешность. Для точных вычислений используется рекурсивная формула, позволяющая скомпенсировать ошибку:

F(a,b,c,x)=y0.
yj=yj+1 x(a+j)(b+j)/(1+j)/(c+j)+1,
      j=0...n-1;
yn=lg n;

Расчет начинается с последнего члена, поэтому для достижения требуемой точности следует произвести вычисления при различных n, наблюдая за устойчивостью результата.

10 INPUT "a",A,"b",B,"c",C,"x",D,"n",F:E=LOG F
20 FOR F=F-1 TO 0 STEP -1:E=E*(A+F)/(C+F)*(B+F)/(1+F)*D+1
30 NEXT F:PRINT E:GOTO 10

Размер: 96, F

Пример: a=1; b=1; c=2; x=0,5; n=30.

Ответ: F(1; 1; 2; 0,5)=1,386294361.

Программа 66. Вырожденная гипергеометрическая функция. Разложение в ряд:

F(a,c,x)=1+ е(i)Х(j)x(a+j)/(1+j)/(c+j), |x|<1< i=0...г, j=0...i.

10 INPUT "a",A,"c",B,"x",C:D=1:E=1:F=0
20 G=D:E=E*(A+F)/(B+F)/(1+F)*C:D=D+E:IF D¦G;F=F+1:GOTO 20
30 PRINT D:GOTO 10

Размер: 93, G

Пример: a=1; c=1; x=0,5.

Ответ: F(1; 1; 0,5)=1,64872127.

5.1.6 Эллиптические функции

Программа 67. Полные эллиптические интегралы.

K(m) = т0p/21/Ц(1-sin2q) dq;
E(m) = т0p/2Ц(1-sin2q) dq

Разложение в ряды:

K(m) = p/2{1+(1/2)2m+(1/2Ч3/4)2m2+(1/2Ч3/4Ч5/6)2m3+...};
E(m) = p/2{1 - (1/2)2m - (1/2Ч3/4)2m2/3 - (1/2Ч3/4Ч5/6)2m3/5 -...};

10 INPUT A:B=-1:C=1:D=1:E=1:F=1
20 G=F:B=B+2:C=C*B/(B+1):D=D*A
30 E=E+C*C*D:F=F-C*C*D/B:IF G¦F THEN 20
40 PRINT "K";E*p/2,"E";F*p/2:GOTO 10

Размер: 115, G

Пример: m=0,5.

Ответ: K(0,5)=1,85407452; E(0,5)=1,35064387.

Программа 68. Эллиптические функции Якоби sn(u,m), cn(u,m), dn(u,m). Для вычисления используется процесс арифметикогеометрического среднего.

Начальные условия: a0=1; b0=Ц(1-m); c0=Ц m.

Пока ai¦bi,
  повторять: i=i+1;
ai=(ai-1+bi-1)/2, bi=Ц(ai-1bi-1),
  ci=(ai-1 - bi-1)/2,

Обратный ход: n=i; jn=2nanu;
ji-1=(arcsin(ci/ai
  sin ji)+ji)/2;
  i=n...1.

Результат: sn(u,m)=sin j0; cn(u,m)=cos j0; dn(u,m)= cos j0/ cos (j1 - j0)

10 MODE 5:A=1:INPUT "m",B:B=SQR (1-B):D=0
20 D=D+1:C=A+B:E(D)=(A-B)/C:B=SQR (A*B):A=C/2:IF A>B THEN 20
30 E=D:INPUT "u",C:C=2-D*A*C
40 B=C:C=(ASN (E(E)*SIN B)+B)/2:E=E-1:IF E>0 THEN 40
50 PRINT "sn";SIN C,"cn";COS C,"dn";COS C/COS (C-B):GOTO 30

Размер: 179, L

Пример: m=1-10-9; u=2.

Ответ: sn(2;1)=0,964027615; cn(2;1)=0,265802135; dn(2;1)=0,265802152.

Программа 69. Тэта-функция Якоби. Для вычисления используется процесс арифметико-геометрического среднего (см. описание программы 68). Результат:

X(u,m)=(cos(j1 - j0)/an/cos j0)1/2(1-m)1/4 Х(j)cos(2ji-1- ji)-(-1/2i+1), i=1...n

10 MODE 5:A=1:INPUT "m",B:B=SQR
  (1-B):F=B:D=0
20 D=D+1:C=A+B:G(D)=(A-B)/C:B=SQR (A*B):A=C/2:IF A>B THEN 20
30 E=D:G=1:INPUT "u",C:C=2-D*A*C
40 B=C:C=(ASN (G(E)*SIN B)+B)/2
50 G=SQR (COS (C+C-B)*G):E=E-1:IF E>0 THEN 40
60 G=COS (B-C)/A/COS C*F/G:PRINT SQR G:GOTO 30

Размер: 195, N

Пример: m=1-10-9; u=2.

Ответ: X(2;1)=0,04880075718.

Программа 70. Дзэта-функция Якоби. Для вычисления используется процесс арифметико-геометрического среднего (см. описание программы 68). Результат:

Z(u,m) = е(i)ci sin ji, i=1...n

10 MODE 5:A=1:INPUT "m",B:B=SQR
  (1-B):D=0
20 D=D+2:C=A+B:E(D)=(A-B)/C:F(D)=(A-B)/2
30 B=SQR (A*B):A=C/2:IF A>B THEN 20
40 E=D:INPUT "u",C:C=2-(D/2)*A*C:F=0
50 B=C:C=SIN B:F=F+F(E)*C
60 C=(ASN (E(E)*C)+B)/2:E=E-2:IF E>0 THEN 50
70 PRINT F:GOTO 40

Размер: 190, T

Пример: m=1-10-9; u=2.

Ответ: Z(2;1)=0,793784783.

Программа 71. Неполные эллиптические интегралы.

F(j,a) = т0j1/Ц(1- sin2a sin2q) dq;
E(j,a) = т0jЦ(1- sin2a sin2q) dq

Процесс арифметико-геометрического среднего.

Начальные условия: a0=1; b0=cos a; c0=sin a; j0=j.

Пока ai¦bi, повторять:
ji+1= Arctg(tg ji bi/ai)+ji;
ai+1=(ai+bi)/2, bi+1=Ц(aibi), ci+1=(ai - bi)/2, i=i+1

Результат: n=i;
K(a)= p/2/an;
E(a)=K(a)(1-
  е(i) 2i-1ci2);
  i=0...n
F(j,a)= jn/2n/an;
Z(j,a) = е(i)ci
  sin ji, i=1...n
E(j,a)= Z(j,a)+
  F(j,a)E(a)/K(a)

10 MODE 5:INPUT "a",B,"ф",C:E=SIN B-2:A=1:B=COS B:G=1:F=0
20 C=ATN (B/A*TAN C)+C+p*INT (C/p+.5):G=G+G:D=(A-B)/2
30 E=E+G*D*D:F=F+D*SIN C:D=A+B:B=SQR (A*B):A=D/2
40 IF A>B THEN 20
50 C=C/G/A:A=p/2/A:B=A-A*E/2:F=F+B/A*C
60 PRINT "K";A,"E";B,"Fv";C,"Ev";F:GOTO 10

Размер: 209, G

Пример: a=80*p/180; j =85*p/180.

Ответ: K(80)=3,153385252; E(80)=1,040114396;
F(85;80)=2,669350448; E(85;80)=1,024363932.
Продолжить
К Оглавлению


Программы для МК-85. Функции, заданные численно

6. Функции в форме таблиц

6.0. Введение

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

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

При интерполяции кривая проходит точно через указанные точки таблицы, при аппроксимации - вблизи точек, но так, что суммарное отклонение в некотором смысле минимально. Аппроксимацию применяют тогда, когда необходимо получить общее выражение для большого множества точек. В библиотеке представлены: интерполяция степенными многочленами, сплайнами и рядами Фурье; алгебраическая и гармоническая аппроксимация по методу наименьших квадратов.

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

Интерполяция рядами Фурье применяется к периодическим функциям, для которых обеспечивает минимальную погрешность на всей оси X.

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

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

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

2. Если границы интервала интегрирования совпадают с узлами таблицы, можно подобрать параметры интегратора так, чтобы для вычисления требовались только точки-узлы таблицы. На эту операцию ориентированы формулы Ньютона - Котеса (программы 86, 87), отличающиеся постоянным шагом интегрирования.

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

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

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

6.1. Численное дифференцирование

Программа 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.

6.2. Интерполяция степенными многочленами

6.2.1 Вычисление значений интерполяционного многочлена

Интерполяционная формула Лагранжа:

f(x)=е(i)yiХ(j)(x-xj)/(xi-xj); i=0...n, j=0...n.

где {xi, yi} - узлы интерполяции.

Программа 94. Общий случай.

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.

6.2.2 Вычисление коэффициентов интерполяционного многочлена

Интерполяционная формула Ньютона:

f(x)=r0+(x-x0){r1+(x-x1){r2+ ... (x-xn-1) rn} ... };

где ri - разделенная разность порядка i.

Программа 96. Общий случай.

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.

6.2.3 Метод Чебышева. Интерполяционный многочлен

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(i p(2k+1)/(n+1)/2), k=0...n

Программа 98. Аппроксимация методом Чебышева. Значение интерполяционного многочлена для произвольной точки x внутри интервала [a,b] определяется алгоритмом:

t=(2x-b-a)/(b-a);
bk=ck-bk+2+2t bk+1,
      k=n...1;
f(x)=c0/2-b2-t b1.

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.

6.3. Сплайн-интерполяция

6.3.1 Локальный кубический сплайн

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.

6.3.2 Глобальный кубический сплайн

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.


Программы для МК-85. Дифференциальные уравнения

7. Дифференциальные уравнения

7.0. Введение

Дифференциальное уравнение (ДУ) выражает зависимость между функцией и ее производными. Решением ДУ является функция. Численные методы позволяют получить функцию в виде таблицы.

В библиотеке рассмотрены только обыкновенные ДУ и их системы. Уравнения в частных производных, как правило, требуют более мощной вычислительной машины.

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

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

7.1. Уравнение первого порядка

y'=f(x,y);

начальные условия: (x0,y0).

Программа 110. Метод Рунге - Кутта 4-го порядка.

k1=sf(xi, yi);
k2=sf(xi+s/2, yi+k1/2);
k3=sf(xi+s/2, yi+k2/2);
k4=sf(xi+s, yi+k3);
yi+1=yi+(k1+2k2+2k3+k4)/6;

где s - шаг интегрирования.

10 INPUT "s/2",A,"x",B,"y",C
20 D=C:GOSUB 50:E=F:B=B+A:D=C+F:GOSUB 50:E=E+F+F
30 D=C+F:GOSUB 50:E=E+F+F:B=B+A:D=C+F+F:GOSUB 50
40 C=C+(E+F)/3:PRINT "x";B,"y";C:GOTO 20
50 F=A*(...f{B,D}...):RETURN

Размер: 142+, F

Пример: y'=x/y; s/2=0,05; x0=0; y0=1.

В программе: 50 F=A*B/D:RETURN

Ответ: ... f(1)=1,41421356 ...

Программа 111. Метод Милна 6-го порядка. Решение устойчиво лишь при y'>1.

Предикатор: ypi+1=yi-5+3s/10{11(y'i+y'i-4)-14(y'i-1+y'i-3)+26y'i-2}.

Корректор: yi+1=yi-3+2s/45{7(y'i+1+y'i-3)+32(y'i+y'i-2)+12y'i-1}.

Для получения первых пяти отсчетов y и y' используется метод Рунге - Кутта 4-го порядка с шагом, уменьшенным в 10 раз.

10 INPUT "s",A,"x",B,"y",C:A=A/20:G=C:GOSUB 90
20 FOR F=.1 TO 5 STEP .1:D=C:E=Q:B=B+A:C=D+Q:GOSUB 90:E=E+Q+Q
30 C=D+Q:GOSUB 90:E=E+Q+Q:B=B+A:C=D+Q+Q:GOSUB 90:C=D+(E+Q)/3
40 GOSUB 90:IF FRAC F=0;PRINT "x";B,"y";C:G(F)=C:L(F)=20*Q
50 NEXT F:A=20*A
60 C=G+.3*(11*(Q+M)-14*(P+N)+26*O):FOR E=0 TO 9:G(E)=H(E)
70 NEXT E:B=B+A:GOSUB 90:L=H+2/45*(7*(Q+M)+32*(P+N)+12*O)
80 PRINT "x";B,"y";L:GOTO 60
90 Q=A*(...f{B,C}...):RETURN

Размер: 313+, Q

Пример: y'=x/y; s=0,1; x0=0; y0=1.

В программе: 90 Q=A*B/C:RETURN

Ответ: ... f(2)=2,23606797 ...

Программа 112. Метод Хемминга 4-го порядка.

Предикатор: ypi+1=yi+3+4s/(32y'i-y'i-1+2y'i-2).

Модификатор: ymi+1=ypi+1+112/121(yci-ypi).

Корректор: yci+1=(9yi-yi-2)/8+3s/8(y'mi+1+2y'i-y'i-1).

yi+1=yci+1-9/121(yci+1-ypi+1).

Для получения первых трех отсчетов y и y' используется метод Рунге - Кутта 4-го порядка с шагом, уменьшенным в 10 раз.

10 INPUT "s",A,"x",B,"y",C:A=A/20:G=C:GOSUB 90
20 FOR F=.1 TO 3 STEP .1:D=C:E=M:B=B+A:C=D+M:GOSUB 90:E=E+M+M
30 C=D+M:GOSUB 90:E=E+M+M:B=B+A:C=D+M+M:GOSUB 90:C=D+(E+M)/3
40 GOSUB 90:IF FRAC F=0;PRINT "x";B,"y";C:G(F)=C:J(F)=20*M
50 NEXT F:A=20*A:D=0
60 F=G+4/3*(M+M-L+K+K):C=F+112/121*D:FOR E=0 TO 5:G(E)=H(E)
70 NEXT E:B=B+A:GOSUB 90:J=(9*I-G)/8+3/8*(M+L+L-K):D=J-F
80 J=J-9/121*D:PRINT "x";B,"y";J:GOSUB 90:GOTO 60
90 M=A*(...f{B,C}...):RETURN

Размер: 334+, M

Пример: y'=x/y; s=0,1; x0=0; y0=1.

В программе: 90 M=A*B/C:RETURN

Ответ: ... f(2)=2,23606797 ...

7.2. Уравнение второго порядка

y"=f(x,y,y');

начальные условия: (x0,y0,y'0).

Программа 113. Метод Рунге - Кутта 4-го порядка.

k1=sf(xi, yi,
      y'i);
k2=sf(xi+s/2, yi+y'is/2,
      y'i+k1/2);
k3=sf(xi+s/2, yi+y'is/2+k1s/2,
      y'i+k2/2);
k4=sf(xi+s, yi+y'is+k2s/2,
      y'i+k3);
y'i+1=y'i+(k1+2k2+2k3+k4)/6
yi+1=yi+s(y'i+(k1+k2+k3)/6)

где s - шаг интегрирования.

10 INPUT "s/2",A,"x",B,"y",C,"y'",D
20 E=C:F=D:GOSUB 60:G=J:B=B+A:C=E+F*A:D=F+G:GOSUB 60:H=J
30 C=C+G*A:D=F+H:GOSUB 60:I=J:B=B+A:C=E+2*A*(F+H):D=F+I+I
40 GOSUB 60:C=E+2*A*(F+(G+H+I)/3):D=F+(G+H+H+I+I+J)/3
50 PRINT "x";B,"y";C,"y'";D:GOTO 20
60 J=A*(...f{B,C,D}...):RETURN

Размер: 211+, J

Пример: y"=-y'/x-y; s/2=0,05; x0=10-30; y0=1; y'0=0.

В программе: 60 J=-A*(D/B+C):RETURN

Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1
y=0,764548247 и y'=-0,439685328.
При интегрировании с шагом s=0,01 y(1)=0,765191315,
при s=0,001 y(1)=0,76519768.
Точное значение J0(1)=0,7651976866.

7.3. Уравнение n-го порядка

y(n)=f(x, y, y', ... , y(n-1));

начальные условия: (x0, y0,y'0, ... , y(n-1)0).

Программа 114. Метод Рунге - Кутта 4-го порядка. Обыкновенное дифференциальное уравнение любого порядка сводится к системе ДУ 1-го порядка посредством замены:

y1=y; yj+1=y'j, j=1...n-1;

где n - порядок уравнения;
y - искомая функция;
yj - переменные системы (см. описание программы 115).

Программа 114 является модификацией 115 с учетом эквивалентности yj+1=y'j. Текст на бейсике записан для уравнения второго порядка. В остальных случаях пользуйтесь представленной ниже таблицей соответствия переменных.

n | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
y0j | G(C) | H(C) | I(C) | J(C) | K(C) | L(C) | M(C) | N(C)
е kij | I(C) | K(C) | M(C) | O(C) | Q(C) | S(C) | U(C) | W(C)
y(n) | G | H | I | J | K | L | M | N
y...y(n-1) | E...F | E...G | E...H | E...I | E...J | E...K | E...L | E...M
Всего | K | N | Q | T | W | Z | 28 | 31

Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом.

10 INPUT "s",A,"x",B:FOR C=1 TO 2:PRINT C;:INPUT D(C):NEXT C
20 FOR C=1 TO 2:G(C)=D(C):I(C)=0:NEXT C
30 FOR D=0 TO 1.5 STEP .5:B=B+A*FRAC D
40 G=...f{B,E,F,...}...
50 FOR C=1 TO 2:I(C)=I(C)+A*E(C)*(1+SGN FRAC (D/1.5))
60 D(C)=G(C)+A*E(C)*(1+INT D)/2:NEXT C:NEXT D
70 FOR C=1 TO 2:D(C)=G(C)+I(C)/6:NEXT C:PRINT B;E:GOTO 20

Размер: 206+, K+

Пример: y"=-y'/x-y; s=0,1; x0=10-30; y0=1; y'0=0.

В программе: 40 G=-F/B-E

Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y=0,764548247.

7.4. Система обыкновенных дифференциальных уравнений

y'j=fj(x, y1, y2, ... , yn);

начальные условия:(x0, y01 ,y02, ... , y0n).

Программа 115. Метод Рунге - Кутта 4-го порядка.

Алгоритм в векторной форме записи:

k1=s F(xi, yi);
k2=s F(xi+s/2, yi+k1/2);
k3=s F(xi+s/2, yi+k2/2);
k4=s F(xi+s, yi+k3);
yi+1=yi+(k1+2k2+2k3+k4)/6;

где s - шаг интегрирования.

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

Имена переменных в зависимости от порядка системы

n | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
y'j | F(C) | G(C) | H(C) | I(C) | J(C) | K(C) | L(C) | M(C)
y0j | H(C) | J(C) | L(C) | N(C) | P(C) | R(C) | T(C) | V(C)
е kij | J(C) | M(C) | P(C) | S(C) | V(C) | Y(C) | Z(2+C) | Z(5+C)
y'1...y'n | G...H | H...J | I...L | J...N | K...P | L...R | M...T | N...V
y1...yn | E...F | E...G | E...H | E...I | E...J | E...K | E...L | E...M
Всего | L | P | T | X | 27 | 31 | 35 | 39

Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом.

10 INPUT "s",A,"x",B:FOR C=1 TO 2:PRINT C;:INPUT D(C):NEXT C
20 FOR C=1 TO 2:H(C)=D(C):J(C)=0:NEXT C
30 FOR D=0 TO 1.5 STEP .5:B=B+A*FRAC D
40 G=...f1{B,E,F,...}...
41 H=...f2{B,E,F,...}...
50 FOR C=1 TO 2:J(C)=J(C)+A*F(C)*(1+SGN FRAC (D/1.5))
60 D(C)=H(C)+A*F(C)*(1+INT D)/2:NEXT C:NEXT D:PRINT "x";B
70 FOR C=1 TO 2:D(C)=H(C)+J(C)/6:PRINT C;D(C):NEXT C:GOTO 20

Размер: 219+, L+

Пример:

y'1=y2
y'2=-(y2/x+y1)
s=0,1; x=10-30; y0=(1, 0)T

В программе:

40 G=F
41 H=-F/B-E

Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y1=0,764548247 и y2=-0,439685328.

Программа 116. Метод Рунге - Кутта - Мерсона. Одношаговый метод 4-го порядка с грубой оценкой локальной ошибки, благодаря которой алгоритм автоматически изменяет шаг интегрирования.

Расчетные формулы в векторной записи:

k0=s F(xi, yi);
k1=s F(xi+s/3, yi+k0/3);
k2=s F(xi+s/3, yi+k0/6+k1/6);
k3=s F(xi+s/2, yi+k0/8+3k2/8);
k4=s F(xi+s, yi+k0/2-3k2/2+2k3);
yi+1=yi+(k0+4k3+k4)/6

Управление шагом интегрирования s:
Локальная ошибка r=-2k0+9k2-8k3+k4.
Если существует j, для которого |rj|>30e, тогда s=s/2;
иначе, если для каждого j |rj|<e, тогда s=s*2;
i=i+1, x=x+s.

Здесь e - допустимая локальная ошибка.

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

10 INPUT "e",E,"s",A,"x",B
20 FOR C=1 TO 2:PRINT C;:INPUT J(C):NEXT C
30 FOR C=1 TO 2:F(C)=J(C):NEXT C:GOSUB 150:FOR C=1 TO 2
40 L(C)=A*H(C):F(C)=J(C)+L(C)/3:NEXT C:B=B+A/3:GOSUB 150
50 FOR C=1 TO 2:F(C)=J(C)+(L(C)+A*H(C))/6:NEXT C:GOSUB 150
60 FOR C=1 TO 2:N(C)=A*H(C):F(C)=J(C)+L(C)/8+.375*N(C)
70 NEXT C:D=0:B=B+A/6:GOSUB 150:FOR C=1 TO 2:P(C)=A*H(C)
80 F(C)=J(C)+L(C)/2-1.5*N(C)+2*P(C):NEXT C:B=B+A/2:GOSUB 150
90 FOR C=1 TO 2:R(C)=A*H(C):F(C)=J(C)+(L(C)+4*P(C)+R(C))/6
100 F=R(C)-2*L(C)+9*N(C)-8*P(C):IF ABS F<E;D=D+1
110 IF ABS F>30*E;B=B-A:A=A/2:GOTO 30
120 NEXT C:IF D=2;A=A+A
130 PRINT "x";B:FOR C=1 TO 2:J(C)=F(C)
140 PRINT C;J(C):NEXT C:GOTO 30
150 I=...f1{B,G,H,...}...
151 J=...f2{B,G,H,...}...
160 RETURN

Размер: 500+, T+

Имена переменных в зависимости от порядка системы

n | 2 | 3 | 4 | 5 | 6 | 7 | 8
y'j | H(C) | I(C) | J(C) | K(C) | L(C) | M(C) | N(C)
y0j | J(C) | L(C) | N(C) | P(C) | R(C) | T(C) | V(C)
k0j | L(C) | O(C) | R(C) | U(C) | X(C) | Z(1+C) | Z(4+C)
k2j | N(C) | R(C) | V(C) | Z(C) | Z(4+C) | Z(8+C) | Z(12+C)
k3j | P(C) | U(C) | Z(C) | Z(5+C) | Z(10+C) | Z(15+C) | Z(20+C)
k4j | R(C) | X(C) | Z(4+C) | Z(10+C) | Z(16+C) | Z(22+C) | Z(28+C)
y'1...y'n | I...J | J...L | K...N | L...P | M...R | N...T | O...V
y1...yn | G...H | G...I | G...J | G...K | G...L | G...M | G...N
Всего | T | 26 | 33 | 40 | 47 | 54 | 61

Изменению подлежат имена и константы, выделенные в программе нежирным шрифтом

Пример:

y'1=y2
y'2=-(y2/x+y1)
e=10-8, s=0,01; x=10-30; y0=(1, 0)T

В программе: 150 I=H:J=-H/B-G:RETURN

Контроль управления шагом (значения x): 6,25 10-4; 1,25 10-3; 2,5 10-3; 5 10-3; 0,01; 0,02; 0,04; ...

Ответ: Отсчеты функции Бесселя J0(x), в частности при x=1 y1=0,76519768 и y2=-0,440050579.


Программы для МК-85. Статистика

8.1. Обработка выборочных данных (Статистика)

8.1.1 Гистограмма

Программа 117. Участок числовой оси разбит на равные отрезки.

Разбиение задается тройкой чисел (m, x0, s) и формулой xl=x0+ls. Значение x принадлежит классу l, если xl≤x<xl+1. Алгоритм подсчитывает число попаданий nj элементов выборки {xi}, i=1...n, в соответствующие классы l=0...m, а по прерывании ввода вычисляет частоты pl=nl/n.

10 VAC:INPUT "n",A,"x0",B,"s",C
20 INPUT D:D=(D-B)/C:IF D<0;D=0
30 E=E+1:IF D>A;D=A
40 F(D)=F(D)+1:GOTO 20
50 FOR D=0 TO A:PRINT D;F(D)/E:NEXT D:GOTO 20

Размер: 109, F(n)

Пример: n=2; x0=1; s=2.

Выборка: 3; 4; 3,2; 2,5; 1,5.

Завершение ввода: [AC] RUN 50[EXE].

Ответ: p0=0,4; p1=0,6; p2=0.

8.1.2 Моменты

Программа 118. При исследовании выборки {xi} сперва вычисляют начальные моменты

ak=1/n е(i)xik;

где k - порядок момента, n - объем выборки; затем из начальных моментов выражают необходимые статистические характеристики, например:

 выборочное среднее m=a1;
смещенное СКО ss=Ц(a2-a12);
несмещенное СКО s= ssЦ(n/(n-1));
асимметрию As=(a3-3a1a2+2a13)/ss3;
эксцесс Ex=(a4-4a3a1+6a2a12-3a14)/
  ss4-3.

10 VAC
20 INPUT A:C=1:FOR B=0 TO 4:D(B)=D(B)+C:C=C*A:NEXT B:GOTO 20
30 A=E/D:B=SQR (F/D-A*A):PRINT "m";A,"s";B/SQR (1-1/D)
40 PRINT "As";(G-3*F*A+2*A*A*E)/D/B-3
50 PRINT "Ex";(H-4*G*A+6*F*A*A-3*A-3*E)/D/B-4-3:GOTO 20

Размер: 166, H

Пример: 1; 2; 3; 4; 5; [AC] RUN 30[EXE].

Ответ: m=3; s=1,58...; As=0; Ex=-1,3.

8.1.3 Корреляция

Программа 119. Коэффициент корреляции

r=(е xy-е x е y/n)/Ц{(е x2-(е x)2/n)(е y2-(е y)2/n)}

служит для оценивания по выборке {xi, yi} степени зависимости случайных величин X и Y.

10 VAC
20 INPUT "x",A,"y",B:C=C+1:D=D+A:E=E+B:F=F+A*A:G=G+B*B:H=H+A*B
30 GOTO 20
40 PRINT (H-D*E/C)/SQR ((F-D*D/C)*(G-E*E/C)):GOTO 20

Размер: 107, H

Пример:

x | 1 | 1 | 3 | 3 | 5 | 5
y | 2 | 4 | 1 | 5 | 2 | 4

Завершение ввода данных: [AC] RUN 40[EXE].

Ответ: r=0.

8.1.4 Тесты

Тест позволяет принять или отвергнуть гипотезу о принадлежности выборки к эталонной популяции.

Программа 120. Критерий согласия c2. Гипотеза: выборка {xi} принадлежит к генеральной совокупности с заданным распределением F(x).

Процедура оценивания базируется на вычислении интеграла c2-распределения. Расчетная формула при числе степеней свободы m≦30:

Fm(t)=(t/2)m/2/G(m/2+1)/et/2(1+е(k)tk/Х(l)(m+2l)), k=0...г, l=1...k.

При m>30 ряд не сходится (существует другая формула). Гипотеза о принадлежности выборки к заданному распределению отвергается, если

1-Fm(t)<a;

где a - заданный уровень значимости;
t - рассогласование.

Для вычисления t необходимо так разбить выборку на классы

l=1...k, чтобы выполнялись условия:

npl≥5, l=2...k-1 и np l≥1, l=1...k;

где n - объем выборки;
pl-вероятность попадания в l-й класс, вычисленная по гипотетическому распределению.

Тогда рассогласование

t=-n+е(l)nl2/n/pl, l=1...k

где nl - число элементов выборки, попавших в l-й класс.

Число степеней свободы c2-распределения

m=k-1-r;

где r - число параметров гипотетического распределения, значения которых были определены по той же выборке.

Если интеграл гипотетического распределения оказывается трудновычислимым, разбиение на классы можно делать "на глаз" с последующим вычислением pl и соответствующим расслоением выборки.

В случае простых распределений можно составить быстрый алгоритм, который проведет разбиение на классы из условия равенства всех pl, что и реализовано в представленной программе.

В качестве примера взято экспоненциальное распределение

F(x)=1-e-x/T.

Предполагается, что значение T=10 вычислено по той же выборке.

Разбиение проводится по правилу pl=1/k, где k=Int(n/5). Квантили экспоненциального распределения, соответствующие значениям вероятности (l-1)/k, l=1...k, находятся по формуле:

xl=-T ln (1-(l-1)/k).

Так как один параметр распределения вычислен по той же выборке, определяем степень свободы как m=k-2.

10 VAC:INPUT "n",A:B=INT (A/5):IF B>30;B=30
20 FOR C=0 TO B-1:F(C+C)=-LN (1-C/B)*10:NEXT C:FOR C=1 TO A
30 INPUT D:FOR E=B+B TO 2 STEP -2:IF D<D(E);NEXT E:STOP
40 E(E)=E(E)+1:NEXT C:D=B/A:A=-A:FOR C=1 TO B:A=A+E(C+C)-2*D
50 NEXT C:B=(B-2)/2:C=p-FRAC B:FOR D=1-FRAC B TO B
60 C=C*D:NEXT D:C=(A/2)-B/C/EXP (A/2):D=0:B=B+B
70 D=D+C:B=B+2:C=C*A/B:IF D+C¦D THEN 70
80 PRINT D

Размер: 261, E(2*B)

Пример: n=15; выборка {11; 12; 8; 7; 1; 2; 19; 6; 11; 14; 1; 2; 12; 33; 11}.

Ответ: 1-F1(2,8)=0,90573569
(при a=0,9 гипотеза будет принята, при a=0,95 - отклонена).

Программа 121. Критерий Стьюдента. Гипотеза: две генеральные совокупности, распределенные нормально с равными дисперсиями, обладают равными матожиданиями. Гипотеза отклоняется, если 1-In(t)<a;

где a - уровень значимости;
In(t) - интеграл по отрезку [-t,t] от распределения Стьюдента.

Рассогласование t определяется по формуле:

t=|1/pе(l)xi-1/qе(j)yj|Ц{pq(p+q-2)/(p+q)/(е(l)xi2-е(j)yj2)}; i=1...p, j=1...q

где p, q - объемы выборок;
xi, yi - выборочные значения.
Расчетные формулы для In(t):

n=p+q-2;
q=arctg(t/Ц n);

при четных n:

In(t)=sin q (1+1/2cos2 q+3/8cos4 q+ ... +1/2 3/4 5/6 ... (n-3)/(n-2)cosn-2 q);

при нечетных n:

In(t)=2/p{ q+sin q (cos q+2/3cos3 q+ ... +2/3 4/5 6/7 ... (n-3)/(n-2)cosn-2 q)};

10 VAC:MODE 5:FOR A=1 TO 2:INPUT "p",A(A):FOR H=1 TO A(A)
20 PRINT H;:INPUT I:C(A)=C(A)+I:E(A)=E(A)+I*I:NEXT H:NEXT A
30 A=ATN (ABS (D/B-E/C)*SQR (B*C/(B+C)/(F+G-D*D/B-E*E/C)))
40 B=B+C-2:C=COS A:D=2*FRAC (B/2):E=C-D*SGN (B-1):F=0
50 F=F+E:D=D+2:IF D<B;E=E*C*C*(D-1)/D:GOTO 50
60 F=F*SIN A:IF FRAC (B/2)¦0;F=(F+A)*2/p
70 PRINT 1-F

Размер: 244, I

Пример: p=3; xi: 0; 1; 2; q=4; yi: 2; 3; 4; 5.

Ответ: 1-I5(2,766)=0,03953282.

Программа 122. Критерий Фишера. Гипотеза: две генеральные совокупности, распределенные нормально, обладают равными дисперсиями. Гипотеза отклоняется, если

F(t,v,w)<a/2;

где a - уровень значимости;
F(t,v,w) - интеграл от распределения Фишера.

F(t,v,w)=т0tG((v+w)/2)/G(v/2)/G(w/2)(v/w)v/2tv/2-1/(1+tv/w)(v+w)/2dt;

где t=S12/S22 - рассогласование;
  S12≤S22 - выборочные дисперсии;
v=n1-1, w=n2-1 - степени свободы;
  n1, n2 - объемы выборок.

Вычисление интеграла в программе производится методом Гаусса 3-го порядка с разбиением отрезка [0; t] на пять интервалов, что обеспечивает в худшем случае (при t=1) три верных знака.

10 FOR A=1 TO 2:F=0:H=0:INPUT "n",J:FOR
  G=1 TO J:PRINT G;
20 INPUT I:H=H+I*I:F=F+I:NEXT G:A(A)=J-1:C(A)=(H-F*F/J)/A(A)
30 PRINT A;C(A):NEXT A
40 A=D/E:IF A>1;A=1/A:D=B:B=C:C=D
50 D=(B+C)/2:E=B/C:G=FRAC D:F=p-G:FOR
  G=1-G TO D-1:F=F*G
60 NEXT G:FOR H=0 TO 1:B(H)=B(H)/2:G=-FRAC B(H):F=F*p-G
70 G=G+1:IF G<B(H);F=F/G:GOTO 70
80 NEXT H:C=F*E-B:F=0:G=SQR .15:H=8:I=5:J=A/5
90 FOR K=J/2 TO A STEP J:FOR L=-1 TO 1:M=K+SGN L*F(ABS L)*J
100F=F+H(ABS L)*M-(B-1)/(1+E*M)-D:NEXTL:NEXT
  K:PRINT F*C*J/18

Размер: 363, M

Пример: n1=3; xi: 1; 2; 3 (S12=1);

n2=4; yi: 1; 3; 5; 7 (S22=0,666...).

Ответ: F(0,15; 2; 3)=0,133215831 (точны 7 знаков).


Программы для МК-85. Графы и сети

9. Графы и сети

9.0. Введение

Сети - математические объекты, обладающие дискретной и нерегулярной структурой. Структура задается графом сети - совокупностью узлов и ветвей. Ветви могут быть однонаправленными (ориентированный граф) или двунаправленными (неориентированный граф).

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

В этой главе приводится несколько примеров, демонстрирующих, как можно программировать на "МК-85" небольшие сетевые задачи.

9.1. Связывающая сеть

Программа 144. Алгоритм Прима получает из полносвязной сети с двунаправленными ветвями связную сеть минимальной стоимости.

Входные данные: n - число узлов сети; cij, (j>i) - стоимость ветви (непосредственной связи) между узлами i и j.

Выходные данные: {ik, jk}, k=1...n-1 - индексы ветвей, обеспечивающих минимум суммарной стоимости сети при сохранении ее связности.

Ограничение размерности: n≤16.

10 INPUT A:DEFM (A-1)*A/2-17
20 FOR B=1 TO A-1:FOR C=B+1 TO A:PRINT B;C;
30 INPUT I((C-2)*(C-1)/2+B):NEXT C:NEXT B:D=1
40 I=1E20:FOR B=1 TO A:IF FRAC (D/2-B)<.5 THEN 90
50 FOR C=1 TO A:F=0:IF C<B;F=(B-2)*(B-1)/2+C
60 IF C>B;F=(C-2)*(C-1)/2+B
70 IF I(F)<I;I=I(F):G=B:H=C:E=F
80 NEXT C
90 NEXT B:I(E)=1E20
100 IF FRAC (D/2-H)<.5;D=D+2-H/2:PRINT G;H;I
110 IF D<2-A-1 THEN 40

Размер: 265, I((n-1)*n/2)

Пример: n=6; C=((10, 24, 66, 91, 15), (..., 17, 22, 51, 43), (..., 52, 37, 29), (...,21, 66), (..., 81)).

Ответ: c12=10; c16=15; c23=17; c =22; c45=21.

9.2. Геодезическая сеть

Программа 145. Уравновешивание невязок углов в свободной геодезической сети. (Другие задачи уравновешивания требуют модификации программы лишь в части вычисления невязок.)

Коррелаты невязок полигонов xi определяются путем решения системы уравнений:

nixi+е(j)(1-nij)xj=180(ni-2)-е(l)uli, i=1...k; j=1...k, l=1...nk

где k - число полигонов;
ni - число точек в i-том полигоне;
nij - число смежных точек полигонов i и j;
uli - отсчет внутреннего угла (в градусах) на l-той
  точке i-го полигона.

Все отсчеты сети - равноточные. Невязки углов можно будет выразить через коррелаты полигонов по следующим формулам:

для внешней точки i-го полигона - D u=xi;
для смежной точки полигонов i и j - D u=xi-xj;
для узловой точки - D u=xi-0,5е xj, j¦i.

(Эти формулы не включены в программу. )

При вводе данных предполагается, что точки сети пронумерованы и их общее число равно m. Ввод данных осуществляется в порядке обхода точек. Для каждой точки указываются номера смежных полигонов и соответствующие внутренние углы. Если точка внутренняя, последний угол может вычисляться автоматически дополнением до 360. Для перехода к новой точке вводите i=0 (если точка внутренняя, можете вводить u=0).

Изменения точности вычислений можно добиться путем изменения порога чувствительности в строке 150.

10 VAC:INPUT "m",A,"k",B:FOR C=1 TO A:D=360:PRINT C;
20 FOR E=1 TO B:G(E)=0:NEXT E
30 INPUT "i",E:IF E=0 THEN 60
40 G(E)=1:E=E+B:INPUT "u",F:IF F>0;D=D-F:G(E)=G(E)+F:GOTO 30
50 G(E)=G(E)+D
60 FOR E=1 TO B:FOR F=E TO B:IF G(E)+G(F)<2 THEN 80
70 GOSUB 180:G(G)=G(G)+1
80 NEXT F:NEXT E:NEXT C:FOR E=1 TO B:FOR F=E TO B:GOSUB 180
90 IF E¦F;G(G)=1-G(G):GOTO 110
100 G(E+B)=(G(G)-2)*180-G(E+B):G(E)=0
110 NEXT F:NEXT E
120 $="":FOR C=1 TO B:A=G(C+B):FOR D=1 TO B:IF C=D THEN 150
130 E=C:F=D:IF F<E;E=D:F=C
140 GOSUB 180:A=A-G(G)*G(D)
150 NEXT D:F=C:GOSUB 180:A=A/G(G):IF ABS (G(C)-A)>1E-10;$=CHR 1
160 G(C)=A:NEXT C:IF $¦"" THEN 120
170 FOR A=1 TO B:PRINT A;G(A):NEXT A:END
180 G=(F-1)*F/2+E+B+B:RETURN

Размер: 475, G((k+5)*k/2)

Пример:

Ввод данных:

6 3
1 40 2 90 0
2 30 3 49 0
3 91 0
1 61 2 36 3 58 0
1 91 0
3 160 2 0

Ответ: x=(-4,11764705; -0,352941175; -0,6470588258)T

9.3. Сеть систем массового обслуживания

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

Рассмотрим конкретный пример:

Число систем в сети n=3. Узел 0 (S)- внешний источник заявок; узлы 1, 2, 3 - системы массового обслуживания. pij вероятности перехода заявки из системы i в систему j.

Программа 146. Оценка средних времен ожидания.

Интенсивности потоков заявок на входах систем

li=ai l0, i=1...n.

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

м(p11-1)a1+p21a2+...+pn1an=-p01
н p12a1+(p22-1)a2+...+pn2an=-p02
п...
о p1na1+p2na2+...+(pnn-1)an=-p0n

Допустимая размерность 2≤n≤9.

Узел сети - многоканальная система массового обслуживания с неограниченным общим буфером. Загрузка такой системы

ri=liJi/mi;

где Ji - среднее время обслуживания;
mi - число каналов обслуживания.

Если ri оказывается больше единицы, программа запрашивает другие Ji и mi.

Среднее время ожидания:

w=J(lJ)m/(m-1)!/(m-lJ)2/{(lJ)m/(m-1)!/(m-lJ)+е(k)(lJ)k/k!}

DEFM n*(n+1)-19
10 INPUT "n",A:FOR B=0 TO A:FOR C=1 TO A:PRINT B;C;
20 INPUT D:IF B=0;D=-D
30 IF C=B;D=D-1
40 G(A*B+C)=D:NEXT C:NEXT B
50 FOR B=1 TO A-1:FOR C=B+1 TO A
60 F=G(B*A+B):G=G(B*A+C):IF F=0;IF G=0 THEN 100
70 E=SQR (F*F+G*G):D=G/E:E=F/E
80 FOR F=0 TO A:G=E*G(F*A+B)+D*G(F*A+C)
90 G(F*A+C)=E*G(F*A+C)-D*G(F*A+B):G(F*A+B)=G:NEXT F
100 NEXT C:NEXT B:FOR B=A TO 1 STEP -1:D=0:IF B=A THEN 120
110 FOR C=B+1 TO A:D=D+G(A*A+C)*G(C*A+B):NEXT C
120 G(A*A+B)=(G(B)-D)/G(B*A+B):NEXT B
130 INPUT "l",C:FOR B=1 TO A:H=G(A*A+B)*C:PRINT B;"l";H
140 INPUT "y",E,"m",F:E=E*H:PRINT "r";E/F:IF E≥F THEN 140
150 I=0:G=1/E:D=0
160 D=D+G:I=I+1:IF I<F;G=G*(F-I)/E:GOTO 160
170 D=E/H/(F-E+(F-E)-2*D):PRINT "w";D:NEXT B:GOTO 130

Размер: 508, G(n*n+n)

Расчет примера.

Ввод: n=3; l0=0,5; P= (1; 0; 0); (0,1; 0,3; 0,6); (0; 0; 0,2); (0; 0; 0))
Вывод: l1=0,5556.
Ввод: J1=10; m1=6.
Вывод: r1=0,9259; w1=18,09; l2=0,1667.
Ввод: J2=7; m2=2.
Вывод: r2=0,5833; w2=3,611; l3=0,3667.
Ввод: J3=1; m3=1.
Вывод: r3=0,3667; w3=0,5789.

9.4. Граф вероятных состояний

Моделируемая динамическая система может находиться в любом из случайных состояний Si, i=1...n. Известны интенсивности lij переходов из состояния i в состояние j, а также исходное состояние Sk (pk(0)=1). Необходимо проследить вероятности состояний pi(t) во времени.

Задача сводится к системе дифференциальных уравнений Колмогорова

dpi(t)/dt = -pi(t)е(j)lij+е(j)pj(t)lij, i=1...n, j=1...n.

Программа 147. Система уравнений Колмогорова. В данном случае для решения используется метод Рунге - Кутта 4-го порядка (см. описание программы 115). Помимо данных, перечисленных выше, необходимо задать шаг интегрирования s. Для завершения списка переходов вводите i=0. В процессе счета программа высвечивает модельное время, не прерываясь. Чтобы вывести на дисплей текущий вектор вероятностей, достаточно нажать любую клавишу (за исключением [EXE]) и немного подождать.

DEFM 4*n+3*m-20
10 VAC:INPUT "n",A:A=A*4:B=A
20 INPUT "i",D:IF D≤0 THEN 40
30 B=B+3:D(B)=D*4:INPUT "j",D,"l",F(B):E(B)=D*4:GOTO 20
40 INPUT "k",C,"s",F:C(4*C)=1
50 E=E+F:PRINT E;:FOR C=4 TO A STEP 4:E(C)=C(C):F(C)=0:NEXT C
60 FOR D=0 TO 1.5 STEP .5:FOR C=4 TO A STEP 4:D(C)=0:NEXT C
70 FOR C=A+3 TO B STEP 3:D(D(C))=D(D(C))-F(C)*C(D(C)):$=$+KEY
80 D(E(C))=D(E(C))+F(C)*C(D(C)):NEXT C:FOR C=4 TO A STEP 4
90 F(C)=F(C)+F*D(C)*(1+SGN FRAC (D/1.5))
100 C(C)=C(C)+F*D(C)*(1+INT D)/2:NEXT C:NEXT D:CSR 0
110FOR C=4 TO A STEP 4:C(C)=E(C)+F(C)/6:IF $¦"";PRINT C/4;C(C)
120 NEXT C:$="":GOTO 50

Размер: 498, F(4*A+3*B)

Пример:

Ввод данных:

DEFM 4*4+3*6-20
SET 7
4
1 2 3.4
2 3 7.2
3 4 2.6
4 1 1.3
1 3 2.1
3 1 3.4
0
1 0.1

Ответ: В момент t=0,2 p=(0,5296883; 0,1313437; 0,2429235; 0,09604452)T

9.5. Герт-сеть

Ветви Герт-сети обозначают действия, характеризуемые матожиданием mi и дисперсией Di времени выполнения. Узлы отображают ветвление процесса, описываемое распределением вероятностей перехода pi. На выходе каждого узла е pi=1.

Программа 148. Правило Мейсона для спектральных коэффициентов передачи ветвей:

Wt=-(е(E)(-1)k(C)Х(C)Wi)/(1+е(L)(-1)k(C)Х(C)Wi);

где Wt - коэффициент передачи всей сети;
E - множество контуров, содержащих фиктивную обратную связь;
L - множество контуров, не содержащих фиктивную обратную связь; (фиктивная обратная
  связь - это условно вводимая ветвь, соединяющая сток сети с ее истоком);
i - номер ветви;
Wi - коэффициент передачи ветви;
C - контур (кольцо ветвей);
k(C) - порядок контура; (контур порядка k - это k
  непересекающихся простых контуров).

Подстановка в уравнение Мейсона вероятностей pi позволяет вычислить вероятность pt перехода от источника к стоку, а двукратное аналитическое дифференцирование уравнения позволяет записать аналогичные выражения для матожидания mt=f{pi,mi} и дисперсии Dt=f{pi,mi,Di}. Эти три выражения заложены в программу.

Программа хранит в памяти список ветвей:

{pi, mi, Di},i=1...n.

Ограничение размерности: n≤41.

Выделять контуры будем визуально (в противном случае программа оказалась бы слишком длинной). Список контуров зависит от того, между какими узлами определяется коэффициент передачи. Между выбранной парой узлов следует замкнуть фиктивную обратную связь. При вводе списка: в ответ на запрос программы "k" указывается порядок контура; в ответ на запрос "e" контуры, содержащие фиктивную ветвь, помечаются цифрой "1", не содержащие - цифрой "0"; затем перечисляются номера ветвей, составляющих контур. Для перехода к следующему контуру вводите i=0, для окончательного завершения ввода - k=0.

10 INPUT "n",A:DEFM A*3-14:FOR B=1 TO A:C=B*3:PRINT B;
20 INPUT "p",J(C),"m",K(C),"D",L(C):NEXT B
30 C=1:D=0:E=0:F=0:G=0:H=0
40 J=1:K=0:L=0:INPUT "k",I:IF I=0 THEN 100
50 INPUT "e",B
60 INPUT "i",A:A=A*3:IF A=0 THEN 80
70 J=J*J(A):K=K+K(A):L=L+L(A):GOTO 60
80 I=-1-I:L=(L+K*K)*J:K=K*J:FOR A=0 TO 2
90 C(A+A+B)=C(A+A+B)+I*J(A):NEXT A:GOTO 40
100 J=-D/C:K=D*E/C/C-F/C:L=D*E/C-F:L=G*D-2*E*L:L=L/C/C-H/C
110 K=K/J:PRINT "p";J,"m";K,"D";L/J-K*K:GOTO 30

Размер: 354, L(3*C)

Пример:

Ввод списка ветвей:

n=5
p1=1; m1=3; D1=6
p2=0,2; m2=1; D2=1
p3=0,1; m3=2; D3=3
p4=0,4; m4=0; D4=0
p5=0,3; m5=4; D5=0

1) Эквивалентная ветвь s - t:

Ввод списка контуров:
1 0 3 0
1 0 2 1 0
1 1 1 4 0
0
Ответ: pt=0,57; mt=4,43; Dt=15,6.

2) Эквивалентная ветвь s - q:

Ввод списка контуров:
1 0 3 0
1 0 1 2 0
1 1 1 5 0
0
Ответ: pq=0,43; mq=8,43; Dq=15,6.

9.6. Перт-сеть

Ветви Перт-сети - это действия, которые должны выполняться в заданной последовательности, включая и параллельное выполнение. Узел сети считается выполненным, если выполнены все входящие в него ветви. Задержки в ветвях описываются матожиданием mi и дисперсией si2. Алгоритм определения времени выполнения сети по методу Перт осуществляет суммирование каждой из оценок по путям от истока к стоку. В узлах, где пути сходятся, выбирается максимальная из сумм. Следует иметь в виду, что операция вычисления максимума m и s2 не имеет строгого математического обоснования, и при наличии в сети параллельных ветвей с близкими характеристиками метод заметно занижает результирующую оценку матожидания и завышает дисперсию.

Программа 149. Оценка раннего времени выполнения. Узлы должны быть пронумерованы подряд, начиная с единицы, так чтобы все ветви были направлены от младших номеров к старшим. Количество узлов не должно превышать 83. Для задания нового узла ответьте "+" на запрос "E" (программа сообщит текущий номер узла n); для завершения ввода сети ответьте E=".". Для текущего узла перечисляются все заходящие в него ветви. На запрос "r" введите номер узла, из которого исходит эта ветвь, далее введите mi и si. Для перехода к новому узлу введите r=0. Программа сообщит ожидемый момент выполнения текущего узла mn и его СКО sn.

DEFM 2k-21 (k - количество узлов в сети)
10 VAC:B=2
20 B=B+2:PRINT "n";B/2;
30 INPUT " r",A:IF A=0 THEN 70
40 INPUT "m",C,"s",D:C=C+C(A+A):IF C>C(B);C(B)=C
50 D=D*D+D(A+A):IF D>D(B);D(B)=D
60 GOTO 30
70 PRINT "m";C(B),"s";SQR D(B):INPUT "E",$:IF $¦"." THEN 20

Размер: 160, D(2*B)

Пример:

Список ветвей

r - n | mi | si
1 - 2 | 15 | 2
1 - 3 | 27 | 1
2 - 3 | 11 | 0
1 - 4 | 25 | 3
2 - 4 | 10 | 3
3 - 4 | 14 | 2
2 - 5 | 38 | 3
4 - 5 | 8.5 | 0.5

Ввод-вывод данных:

n=1
1; 15; 2
0 [+]
m1=15; s1=2
n=2
1; 27; 1
2; 11; 0
0 [+]
m2=27; s2=2
n=3
1; 25; 3;
2; 10; 3;
3; 14; 2;
0 [+]
m3=41; s3=3.6
n=4
2; 38; 3;
4; 8.5; 0.5
0 [.]
m4=53; s4=3.64

9.7. Дерево решений

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

Программа 150. Задача об инвестициях. Адаптированная версия.

Вы - владелец коммерческой фирмы. Через надежные каналы до Вас дошла информация о том, что через неделю все излишки со счетов предприятий изымут в централизованный фонд помощи перестройке. У Вас горят R рублей. Вы озадачили друзей - и начался мозговой штурм.

Когда шторм утих, Вы изучили его дары.

В общей сложности друзья, а также приятели их знакомых предложили Вам n выгодных способов пожарного вложения денег, каждый из которых требует ровно si рублей и обещает ровно через год ci дохода в твердой валюте, i=1...n. Но получился небольшой перебор - оказалось, что ваших "деревянных" не хватает на реализацию всех предложений. Приходится отбирать то, что по зубам. Разумеется, любой законопослушный гражданин будет сортировать не друзей, а предложения. Таким образом, задача может быть строго формализована:

Найти такое I={i} - подмножество {1...n}, для которого

C=е(i)ci !- max при соблюдении r=R-е(i)si≥ 0.

Решение такой задачи "в лоб" требует перебора всех комбинаций. Предположим, что n=30. Тогда всего комбинаций 230, то есть миллиард - и у Вас мало шансов получить ответ при жизни. Методы дискретной оптимизации дают рецепт того, как уменьшить число просматриваемых вариантов, не упустив оптимальный. Так, не следует просматривать нереализуемые и заведомо неперспективные направления поиска.

Предварительно расположим предложения в порядке уменьшения их рентабельности - так чтобы выполнялось:

ci/si ≥ ci+1/si+1, i=1...n-1

Такое упорядочивание обеспечивает скорейшее нахождение оптимума. Скорость зависит еще и от того, насколько точна оценка перспективности направления поиска. Здесь мы примем самую близорукую оценку - текущую сумму потерь дохода. В начале поиска примем потери субоптимального варианта большими (P=1099), а текущие потери нулевыми (p0=0).

Процесс поиска представляет собой чередование прямого и обратного хода. Прямой ход - i - i+1, обратный - i - i-1. Начинаем с i=1 и движемся вперед.

На прямом ходе будем принимать/отклонять предложения последовательно, начиная с первого. Мы можем принять текущее предложение, если на него хватает денег, то есть ri-1≥si. В этом случае остаток средств уменьшится:

ri=ri-1-si.

В противном случае будем отказываться от предложения. Но тогда увеличатся накопленные потери:

pi=pi-1+ci.

Если накопленные потери сравнялись с тем, что мы имеем для уже найденного претендента на лучший вариант,

pi≥P,

дальнейшее движение вперед в этом направлении бесперспективно - нужно возвращаться. Но этого может и не произойти. Тогда, рассмотрев последнее, n-е предложение, мы записываем новый рекорд P и запоминаем тот выбор, который к нему привел. Оптимальный выбор записывается в виде бинарного вектора x. xi=1, если i-е предложение было принято, и xi=0, если отклонено.

Разумеется, двигаясь вперед, мы запоминали текущий выбор в бинарном векторе t. Таким образом, запоминание нового рекорда сводится к следующим присвоениям:

x - t, P - pn.

Как вести себя, двигаясь обратно?

Если текущий выбор был ti=0, продолжаем двигаться обратно.

Если же ti=1, возможно, что pi<P. Если это действительно так, переключаемся на ti=0 и пытаемся двигаться вперед. Если же, двигаясь назад, мы вернулись в i=0, поиск закончен.

Ограничение памяти - n≤22. Распределение переменных:

A | B | C | D | E | F | G | H | I | J | K | L
n*6-6 | P | R | i*6-6, p0 | s1 | c1 | t1 | x1 | r1 | p1 | s2 | c2

Для пересчета задачи с новым R можете использовать команду RUN 30.

10 INPUT A:DEFM 6*A-22:FOR B=1 TO A:D=B*6-6:PRINT "c";B;
20 INPUT F(D),"s",E(D):NEXT B:A=A*6-6
30 INPUT "R",C:B=1E99:D=-6
40 D=D+6:IF D>A THEN 100
50 IF E(D)≤C(D);G(D)=1:I(D)=C(D)-E(D):J(D)=D(D):GOTO 40
60 G(D)=0:I(D)=C(D):J(D)=D(D)+F(D):IF J(D)<B THEN 40
70 D=D-6:IF D<0 THEN 110
80 IF G(D)=1 THEN 60
90 GOTO 70
100 B=D(D):FOR D=0 TO A STEP 6:H(D)=G(D):NEXT D:GOTO 70
110 B=0:FOR D=0 TO A STEP 6
120 IF H(D)>0;B=B+F(D):C=C-E(D):PRINT D/6+1;
130 NEXT D:PRINT CHR 0,"C";B,"r";C

Размер: 339, D(n*6)

Пример: n=10;

i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
ci | 87 | 62 | 21 | 81 | 50 | 87 | 46 | 33 | 14 | 11
si | 15 | 29 | 12 | 56 | 35 | 98 | 63 | 50 | 29 | 75

R=300.

Ответ: I={1,2,3,4,5,6,8}; C=421; r=5.

К Оглавлению


Программы для МК-85. Численные методы

Численные методы

5.2. Поиск нулей

Раздел содержит численные методы решения нелинейных алгебраических уравнений

F(x)=0

Для половины методов поиск корня представляет собой итерационный процесс движения по оси X от некоторой точки x0, выбранной в качестве нулевого приближения, к точному значению корня x*. Новые приближения вычисляются по рекуррентной формуле xi+1=G(xi). Для другой половины методов задается начальный интервал поиска [a,b], сужаемый затем процедурой

[ai+1,bi+1]=G([ai,bi]).

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

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

5.2.1 Прямые методы

Уравнение должно быть выражено относительно x:

x=f(x).

Процесс сходится, если в области поиска корня |f'(x)|<1.

Программа 72. Метод простых итераций:

xi+1=f(xi)

10 INPUT A
20 B=A:A=...f{B}...:IF A¦B THEN 20
30 PRINT A

Размер: 28+, B

Пример: x=sin x+1/4; x0=1,2.

В программе:

MODE 5
20 B=A:A=SIN A+.25: ...

Ответ: x*=1,17122966.

Программа 73. Метод Эйткена - Стеффенсона:

x'=f(xi); x"=f(x'); xi+1= (xix" - x'x')/(xi - 2x'+x").

10 INPUT A
20 B=A:GOSUB 50:C=A:GOSUB 50:D=A
30 E=B-C-C+D:IF ABS E<1E-5;PRINT D:GOTO 10
40 A=B*D/E-C*C/E:GOTO 20
50 A=...f{A}...:RETURN

Размер: 79+, E

Пример: x=sin x+1/4; x0=1,2.

В программе:

MODE 5
50 A=SIN A+.25:RETURN

Ответ: x*=1,17122966.

5.2.2 Методы деления отрезка

Уравнение записывается в канонической форме:

F(x)=0.

Методы деления отрезка накладывают минимум ограничений на вид F(x): если F(x) непрерывна на отрезке поиска [a,b], и Sign F(a) ¦ Sign F(b), то существует корень x внутри [a,b]. Наилучшая сходимость наблюдается при делении отрезка пополам. Требуемое число итераций

Н=log2{(b - a)/e}, где e - допустимая абсолютная ошибка.

Программа 74. Метод дихотомии.

Пока a¦b, повторять:
x=(a+b)/2; если Sign F(x)=Sign F(a), тогда a=x, иначе b=x.

10 INPUT "a",A,"b",B,"e",H:C=A:GOSUB 50:E=D:C=B:GOSUB 50:F=D
20 IF B-A<H;PRINT C:GOTO 10
30 C=(A+B)/2:GOSUB 50:G=1:IF SGN D=SGN E;G=0
40 A(G)=C:E(G)=D:GOTO 20
50 D=...f{C}...:RETURN

Размер: 120+, H

Пример: F(x)=1/(x - 1) - 1; a=-0,5; b=1,5; e=10-5.

В программе: 50 D=1/(C-1)-2-1:RETURN

Ответ: x*=-7,6 10-6.

5.2.3 Метод Ньютона и его модификации

Согласно методу Ньютона, новое приближение xi+1 определяется как пересечение оси X с прямой, которая является касательной кривой F(x) в точке xi.

Программа 75. Метод Ньютона:

x i+1= xi-F(xi)/F'(xi).

Процесс сходится, если в области поиска корня F'(x) ¦ 0 и F(x)F"(x)>0. Метод следует применять тогда, когда выражение F(x)/F'(x) вычисляется быстрее, чем дважды F(x).

10 INPUT A
20 B=A:A=B-...f{B}/f'{B}...:IF A¦B THEN 20
30 PRINT A

Размер: 30+, B

Пример: F(x)=x-sin x-1/4; (F'(x)=1-cos x); x0=1,2.

В программе:

MODE 5
20 B=A:A=A-(A-SIN A-.25)/(1-COS A): ...

Ответ: x*=1,17122966.

Программа 76. Модифицированный метод Ньютона:

xi+1=xi - D F(xi)/(F(xi+D) - F(xi))

В программе принято D=x0Ч10-6.

Модификация заключается в том, что производная берется не аналитически, а численно.

10 INPUT A:B=A*1E-6
20 GOSUB 50:D=C:A=A+B:GOSUB 50
30 D=D*B/(C-D):A=A-D-B:IF D+A¦A THEN 20
40 PRINT A:GOTO 10
50 C=...f{A}...:RETURN

Размер: 82+, D

Пример: F(x)=x-sin x-1/4; x0=1,2.

В программе:

MODE 5
50 C=A-SIN A-.25:RETURN

Ответ: x*=1,17122966.

Программа 77. Метод Рыбакова. Алгоритм осуществляет поиск корней на заданном промежутке [a,b]. Чтобы условия сходимости были обеспечены на всем интервале, F'(xi) заменена некоторой константой M, которая выбирается из условия M ≥ | F'(x)|, a<x<b.

Не следует задавать M с большим запасом, так как в этом случае поиск пойдет слишком медленно.

Чтобы процесс не увязал в точках, где F(x)=0, в алгоритм вводится интервал различимости корней e (в программе задано e =10-3).

Положить x=a.

Пока x≤b, повторять:

пока F(x)¦0, повторять: x=x+|F(x)|/M;

отпечатать x; положить x=x+e.

10 INPUT "a",A,"b",B,"M",C
20 D=...f{A}...:A=A+ABS D/C:IF A>B THEN 10
30 IF C+D=C;PRINT A:A=A+1E-3
40 GOTO 20

Размер: 70+, D

Пример: F(x)=sin(x-0,3)+0,2x; a=-1; b=1; M=1,2.

В программе:

MODE 5
20 D=SIN (A-.3)+.2*A: ...

Ответ: x*=0,249982623.

5.2.4 Интерполяционные методы

Интерполяционный метод порядка n заменяет функцию F(x) на интервале [a,b] интерполяционным многочленом степени n и вычисляет его корень аналитически, после чего конструирует новый многочлен.

Программа 78. Комбинированный метод секущих - хорд:

xi+1= xi - (xi - xi-1)/( F(xi) - F(xi-1) ) F(xi-1).

Процесс сходится, если в области поиска корня производная не меняет знак.

10 INPUT "a",A,"b",B:C=A:GOSUB 40:E=D:C=B:GOSUB
  40
20 C=B-(B-A)/(D-E)*D:IF B=C;PRINT C:GOTO 10
30 A=B:B=C:E=D:GOSUB 40:GOTO 20
40 D=...f{C}...:RETURN

Размер: 97+, E

Пример: F(x)=x-sin x-1/4; x0=1,2; x1=1,4.

В программе:

MODE 5
40 D=C-SIN C-.25:RETURN

Ответ: x*=1,171229652.

Программа 79. Метод обратной квадратичной интерполяции:

xi+1= a - F(a){(b-a)/(F(b)-F(a)) + ((b-a)/(F(b) - F(a))-(xi - a)/(F(xi) - F(a)))/ (F(xi) - F(a)) F(b)}

Здесь a и b - две неподвижные опорные точки.

Процесс сходится, если F(x) непрерывна. Сходимость - быстрая.

10 INPUT "a",A,"b",B:C=A:GOSUB 60:E=D:C=B:GOSUB 60:F=D
20 C=(A+B)/2:H=(B-A)/(F-E)
30 G=C:GOSUB 60:I=(C-A)/(D-E)
40 C=A-E*(I+D*(I-H)/(F-D)):IF C¦G THEN 30
50 PRINT C:GOTO 10
60 D=...f{C}...:RETURN

Размер: 140+, I

Пример: F(x)=x-sin x-1/4; x0=1,2; x1=3,4.

В программе:

MODE 5
60 D=C-SIN C-.25:RETURN

Ответ: x*=1,171229653.

5.2.5 Методы решения систем нелинейных уравнений

Представлены два популярных метода решения систем - метод Зейделя и метод Ньютона - Канторовича. Обе программы снабжены минимумом сервиса. Для контроля сходимости процесса приостанавливайте счет и просматривайте промежуточные результаты в калькуляторном режиме.

Программа 80. Метод Зейделя. Система записана в нормальной форме:

xj=fj(x1,...,xn), j=1...n.

Условие сходимости:

|dfj/dxj|<1, j=1...n.

Рекуррентная формула:

xj(i+1)=fj(x1(i+1),...,xj-1(i+1),xj(i),...,xn(i)), j=1...n.

Условие окончания счета:

|xj(i+1)-xj(i)| < e, j=1...n.

10 A=n:FOR B=1 TO A:PRINT B;:INPUT D(B):NEXT B
20 INPUT "e",C
30 $=CHR 1:FOR B=1 TO A:D=D(B):GOSUB 60+10*B
40 IF ABS (D(B)-D)>C;$=""
50 NEXT B:IF $="" THEN 30
60 FOR B=1 TO A:PRINT B;D(B):NEXT B:GOTO 20
70 E=...f1{E,F,...,D(A)}...:RETURN
80 F=...f2{E,F,...,D(A)}...:RETURN
90 ...
... D(A)=...fn{E,F,...,D(A)}...:RETURN

Размер: 131+, D(n)

Пример:

x1=1-exp(-x2)
x2=Ц x1
x(0)=(0.01, 0.01)T
e=10-10.

В программе:

10 A=2: ...
70 E=1-EXP -F:RETURN
80 F=SQR E:RETURN

Ответ: x1*=0,5105908269; x2*=0,7145563847.

Программа 81. Метод Ньютона - Канторовича. Система записана в канонической форме:

Fj(x1,...,xn)=0, j=1...n.

Линеаризованная система:

W(x(i))D x=F(x);

где W(x) - матрица частных производных векторной функции F(x(i)) в точке x(i).

Рекуррентная формула:

x(i+1)=x(i) - Dx(i).

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

Условие окончания счета:

|Dxj(i+1)|<e, j=1...n.

A=n
I=A*A+A
DEFM A+I-17
10 FOR B=1 TO A:PRINT B;:INPUT I(I+B):NEXT B
20 INPUT "e",H
30 FOR B=1 TO A:I(I+B)=I(I+B)+H:GOSUB 190:FOR C=1 TO A
40 I(A*C+B)=I(C):NEXT C:I(I+B)=I(I+B)-H:NEXT B:GOSUB 190
50 FOR B=1 TO A:FOR C=1 TO A:I(A*B+C)=(I(A*B+C)-I(B))/H:NEXT C
60 NEXT B:FOR B=1 TO A-1:FOR C=B+1 TO A
70 F=I(B*A+B):G=I(C*A+B):IF F=0;IF G=0 THEN 120
80 E=SQR (F*F+G*G):D=G/E:E=F/E
90 FOR F=1 TO A:G=E*I(B*A+F)+D*I(C*A+F)
100 I(C*A+F)=E*I(C*A+F)-D*I(B*A+F):I(B*A+F)=G:NEXT F
110 G=E*I(B)+D*I(C):I(C)=E*I(C)-D*I(B):I(B)=G
120 NEXT C:NEXT B
130 $=CHR 1:FOR B=A TO 1 STEP -1:D=0:IF B=A THEN 150
140 FOR C=B+1 TO A:D=D+I(I-A+C)*I(B*A+C):NEXT C
150 I(I-A+B)=(I(B)-D)/I(B*A+B):I(I+B)=I(I+B)-I(I-A+B)
160 IF ABS I(I-A+B)>H;$=""
170 NEXT B:IF $="" THEN 30
180 FOR B=1 TO A:PRINT B;I(I+B):NEXT B:GOTO 20
190 J=...f1{J(I),K(I),...,I(I+A)}...
200 K=...f2{J(I),K(I),...,I(I+A)}...
... I(A)=...fn{J(I),K(I),...,I(I+A)}...:RETURN

Размер: 566+, I(n*(n+2))

Пример:

F1=Цx1-x2
F2=x1+exp(-x2)-1
x(0)=(1, 1);

e=10-6

В программе:

A=2
I=A*A+A
190 J=SQR P-Q:K=P+EXP -Q-1:RETURN

Ответ: x1*=0,5105908269; x2*=0,7145563847.

5.3. Поиск экстремальных точек

Поиск экстремальных точек функции - операция, применяемая при оптимизации решений. Соответствующие прикладные программы получили чрезвычайно широкое распространение. Ресурсы "МК 85" не позволяют реализовывать оптимизационные методы в приемлемом для использования виде. Исключение составляют лишь некоторые алгоритмы непрерывной безусловной оптимизации.

Одномерная задача представлена двумя высокоэффективными методами, первый из которых основан на процедуре деления отрезка, второй - на замене функции интерполяционным многочленом.

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

5.3.1 Одномерная задача

Программа 82. Метод золотого сечения. Процесс сходится к локальному минимуму. Для определения максимума следует изменить знак алгебраического отношения в строке 30. В качестве исходных данных указывается интервал поиска [a,b]. Деление отрезка в пропорции золотого сечения позволяет обойтись минимальным количеством вычислений F(x). Требуемое число итераций:

ln((b-a)/e) / ln (1/d), где d=(Ц5-1)/2 - пропорция золотого сечения.

Положить: xb=a+d(b-a); xa=b-d(b-a);
Пока xb-xa> e, повторять:
если F(xb)>F(xa), тогда b=xb; xb=xa; xa=b-d(b-a);
иначе a=xa; xa=xb; xb=a+d(b-a).

10 INPUT "a",A,"b",B:G=SQR 5/2-.5:H=0:GOSUB 50:H=1:GOSUB 50
20 IF C=D;PRINT "x";C,"F";E:GOTO 10
30 H=0:IF E>F;H=1
40 A(H)=D(-H):D(-H)=C(H):F(-H)=E(H):GOSUB 50:GOTO 20
50 C(H)=A(H)+G*(B(-H)-A(H)):E(H)=...f{C(H)}...:RETURN

Размер: 162+, H

Пример: F(x)=|x-1|; a=0; b=2.

В программе: 50 ... :E(H)=ABS (C(H)-1):RETURN

Ответ: x=1; F(1)=0.

Программа 83. Метод квадратичной интерполяции:

xi+1={F(xi-s)(2xi+s)+4F(xi)xi+F(xi+s)(2xi-s)}/{F(xi-s)-2F(xi)+F(xi+s)}/2;

где s - шаг опорных точек.

Если F(x) - гладкая, процесс сходится; при F"(x0)>0 - к минимуму, при F"(x0)<0 - к максимуму.

10 INPUT "s",A,"x",B
20 E=0:C=B:GOSUB 50:E=C:C=B-A:GOSUB 50:D=C:C=B+A:GOSUB 50
30 C=(D-C)/(D+C)/2:B=B+C*A:IF ABS C>1E-9 THEN 20
40 PRINT "x";B,"F";E:GOTO 10
50 C=...f{C}-E...:RETURN

Размер: 123+, E

Пример: F(x)=exp(-x2); s=1; x0=1/2.

В программе: 50 C=EXP(-C*C)-E:RETURN

Ответ: x=0; F(0)=1.

5.3.2 Многомерная задача

Программа 84. Метод покоординатного спуска. Требует задания начального приближения r0, начального шага поиска s и допустимой линейной погрешности e. По каждой из n координат проводится грубый одномерный поиск минимума методом поразрядного приближения с шагом s и s=-s/3, затем шаг уменьшается в той же пропорции и процедура повторяется. Условие завершения - s<e.

10 A=n:FOR D=1 TO A:PRINT D;:INPUT G(D):NEXT D:GOSUB 80
20 INPUT "s",C,"e",E
30 FOR D=1 TO A:B=C
40 G(D)=G(D)+C:G=F:GOSUB 80:IF F<G THEN 40
50 IF B=C;C=-C/3:GOTO 40
60 NEXT D:IF ABS C>E THEN 30
70 FOR D=1 TO A:PRINT D;G(D):NEXT D:PRINT "F";F:GOTO 10
80 F=...f{H,I,...,G(A)}...:RETURN

Размер: 159+, G(n)

Пример: F(x,y,z)= ex+y+z/x/y2/z3
x0=0,5; y0=0,5; z0=0,5; s=0,5; e= 10-4.

В программе:

10 A=3: ...
80 F=EXP (H+I+J)/H/I/I/J-3:RETURN

Ответ: x*=0,9997715; y*=1,99961875; z*=3,00071107;
Fmin=3,73545085.

Программа 85. Метод наискорейшего спуска. Требует задания начального приближения r0, начального шага поиска s и допустимой линейной погрешности e. В текущей точке ri численно определяется направление, противоположное местному градиенту минимизируемой функции. Вдоль этого направления проводится точный одномерный поиск экстремума поразрядным приближением с кратностью шага 3. Найденная точка принимается за новое приближение ri+1, для которого процедура повторяется с начальным шагом s=|ri+1 - ri|. Условие окончания счета - s<e.

10 A=n:FOR B=1 TO A:PRINT B;:INPUT H(B):NEXT B
20 INPUT "s",G,"e",E
30 IF G<E;FOR B=1 TO A:PRINT B;H(B):NEXT B:PRINT "F";F:GOTO 20
40 C=G:H=0:GOSUB 90:D=F:FOR B=1 TO A:H(B)=H(B)-E:GOSUB 90
50 H(B+A)=(F-D)/E:H(B)=H(B)+E:H=H+H(B+A)-2:NEXT B:H=SQR H:G=0
60 FOR B=1 TO A:H(B)=H(B)+C/H*H(B+A):NEXT B:G=G+C
70 GOSUB 90:IF F≥D;C=-C/3:IF ABS C<E/5 THEN 30
80 D=F:GOTO 60
90 F=...f{I,J,...,H(A)}...:RETURN

Пример: 258+, H(2*n)

Пример: F(x,y)=x2+y2; x0=1; y0=1; s=1,5; e =10-5.

В программе:

10 A=2: ...
90 F=I*I+J*J:RETURN

Ответ: x*=0; y*=0; Fmin=0.

5.4. Численное интегрирование

Методы численного интегрирования имеют дело с определенным интегралом вида:

I=тabf(x)dx.

Несобственные интегралы заменяются определенным интегралом путем выбора соответствующих пределов a, b и деления интервала интегрирования на части. Возможны следующие варианты:

1. Если функция содержит разрыв второго рода (f(x)=г), то выделяется e-окрестность разрыва, в которой интеграл вычисляется многократно, с поэтапным уменьшением шага интегрирования, до тех пор, пока результат не сойдется (функция должна быть интегрируемой). Затем с обычным шагом вычисляются интегралы вокруг окрестности разрыва. Суммарный интеграл:

I=тax-ef(x)dx+тx-ex+ef(x)dx +тx+ebf(x)dx.

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

I=т0гf(x)dx+т0bf(x)dx +тbb+sf(x)dx+тb+sb+2sf(x)dx+...

где s - шаг интегрирования, не сойдется. f(x) должна быть абсолютно интегрируемой - быстро затухать при стремлении x к бесконечности. Шаг интегрирования увеличивают по ходу вычислений, используя правило: если (1-Ii-1/Ii)<e, тогда s=s+s. Здесь e - порог, задаваемый в зависимости от требуемой точности.

5.4.1 Квадратурные формулы Ньютона - Котеса

Интегрирование ведется по равноотстоящим отсчетам. Представлены формулы 2-го и 6-го порядков. Интеграторы других порядков могут быть получены путем незначительной модификации программы 87.

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

Программа 86. Формула Симпсона.

I=тx0x2f(x)dx = D x(4f(x1)+f(x0)+f(x2))/3

Формула точна для многочлена третьей степени.

10 INPUT "a",A,"b",B,"n",C:C=(B-A)/C/2:D=0:GOSUB 40:E=F
20 IF A≤B;GOSUB 40:D=2-D:E=E+(2+D)*F:GOTO 20
30 PRINT (E-F)*C/3:GOTO 10
40 F=...f{A}...:A=A+C:RETURN

Размер: 109+, F

Пример: f(x)=x3; a=0; b=1; n=2.

В программе: 40 F=A-3: ...

Ответ: I=0,25.

Программа 87. Формула Ньютона - Котеса.

I=тx0x6f(x)dx
      = D x{272f(x3)+27(f(x2)+f(x4))
      +216(f(x1)+f(x5)) +41(f(x0)+f(x6))}/140

Формула точна для многочлена седьмой степени и применяется к громоздким зависимостям.

10 INPUT "a",A,"b",B,"n",C:G=272:H=27:I=216:J=82
20 C=(B-A)/C/6:GOSUB 50:E=41*F
30 IF A>B;PRINT (E-41*F)*C/140:GOTO 10
40 FOR D=-2 TO 3:GOSUB 50:E=E+G(ABS D)*F:NEXT D:GOTO 30
50 F=...f{A}...:A=A+C:RETURN

Размер: 142+, J

Пример: f(x)=x7; a=0; b=2; n=2.

В программе: 50 F=A-7: ...

Ответ: I=32.

5.4.2 Квадратурные формулы Чебышева и Гаусса

Интегрирование ведется по неравноотстоящим отсчетам. Общая формула:

т-11f(t)dt = е(j)Ajf(tj), j=1...k;

где k - число отсчетов и порядок интегратора;
Aj - удельный вес отсчета;
f(tj) - отсчет интегрируемой функции f(t) в точке -1 ≤ tj ≤ 1.

В формулах Чебышева все Aj=2/k.

Программа 88. Формула Чебышева 3-го порядка.

t1=-1/Ц2; t2=0; t3=1/Ц2.

10 G=SQR 2/4:INPUT "a",A,"b",B,"n",C:C=(B-A)/C:F=0
20 FOR A=A+C/2 TO B STEP C:FOR D=-1 TO 1:E=A+SGN D*F(ABS D)*C
30 E=...f{E}...
40 F=F+E:NEXT D:NEXT A:PRINT F*C/3:GOTO 10

Размер: 105+, G

Пример: f(x)=x3; a=0; b=1; n=2.

В программе: 30 F=F+E-3:NEXT D:NEXT A:PRINT F*C/3:GOTO 10

Ответ: I=0,25.

Программа 89. Формула Чебышева 4-го порядка.

t1=-0,7946544723; t2=-0,1875924741.
t4=0,7946544723; t3=0,1875924741.

10 INPUT "a",A,"b",B,"n",C:C=(B-A)/C:F=0
20 FOR A=A+C/2 TO B STEP C:D=C
30 E=.0937962:GOSUB 50:E=.3972724:GOSUB 50:IF D=C;D=-D:GOTO 30
40 NEXT A:PRINT F*C/4:GOTO 10
50 E=A+D*E:F=F+...f{E}...:RETURN

Размер: 134+, H

Пример: f(x)=x4; a=0; b=1; n=2.

В программе: 50 E=A+D*E:F=F+E-4:RETURN

Ответ: I=0,19999.

Программа 90. Формула Гаусса 3-го порядка.

t1=-Ц(3/5); t2=0; t3=Ц(3/5); A1,3=5/9; A2=8/9.

10 G=SQR .15:H=8:I=5:INPUT "a",A,"b",B,"n",C
20 C=(B-A)/C:F=0:FOR A=A+C/2 TO B STEP C
30 FOR D=-1 TO 1:E=A+SGN D*F(ABS D)*C
40 E=...f{E}...
50 F=F+H(ABS D)*E:NEXT D:NEXT A:PRINT F*C/18:GOTO 10

Размер: 128+, I

Пример: f(x)=x5; a=0; b=1; n=2.

В программе: 40 E=E-5: ...

Ответ: I=0,1666666666.

Программа 91. Сложная квадратура Гаусса для двумерного интеграла. Вычисляется интеграл по прямоугольной области:

I=тa1b1тa2b2 f(x1, x2)dx1dx2.

Интервал [a1,b1] делится на n1 отрезков, [a2,b2] - на n2 отрезков. В результате получается n1n2 элементарных прямоугольников, к каждому из которых применяется квадратура, включающая четыре точки с координатами (¦Ц(1/3), ¦Ц(1/3)) и равными весами.

10 INPUT "a1",A,"b1",B,"n1",C,"a2",D,"b2",E,"n2",F:G=.5/SQR 3
20 H=0:C=(B-A)/C:F=(E-D)/F:FOR A=A+C/2 TO B STEP C
30 FOR L=D+F/2 TO E STEP F:FOR I=-1.5 TO 1.5
40 J=A+SGN I*G*C:K=L+SGN (ABS I-1)*G*F
50 J=...f{J,K}...
60 H=H+J:NEXT I:NEXT L:NEXT A:PRINT H*C*F/4:GOTO 10

Размер: 186+, L

Пример: f(x)=x1+x2-x1x2; a1=0; b1=2; n1=4; a2=0; b2=3; n2=5.

В программе:

50 H=H+J*J+K*K-K*J
60 NEXT I:NEXT L:NEXT A:PRINT H*C*F/4:GOTO 10

Ответ: I=17.


Программы для МК-85. Метод наименьших квадратов

Метод наименьших квадратов

6.4. Метод наименьших квадратов

Аппроксимация по методу наименьших квадратов (МНК) доставляет минимум функционалу

е(i)(yi-f(xi))2; i=1...n

где (xi,yi) - исходные данные,
f(x) - аппроксимирующая функция.

Если число определяемых параметров функции равно числу отсчетов n, аппроксимация вырождается в интерполяцию.

Линейная регрессия:

y=ax+b.

Расчетные формулы:

d=е(i)xi;
      e=е(i)yi;
      f=е(i)xi2;
      g=е(i)yi2;
      h=е(i)xiyi;
      i=1...n
a=(de-nh)/(d2-nf); b=(e-ad)/n;
      s=Ц((g-be-ah)/n)

a, b - параметры аппроксимирующей прямой;
s - среднеквадратическая ошибка.

Программа 103. Линейная регрессия

10 VAC
20 PRINT C+1;:INPUT "x",A,"y",B
30 C=C+1:D=D+A:E=E+B:F=F+A*A:G=G+B*B:H=H+A*B:GOTO 20
40A=(D*E-C*H)/(D*D-C*F):B=(E-A*D)/C
50 PRINT "a";A,"b";B,"s";SQR ((G-B*E-A*H)/C):GOTO 20

Размер: 149, H

Пример:

xi | 1 | 2 | 3 | 4 | 5
yi | 1 | 2.02 | 3 | 4.1 | 4.96

Завершение ввода данных: [AC] RUN 40[EXE].

Ответ: a=1; b=0,016; s=0,0463.

Программа 104. Регрессия, сводимая к линейной. Трансцендентная двухпараметрическая зависимость y=f(a,b,x) путем замены переменных сводится к линейной

fy(x,y)=fa(a,b)fx(x,y)+fb(a,b).

Рассмотрим частный случай y=px/(q+x).

Линеаризация: 1/y=(q/p)(1/x)+1/p;

откуда p=1/b, q=a/b, y=1/y*, x=1/x*.

10 VAC
20 PRINT C+1;:INPUT "x",A,"y",B:A=1/A:B=1/B
30 C=C+1:D=D+A:E=E+B:F=F+A*A:G=G+B*B:H=H+A*B:GOTO 20
40A=(D*E-C*H)/(D*D-C*F):B=(E-A*D)/C
50 PRINT "p";1/B,"q";A/B,"s";SQR ((G-B*E-A*H)/C):GOTO 20

Размер: 165, H

Пример:

xi | 1 | 2 | 3 | 4 | 5
yi | 1 | 1.5 | 1.8 | 2 | 2.15

Завершение ввода данных: [AC] RUN 40[EXE].

Ответ: p=3,007...; q=2,007...; s=5,4 10-4.

Вычисленная оценка погрешности s относится к линеаризованной зависимости.

Программа 105. Нелинейная регрессия. Представленный здесь упрощенный алгоритм МНК принимает все отсчеты равноточными, кроме того, не оценивает погрешность.

В качестве примера рассмотрим зависимость:

y=ax/(1-bx)+cx+d.

Эквивалентное уравнение, линейное по параметрам:

y=bxy+(-bc)x2+(a+c-db)x+d;

в общем виде: y=a1y+a2x2+a3x+a4;
откуда: b=a1, d=a4, c=-a2/b, a=a3-c-db.

Уравнения для частных производных функционала МНК по неизвестным параметрам

T/T aj е(i)(y-a1xy-a2x2-a3x-a4)2=0, i=1...4

образуют систему:

Sa=r

Формулы вычисления элементов матриц:

S |  |  |  | r
е x2y2 | е x3y | е x2y | е xy | е xy2
е x3y | е x4 | е x3 | е x2 | е x2y
е x2y | е x3 | е x2 | е x | е xy
е xy | е x2 | е x | n | е y

Для решения системы использован компактный алгоритм метода Гаусса для симметричных матриц.

Размещение в памяти:

I J K L E
M N O P F
Q R S T G
U V W X H

В программной области P0:

10 VAC
20 INPUT "x",C,"y",D:A=C*C:B=D*D:E=E+C*B:F=F+A*D:G=G+C*D:H=H+D
30 I=I+A*B:M=M+A*C*D:N=N+A*A:R=R+A*C:S=S+A:W=W+C:X=X+1:GOTO 20

В программной области P1:

10 Q=F:U=G:V=S:A=4
20 FOR D=2 TO A:D(A+D)=E(A*D)/E(A)
30 FOR B=2 TO A:IF D>B;D(A*B+D)=D(A*D+B)/D(A*B+B):GOTO 50
40 FOR C=1 TO D-1:D(A*B+D)=D(A*B+D)-D(A*B+C)*D(A*C+D):NEXT C
50 NEXT B:NEXT D:E=E/E(A):FOR B=2 TO A:FOR C=1 TO B-1
60 D(B)=D(B)-D(A*B+C)*D(C):NEXT C:D(B)=D(B)/D(A*B+B):NEXT B
70 FOR B=A-1 TO 1 STEP -1:FOR C=B+1 TO A
80 D(B)=D(B)-D(A*B+C)*D(C):NEXT C:NEXT B
90 J=E:L=H:K=-F/J:I=G-K+J*L
100 FOR B=1 TO A:PRINT CHR (64+B);H(B):NEXT B

Пример:

xi | 0.5 | 0.8 | 1 | 1.5
yi | 14 | 17.2 | 20 | 34

Ввод данных: [P0] .5[EXE] 14[EXE] .8[EXE] 17.2[EXE] 1[EXE] 20[EXE] 1.5[EXE] 34[EXE] [AC] [P1]

Ответ: a=3,0027...; b=0,4999...; c=3,9948; d=10,0009;
то есть y=3x/(1-0,5x)+4x+10.

6.5. Гармоническая аппроксимация

Программа 106. Вычисление коэффициентов Фурье

ak=2/m е(i)yicos(2p ik/m),
      i=0...m-1, k=0...na;
bk=2/m е(i)yisin(2p ik/m),
      i=0...m-1, k=0...nb.

Отсчеты yi аппроксимируемой функции делят ее период на m равных интервалов. Порядок аппроксимации n=na+nb.

10 VAC:MODE 5:INPUT "n",A,"m",B
20 FOR D=0 TO B-1:PRINT D;:INPUT E:FOR C=1 TO A
30 E(C)=E(C)+E*SIN (2*p*D/B*INT (C/2)+p*FRAC (C/2)):NEXT C
40 NEXT D:FOR C=1 TO A
50 PRINT CHR (66-2*FRAC (C/2));INT (C/2);2/B*E(C):NEXT C

Размер: 138, E(n)

Пример: n=2; m=3; y0=0; y1=0.8660254; y2=-0.8660254.

Ответ: a0=0; b1=1.

Программа 107. Гармоническая интерполяция. В отличие от предыдущей программы предполагается n=m. Интерполяционная формула:

f(j)=a0/2+е(k)bksin j k+е(k)akcos j k; k=1..nb,na

где ak, bk - коэффициенты Фурье (см. описание программы 106);

j=2p(x-x0)/T
      - фаза;
T - период;
x0 - левый край интервала-периода;
x - точка, в которой определяется значение функции.

10 VAC:MODE 5:INPUT "m",B
20 FOR D=0 TO B-1:PRINT D;:INPUT E:FOR C=1 TO B
30 F(C)=F(C)+E/B*SIN (2*p*D/B*INT (C/2)+p*FRAC (C/2))
40 NEXT C:NEXT D:INPUT "x0",A,"T",E
50 INPUT "x",D:D=(D-A)/E*2*p:F=G:FOR C=2 TO B
60 F=F+2*F(C)*SIN (D*INT (C/2)+p*FRAC (C/2))
70 NEXT C:PRINT F:GOTO 50

Размер: 184, F(m)

Пример: m=3; y0=0; y1=0.8660254; y2=-y1;
  x0=2p; T=2p;
  x=p/4.

Ответ: f(p/4)=0,7071067812.

Программа 108. Интерполяция спектра кусочно-линейной функции с разрывами. Значения действительной (af) и мнимой (bf) составляющих спектра на указанной частоте f определяются по формулам:

af=1/x/n{е(i)(yi-sin 2xi-y+i-1sin 2x(i-1) ) - sin x /x е(i)(yi--y+i-1sin x(2i-1) )};
bf=1/x/n{е(i)(y+i-1cos 2x(i-1) - yi-cos 2xi) - sin x /x е(i)(yi--y+i-1cos x(2i-1) )};
i=1...n

где x=p f D t, f - частота;

D t =ti-ti-1
  - длина линейного участка аппроксимируемой функции;
yi-=y(ti-0) --значение
  функции слева от разрыва или излома;
y+i=y(ti +0) --значение
  функции справа от разрыва или излома;
n -число интервалов разбиения периода функции.

10 INPUT "n",A:FOR B=1 TO A:PRINT
  B;
20 INPUT "t",J(3*B),"y-",I(3*B),"y+",K(3*B):NEXT B:K=K(3*A)
30 MODE 5:I=0:J=0:INPUT "f",C:C=2*p*C:FOR
  B=3 TO 3*A STEP 3
40 D=C*G(B):E=C*J(B):F=(I(B)-H(B))/(E-D):G=SIN D:D=COS D
50 H=SIN E:E=COS E:I=I+F*(E-D)+I(B)*H-H(B)*G
60 J=J+F*(H-G)-I(B)*E+H(B)*D:NEXT B:C=C/2*G(B)
70 PRINT "a";I/C,"b";J/C:GOTO 30

Размер: 262, K(3*n)

При вводе данных y+0 запрашивается в конце, как yn-.

Пример: n=2; t1=1 10-6; y1- =1; y+1=0; t2=4 10-6; y2- =0; y+2=1.

Задание 1: f=0,01. Ответ: af=0,5; bf=0.

Задание 2: f=2,5 105. Ответ: af=0,3183098862; bf=0,3183098862.

6.6. Преобразование Фурье

Прямое:

zk=е(i)uie-j2pik/n, i=0...n-1, k=0...n-1;

Обратное:

zi=е(k)ukej2pik/n, k=0...n-1, i=0...n-1.

Здесь zi=ui=xi+jyi
  - комплексные отсчеты функции времени;
zk=uk=xk+jyk
  - комплексные отсчеты функции частоты.

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

Программа 109. Программа допускает разбиение периода функции только на n=2m интервалов, где m - натуральное.

Ограничение размерности: m≤5.

10 INPUT "m",A,"П",B:C=2-A:DEFM 2*C-12:FOR D=1 TO C:PRINT D;
20 INPUT "x",M(D+D),"y",N(D+D):NEXT D:MODE 5
30 FOR D=1 TO A:E=2-(A-D):G=p/E:F=COS G:G=B*SIN G:H=1:I=0
40 FOR J=1 TO E:FOR K=J+J TO C+C STEP 4*E:L=K+E+E:M=M(K)-M(L)
50 N=N(K)-N(L):M(K)=M(K)+M(L):N(K)=N(K)+N(L)
60 M(L)=M*H-N*I:N(L)=N*H+M*I:NEXT K:K=H*F+I*G:I=I*F-H*G
70 H=K:NEXT J:NEXT D:FOR D=1 TO C-1:IF D≥E THEN 100
80 FOR H=0 TO 1:M(H)=M(E+E+H):M(E+E+H)=M(D+D+H)
90 M(D+D+H)=M(H):NEXT H
100 I=C/2
110 IF I<E;E=E-I:I=I/2:GOTO 110
120 E=E+I:NEXT D:IF B>0;D=1
130 FOR E=1 TO C:PRINT E,"x";M(E+E)/D,"y";N(E+E)/D:NEXT E

Размер: 438, N(2*2-m)

Для прямого преобразования необходимо задать "П"=1, для обратного - "П"=-1.

Пример: m=3; "П"=1;

xi | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0
yi | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0

Проверка по свойству обратимости. Применение обратного преобразования к результатам прямого даст исходную функцию с двойной вычислительной погрешностью. Достаточно подать команды B=-B и RUN 30.

Ниже представлена модификация программы 109, рассчитанная на число отсчетов n=64. Ввод и вывод данных ведется в режиме калькулятора.

10 FOR D=1 TO 6:E=64-2-D:G=p/E:F=COS G:G=SIN G:H=1:I=0
20 FOR J=1 TO E:FOR K=J+J TO 128 STEP 4*E:A=K+E+E:B=J(K)-J(A)
30 C=K(K)-K(A):J(K)=J(K)+J(A):K(K)=K(K)+K(A)
40 J(A)=B*H-C*I:K(A)=C*H+B*I:NEXT K:K=H*F+I*G:I=I*F-H*G
50 H=K:NEXT J:NEXT D:FOR D=1 TO 63:IF D≥E THEN 80
60 FOR H=0 TO 1:J(H)=J(E+E+H):J(E+E+H)=J(D+D+H)
70 J(D+D+H)=J(H):NEXT H
80 I=C/2
90 IF I<E;E=E-I:I=I/2:GOTO 90
100 E=E+I:NEXT D

Размер: 1212

Ввод/вывод данных

Прямое

40 ... :K=H*F+I*G:I=I*F-H*G
MODE 5
DEFM 113
L=x1
M=y1
N=x2
...
K(128)=y64
RUN
x1- L
y1- M
x2- N
...
y64- K(128)

Обратное

40 ... :K=H*F-I*G:I=I*F+H*G
MODE 5
DEFM 113
L=x1
M=y1
N=x2
...
K(128)=y64
RUN
x1- L/64
y1- M/64
x2- N/64
...
y64- K(128)/64


Программы для МК-85. Генераторы случайных чисел

8.2 Генераторы случайных чисел

Генераторы случайных чисел преобразуют стандартное равномерное распределение к требуемому виду, оставляя числа независимыми. В подразделе "Получение случайных чисел" представлены наиболее важные распределения.

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

В бейсике "МК 85" датчик псевдослучайных чисел вызывается функцией RAN#. Она возвращает числа, распределенные равномерно в интервале [0;1[. Рандомизация входа в последовательность выполняется автоматически и обеспечивается запоминанием последнего полученного числа. Выполнение вшитого теста (команда TEST) переводит датчик в исходное состояние.

8.2.1 Распределения дискретных величин

Программа 123. Распределение Бернулли

P(0)=1-p, P(1)=p, 0≤p≤1.

Размещение переменных. Здесь и далее результат выделен жирным шрифтом.

A | B
X | p

100 A=0:IF RAN#≤B;A=1
110 RETURN

Размер: 19, B

Программа 124. Биномиальное распределение

P(x)=n!/x!/(n-x)!px(1-p)x, x=0...n.

100 A=0:FOR C=1 TO D:IF RAN#≤B;A=A+1
110 NEXT C:RETURN

Размер: 31, D

Программа 125. Биномиальное распределение с малой вероятностью успеха

A | B | С | D
X | p |  | n

10 B=LN (1-B)
100 A=0:C=0
110 C=C+1+INT (LN RAN#/B):IF C>D;RETURN
120 A=A+1:GOTO 110

Размер: 57, D

Программа 126. Полиномиальное распределение

P(x1, ... , xm)= n!/x1!/.../xm!p1...pmxm, е(i)xi=n., е(i)pi=1, i=1...n.

A | B | G | F(A) | G(A) | F(A+A)
m | n | p1 | pm | X1 | Xm

10 F=0:FOR C=1 TO A:F(C)=F(C)+F:F=F(C):NEXT C
100 FOR C=1 TO A:F(C+A)=0:NEXT C
110 FOR E=1 TO B:D=RAN#:FOR C=1 TO A
120 IF D<F(C);F(C+A)=F(C+A)+1:GOTO 140
130 NEXT C
140 NEXT E:RETURN

Размер: 119, F(2*m)

Программа 127. Геометрическое распределение

P(x)=p(1-p)x-1, x=1,2, ...

A | B
X | p

100 A=0
110 A=A+1:IF RAN#>B THEN 110
120 RETURN

Размер: 27, B

Программа 128. Геометрическое распределение с малой вероятностью успеха

A | B
X | p

10 B=LN (1-B)
100 A=1+INT (LN RAN#/B):RETURN

Размер: 27, B

Программа 129. Гипергеометрическое распределение

P(m)=CmM Cn-mN-M/CnN

A | B | С | D
N | M | m | n

100 C=0:E=A:FOR F=1 TO D:IF RAN#≤B/E;C=C+1:B=B-1
110 E=E-1:NEXT F:RETURN

Размер: 49, F

Программа 130. Распределение Паскаля

P(x)= Cn-xn-1 px(1-p)n-x, x=1...n

B | С | D | E
p | X |  | n

10 B=LN (1-B)
100 C=0:FOR D=1 TO E:C=C+INT (LN RAN#/B):NEXT D:RETURN

Размер: 41, E

Программа 131. Распределение Пуассона

P(x)=ax/x!e-a, a>0, x=1,2, ...

A | B | С
a |  | X

10 A=EXP -A
100 B=1:C=0
110 B=B*RAN#:IF B≥A;C=C+1:GOTO 110
120 RETURN

Размер: 46, C

8.2.2 Распределения непрерывных величин

Программа 132. Равномерное распределение

j(x)=1/(b-a), a≤x<b.

A | B | С
a | b | X

Размер: 21, C

10 B=B-A
100 C=RAN#*B+A:RETURN

Гамма-распределение и его частные случаи

j(x)=lcxc-1e-lx/G(c), l>0, c>0, x>0.

Частный случай - c=1.

Программа 133. Экспоненциальное распределение

j(x)=le-lx

A | B
l | X

100 B=-LN RAN#/A:RETURN

Размер: 12, B

Частный случай - c натуральное

Программа 134. Распределение Эрланга

A | B | С
l | X | c

100 B=1:FOR D=1 TO C:B=B*RAN#:NEXT D:B=-LN B/A:RETURN

Размер: 32, D

В программах 135, 136 реализован метод браковки. Более быстрым является алгоритм, использующий полиномиальную аппроксимацию, однако он занял бы значительно больше места.

Программа 135. Гамма-распределение с параметром формы c<1

A | B | С | D | E | F | G
l |  | c | X |  |  |

10 B=1/C:C=1/(1-C):A=-2/A:F=0
100 F=F+.5:IF FRAC F=0;D=E:RETURN
110 G=RAN#:D=G-B:E=G-C:IF
  D+E>1 THEN 110
120 G=A*LN RAN#/(D+E):D=D*G:E=E*G:RETURN

Размер: 108, G

Программа 136. Гамма-распределение с произвольным параметром формы. Числа генерируются парами.

A | B | С | D | E | F | G | H | I | J
l |  | c |  | X2i | X2i+1 |  |  |  |

10 D=INT C+2:B=C/D:C=1/(1-B):B=1/B:A=-2/A
100 E=0:F=0:FOR G=1 TO D
110 H=RAN#:I=H-B:J=H-C:IF
  I+J>1 THEN 110
120 H=A*LN RAN#/(I+J):E=E+I*H:F=F+J*H:NEXT G:RETURN

Размер: 122, J

Программа 137. Бета-распределение

j(x)=xv-1(1-x)w-1G(v+w)/G(v)/G(w), 0≤x≤1, v>0, w>0

Метод браковки для бета-распределения основан на получении чисел, распределенных по гамма-, с использованием соотношения:

b(v, w)=g(1, v)/(g(1, v)+g(1, w))

Программа записана для частного случая, когда v и w - целые.

A | B | С | D | E | F
v | w |  |  |  | X

100 C=A:GOSUB 110:F=-LN D:C=B:GOSUB 110:F=F/(F-LN D):RETURN
110 D=1:FOR E=1 TO C:D=D*RAN#:NEXT E:RETURN

Размер: 63, F

Программа 138. Распределение Вейбулла

j(x)=cxc-1/bcexp(-(x/b)c), b>0, c>0, x≥0.

100 A=B*(-LN RAN#)-(1/B):RETURN

Размер: 20, C

Нормальное распределение

j(x)=1/s/Ц(2p) exp{-1/2 ((x-m)/s)2}, s>0.

Ниже приведены генераторы, основанные на разных принципах. Какой из них быстрее, определяется особенностями ЭВМ.

Программа 139. Точный генератор. Создает числа попарно:

X2i=Ц(-2ln R2i) sin(2p R2i+1)s+m;
X2i+1=Ц(-2ln R2i) cos(2p R2i+1)s+m;

где R2i- элемент равномерной последовательности;

X2i- элемент нормальной последовательности.

A | B | С | D | E
m | s | X2i | X2i+1 |

10 MODE 5
100 D=SQR (-2*LN RAN#):E=2*p*RAN#
110 C=D*SIN E*B+A:D=D*COS E*B+A:RETURN

Размер: 47, E

Программа 140. Использование закона больших чисел

X=(е(i)Ri-6)s+m, i=1...12

A | B | С | D
m | s | X |

Размер: 35, D

100 C=-6:FOR D=1 TO 12:C=C+RAN#:NEXT D:C=C*B+A:RETURN

Программа 141. Использование полиномиальной аппроксимации

X={Y-41/13440/n2(Y5-10Y3+15Y)}s+m,

где Yi=е(j)Rni+j, j=1...n, n=2

A | B | С
m | s | X

100 C=(RAN#-RAN#)*1.224744871
110 C=(C-7.626488E-4*C*(15+C*C*(C*C-10)))*B+A:RETURN

Размер: 67, C

Программа 142. Распределение Релея

j(x)=x/a exp(-x2/2/a), a>0, x>0.

A | B
a | X

Размер: 17, B

100 B=A*SQR (-2*LN RAN#):RETURN

Программа 143. Логистическое распределение

j(x)= p/s/Ц3 exp((x-m)p/s/Ц3)/{1+exp((x-m)p/s/Ц3)}2

A | B | С
m | s | X

10 B=SQR 3/p*B
100 C=RAN#:C=LN (C/(1-C))*B+A:RETURN

Размер: 36, C


Программы для МК-85. Оглавление

150 ВЫЧИСЛИТЕЛЬНЫХ ПРОГРАММ
для карманного компьютера
"Электроника МК 85"

Книга содержит тексты, краткие описания и примеры использования 150 бейсик-программ для "Электроники МК 85". Программы охватывают основные разделы вычислительной математики. Материал может быть использован как в практических, так и в учебных целях.

Для студентов, инженерно-технических и научных работников - обладателей "Электроники МК 85".

Издание НПП "Вектор", Москва, 1992

Оглавление

  • 0. Предварительные замечания
    • 0.1. О пользователе
      • 0.1.1 Кому адресована библиотека
      • 0.1.2 Что нужно иметь и уметь
      • 0.1.3 Где искать совместимые программы
    • 0.2. О библиотеке
      • 0.2.1 Классификация программ
      • 0.2.2 Стиль программирования
      • 0.2.3 Стиль описания программ
    • 0.3. О творческом использовании библиотеки
      • 0.3.1 Внесение изменений в программу
      • 0.3.2 Средства обнаружения ошибок
      • 0.3.3 Средства диалога
  • 1. Числа
    • 1.0. Введение
      • 1.0.1 Действительные числа
      • 1.0.2 Логические переменные
      • 1.0.3 Бинарные числа
      • 1.0.4 Комплексные числа
    • 1.1. Числа с произвольным основанием
      • Программа 1. Перевод произвольного числа в десятеричное
      • Программа 2. Перевод десятеричного числа в число с требуемым основанием
      • Программа 3. Перевод десятеричного числа в число с основанием от 2 до 16
    • 1.2. Делимость натуральных чисел
      • Программа 4. Наибольший общий делитель двух чисел
      • Программа 5. Разложение числа на простые сомножители
      • Программа 6. Наименьшее общее кратное двух чисел
      • Программа 7. Наибольший общий делитель нескольких чисел
      • Программа 8. Наименьшее общее кратное нескольких чисел
    • 1.3. Комбинаторика
      • Программа 9. Число перестановок
      • Программа 10. Число размещений
      • Программа 11. Число сочетаний
    • 1.4. Числовые последовательности
      • Программа 12. Последовательность простых чисел
      • Программа 13. Последовательность чисел Бернулли
      • Программа 14. Последовательность s-чисел Фибоначчи
    • 1.5. Сортировка чисел
      • Программа 15. Метод пузырька
  • 2. Матрицы
    • 2.0. Введение
    • 2.1. Матричное преобразование
      • Программа 16. Матрица вращения в трехмерном пространстве
    • 2.2. Произведение матриц
      • Программа 17
    • 2.3. Детерминант. Обратная матрица
      • Программа 18. Модифицированный метод Гаусса - Жордана
    • 2.4. Характеристический многочлен
      • Программа 19. Метод Данилевского
    • 2.5. Собственные значения
      • Программа 20. Метод вращений с преградами
    • 2.6. Системы линейных уравнений
      • Программа 21. Метод вращений
      • Программа 151. Метод Гаусса с выбором максимального элемента.
      • Программа 22. Метод Зейделя для разреженных матриц
  • 3. Степенные многочлены
    • 3.0. Введение
    • 3.1. Вычисление значений
      • Программа 23. Случай действительного аргумента
      • Программа 24. Случай комплексного аргумента
    • 3.2. Оценка устойчивости
      • Программа 25. Критерий Рауса - Гурвица
    • 3.3. Арифметические операции над многочленами
      • Программа 26. Перемножение многочленов
      • Программа 27. Деление многочленов
      • Программа 28. Общий множитель двух многочленов
    • 3.4. Интегрирование и дифференцирование
      • Программа 29. Интегрирование и дифференцирование
    • 3.5. Замена переменной
      • Программа 30. Линейная замена переменной
    • 3.6. Разложение на сомножители
      • Без программы. Квадратное уравнение
      • Программа 31. Кубическое уравнение
      • Программа 32. Метод Хичкока
  • 4. Степенные ряды
    • 4.0. Введение
    • 4.1. Вычисление систем ортогональных многочленов
      • Программа 33. Многочлены Чебышева
      • Программа 34. Многочлены Эрмита
      • Программа 35. Многочлены Лежандра
      • Программа 36. Многочлены Лагерра
    • 4.2. Вычисление значений ортогональных многочленов
      • Без программы. Многочлены Чебышева
      • Программа 37. Многочлены Лежандра
      • Программа 38. Многочлены Лагерра
      • Программа 39. Многочлены Эрмита
    • 4.3. Вычисление функций с помощью рядов
      • Программа 40. Интеграл вероятностей. Вычисление коэффициентов ряда
      • Программа 41. Интеграл вероятностей. Вычисление значений
      • Программа 42. Интегральный синус
      • Программа 43. Интегральный косинус
      • Программа 44. Совместное вычисление Ci(x) и Si(x)
      • Интегралы Френеля
        • Программа 45. Интегралы Френеля. Случай x≤3
        • Программа 46. Интегралы Френеля. Случай x>3
      • Программа 47. Функция Лапласа
      • Программа 48. Дилогарифм
  • 5. Функции в форме выражений
    • 5.0. Введение
      • 5.0.1 Понятие функции
      • 5.0.2 Встроенные функции бейсика "Электроники МК 85"
      • 5.0.3 Запись функций-выражений на бейсике
      • 5.0.4 Использование циклов и подпрограмм
    • 5.1. Вычисление значений
      • 5.1.1 Интегральные показательные функции
        • Программа 49. Функция En(x)
        • Программа 50. Функция Ei(x)
        • Программа 51. Функция an(x)
        • Программа 52. Функция bn(x)
      • 5.1.2 Факториальная функция и ее производные
        • Программа 53. Гамма-функция. Аппроксимация многочленом
        • Программа 54. Гамма-функция. Аппроксимация формулой Стирлинга
        • Программа 55. Дигамма-функция
        • Программа 56. Тригамма-функция
        • Программа 57. Неполная Гамма-функция
      • 5.1.3 Функции Бесселя
        • Программа 58. Функции Бесселя дробного порядка
        • Программа 59. Функции Бесселя натурального порядка
        • Программа 60. Модифицированные функции Бесселя дробного порядка
        • Программа 61. Модифицированные функции Бесселя натурального порядка
        • Программа 62. Сферические функции Бесселя
      • 5.1.4 Функции Кельвина
        • Программа 63. Функции Кельвина ber x и bei x
        • Программа 64. Функции Кельвина bernx и beinx, kernx и keinx
        • Гипергеометрические функции
          • Программа 65. Гипергеометрическая функция
          • Программа 66. Вырожденная гипергеометрическая функция
      • 5.1.5 Эллиптические функции
        • Программа 67. Полные эллиптические интегралы
        • Программа 68. Эллиптические функции Якоби sn(u,m), cn(u,m) и dn(u,m)
        • Программа 69. Тэта-функция Якоби
        • Программа 70. Дзэта-функция Якоби
        • Программа 71. Неполные эллиптические интегралы
    • 5.2. Поиск нулей
      • 5.2.1 Прямые методы
        • Программа 72. Метод простых итераций
        • Программа 73. Метод Эйткена - Стеффенсона
      • 5.2.2 Методы деления отрезка
        • Программа 74. Метод дихотомии
      • 5.2.3 Метод Ньютона и его модификации
        • Программа 75. Метод Ньютона
        • Программа 76. Модифицированный метод Ньютона
        • Программа 77. Метод Рыбакова
      • 5.2.4 Интерполяционные методы
        • Программа 78. Комбинированный метод секущих - хорд
        • Программа 79. Метод обратной квадратичной интерполяции
      • 5.2.5 Методы решения систем нелинейных уравнений
        • Программа 80. Метод Зейделя
        • Программа 81. Метод Ньютона - Канторовича
    • 5.3. Поиск экстремальных точек
      • 5.3.1 Одномерная задача
        • Программа 82. Метод золотого сечения
        • Программа 83. Метод квадратичной интерполяции
      • 5.3.2 Многомерная задача
        • Программа 84. Метод покоординатного спуска
        • Программа 85. Метод наискорейшего спуска
    • 5.4. Численное интегрирование
      • 5.4.1 Квадратурные формулы Ньютона - Котеса
        • Программа 86. Формула Симпсона
        • Программа 87. Формула Ньютона - Котеса
      • 5.4.2 Квадратурные формулы Чебышева и Гаусса
        • Программа 88. Формула Чебышева 3-го порядка
        • Программа 89. Формула Чебышева 4-го порядка
        • Программа 90. Формула Гаусса 3-го порядка
        • Программа 91. Сложная квадратура Гаусса для двумерного интеграла
  • 6. Функции в форме таблиц
    • 6.0. Введение
    • 6.1. Численное дифференцирование
      • Программа 92. Вычисление трех производных по семи точкам
      • Программа 93. Вычисление первой производной с интерполяцией по пяти точкам
    • 6.2. Интерполяция степенными многочленами
      • 6.2.1 Вычисление значений интерполяционного многочлена
        • Программа 94. Общий случай
        • Программа 95. Случай регулярной решетки
      • 6.2.2 Вычисление коэффициентов интерполяционного многочлена
        • Программа 96. Общий случай
        • Программа 97. Случай регулярной решетки
      • 6.2.3 Метод Чебышева
        • Программа 98. Аппроксимация методом Чебышева
        • Программа 99. Вычисление коэффициентов аппроксимирующего многочлена
    • 6.3. Сплайн-интерполяция
      • 6.3.1 Локальный кубический сплайн
        • Программа 100. Интерполяция локальным сплайном
      • 6.3.2 Глобальный кубический сплайн
        • Программа 101. Вычисление вторых производных
        • Программа 102. Интерполяция глобальным сплайном
    • 6.4. Метод наименьших квадратов
      • Программа 103. Линейная регрессия
      • Программа 104. Регрессия, сводимая к линейной
      • Программа 105. Нелинейная регрессия
    • 6.5. Гармоническая аппроксимация
      • Программа 106. Вычисление коэффициентов Фурье
      • Программа 107. Гармоническая интерполяция
      • Программа 108. Интерполяция спектра кусочно-линейной функции с разрывами
    • 6.6. Преобразование Фурье
      • Программа 109. Преобразование Фурье.
  • 7. Дифференциальные уравнения
    • 7.0 Введение
    • 7.1. Уравнение первого порядка
      • Программа 110. Метод Рунге - Кутта 4-го порядка
      • Программа 111. Метод Милна 6-го порядка
      • Программа 112. Метод Хемминга 4-го порядка
    • 7.2. Уравнение второго порядка
      • Программа 113. Метод Рунге - Кутта 4-го порядка
    • 7.3. Уравнение n-го порядка
      • Программа 114. Метод Рунге - Кутта 4-го порядка
    • 7.4. Система обыкновенных дифференциальных уравнений
      • Программа 115. Метод Рунге - Кутта 4-го порядка
      • Программа 116. Метод Рунге - Кутта - Мерсона
  • 8. Случайные числа
    • 8.1. Обработка выборочных данных
      • 8.1.1 Гистограмма. Программа 117
      • 8.1.2 Моменты. Программа 118
      • 8.1.3 Корреляция. Программа 119
      • 8.1.4 Тесты
        • Программа 120. Критерий согласия c2
        • Программа 121. Критерий Стьюдента
        • Программа 122. Критерий Фишера
    • 8.2. Получение случайных чисел
      • 8.2.1 Распределения дискретных величин
        • Программа 123. Распределение Бернулли
        • Программы 124, 125. Биномиальное распределение
        • Программа 126. Полиномиальное распределение
        • Программы 127, 128. Геометрическое распределение
        • Программа 129. Гипергеометрическое распределение
        • Программа 130. Распределение Паскаля
        • Программа 131. Распределение Пуассона
      • 8.2.2 Распределения непрерывных величин
        • Программа 132. Равномерное распределение
        • Гамма-распределение и его частные случаи
          • Программа 133. Экспоненциальное распределение
          • Программа 134. Распределение Эрланга
          • Программы 135, 136. Гамма-распределение
        • Программа 137. Бета-распределение
        • Программа 138. Распределение Вейбулла
        • Нормальное распределение
          • Программа 139. Точный генератор
          • Программа 140. Использование закона больших чисел
          • Программа 141. Использование полиномиальной аппроксимации
        • Программа 142. Распределение Релея
        • Программа 143. Логистическое распределение
  • 9. Графы и сети
    • 9.0. Введение
    • 9.1. Связывающая сеть.
      • Программа 144. Алгоритм Прима
    • 9.2. Геодезическая сеть.
      • Программа 145. Уравновешивание невязок углов
    • 9.3. Сеть систем массового обслуживания. Программа 146. Оценка средних времен ожидания
    • 9.4. Граф вероятных состояний.
      • Программа 147. Система уравнений Колмогорова
    • 9.5. Герт-сеть.
      • Программа 148. Правило Мейсона
    • 9.6. Перт-сеть.
      • Программа 149. Оценка раннего времени выполнения
    • 9.7. Дерево решений.
      • Программа 150. Задача об инвестициях


Created using
htm2chm 3.0

Download it NOW!