close

Вход

Забыли?

вход по аккаунту

?

FreeBASIC23

код для вставкиСкачать
Пособие для старшеклассников, младшекурсников и преподавателей по вычислительным методам на языке FreeBASIC.
1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из множества фрагментов. Этот фрагмент, надеюсь, тоже будет
интересен для старшеклассников, младшекурсников и их преподавателей.
Фрагмент 23. "Лапидарий: Долг платежом красен"
"Долг платежом красен",
"Всякие займы платежом красны",
"Долг платежом красен, а займы отдачею",
"Бери да помни! Не штука занять, штука отдать".
Даль В.И. Сборник "Пословицы русского народа" [01]
Пословицы русского народа. Сборник Владимира Даля, 2 тома.
CПб.-М.: Издание книгопродавца-типографа М.О.Вольфа. 1879. 750+642 с.
можно скачать, например, по ссылке:
http://bukvy.net/books/hobby/90312-Poslovitsi-russkogo-naroda.html
Размер архива: Dal.rar ~13,32 Мб, 2 файла формата DjVu
2
или почитать в режиме online на сайте:
http://hobbitaniya.ru/
"Хранители сказок" где собраны сказки для детей и взрослых.
http://hobbitaniya.ru/dal/
Даль Владимир Иванович (1801 — 1872). Пословицы русского народа.
http://hobbitaniya.ru/dal/dal99.php
Даль В.И. Пословицы русского народа > Займы.
Сборник "Пословицы русского народа", включающий более 30 тысяч пословиц,
поговорок и прибауток и впервые опубликованный в 1861-1862 годах, является одним
из основных трудов великого русского этнографа и лексикографа Владимира
Ивановича Даля (1801-1872). Как справедливо писал сам Даль в "Напутном" к первому
изданию: "Разбор и оценка других изданий должны окончиться прямым или
косвенным, скромным признанием, что наше всех лучше".
Предисловие
...и все бы ничего, но получил на свои опусы "отклики" читателей, точнее "рекламации", точнее – "претензии" читателей по поводу "ненадлежащего качества
поставляемого товара".
Клиент пошел грамотный и разборчивый!
И кроме брошюрок, читающий еще и
Законы
Российской
Федерации!
Например, Закон РФ от 07.02.1992 №
2300-1 (ред. от 13.07.2015) "О защите
прав
потребителей",
который
устанавливает права потребителя на
покупку товаров, услуг или работ
надлежащего качества и не несущих
какой-либо опасности для жизни и
здоровья человека, принадлежащего ему
имущества и окружающей среде. Кроме
того согласно этому закону, потребитель
вправе получить всю интересующую его
информацию
относительно
товаров,
услуг, работ и их изготовителей,
продавцов и исполнителей...
Вот и потребовали, а теперь "изготовитель" постарается предоставить читателям
"всю интересующую его информацию относительно товаров". Это можно считать
долгом автора, а "долг платежом красен". Правда, есть в этой пословице скрытая
двусмысленность: давать в долг с надеждой на навар от платежа или смиренно
ожидать желаемую отдачу долга... Это изречение обычно используется, когда
подразумевается ответная благодарность человека взявшего в долг. В рассматриваемом
случае возникает некоторая непонятность, несоответствие, например, ситуациям,
приведенным на картинках ниже...
3
Александр Сергеевич Пушкин
(1799—1837), повесть "Капитанская
дочка" (1836):
...Гринев поблагодарил Пугачева за
подаренного тем коня и тулуп:
Пугачев развеселился. "Долг
платежом красен!" - сказал он,
мигая и прищуриваясь.
Европа 2016: Бездомный и голодный...
4
Слава Богу, брошюрки читают не бездомные и совсем не голодные люди, но
попытаюсь смягчить удар "потребителей" по трем пунктам:
1. не сказано ничего о функциях Бесселя;
2. не приведена подпрограмма приведения матрицы к форме Хессенберга;
3. не сказано про алгоритмы управления промышленными объектами.
Все указанные претензии, так или иначе, могут быть истолкованы, как
относящиеся к большой теме "Теория Автоматического Управления". Об этом, конечно
стоит написать, но в следующей брошюрке (или брошюрках), а пока только
"лапидарий" — описание различных примечательных камней...
"Никто не обнимет необъятного" - Козьма Прутков...
Содержание
1.
Уравнение и функции Бесселя . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1. Цилиндрические функции Бесселя . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Решение дифференциального уравнения Бесселя . . . . . . . . . . . . . . . . .
5
5
13
2.
Приведение матрицы к верхней форме Хессенберга . . . . . . . . . . . . . . .
18
3.
Идентификация объектов управления . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1. Трогательное введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2. Предложение для заинтересовавшихся . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3. Определение параметров методом Симою . . . . . . . . . . . . . . . . . . . . . . .
26
26
30
32
Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
5
1. Уравнение и функции Бесселя.
На самом деле, цель энциклопедии — собрать знания,
рассеянные по свету, привести их в систему, понятную
для людей ныне живущих, и передать тем, кто придёт
после нас, с тем, чтобы труд предшествующих веков
не стал бесполезным для веков последующих, и чтобы
наши потомки, обогащённые знаниями, стали добрее и
счастливее, и чтобы мы не канули в вечность, не сумев
послужить грядущим поколениям.
Дени Дидро (1713-1784)
французский философ-просветитель.
Серьезный читатель обиделся на то, что в брошюре "Фрагмент 13: Решение
дифференциальных уравнений" слишком мало внимания было уделено
дифференциальным уравнениям, в частности, дифференциальному уравнению Бесселя
и цилиндрическим функциям Бесселя. Постараюсь скрасить обиду... а "спасибо" нужно
сказать Дени Дидро за великолепные слова: "...чтобы труд предшествующих веков не
стал бесполезным для веков последующих"!
1.1. Цилиндрические функции Бесселя.
Исходя из того, что читающих материалы по вынесенной в заголовок теме в
любой из известных мне книг на несколько порядков меньше, чем обращающихся за
справкой к "ВикипЕдии", решил начать с цитирования этого неиссякаемого источника
знаний. Это вполне объяснимо, если учесть, что слово "Wikipedia" получено слиянием
двух понятий: "wiki" - быстро и "encyclopedia" – справочник. Сегодня она является
общедоступным, универсальным, пополняемым мультиязычным интернет-ресурсом
(принадлежит американской некоммерческой организации "Фонд Викимедиа") с
бесплатным и безлимитным доступом. Начнем:
https://ru.wikipedia.org/wiki/Функции_Бесселя
Функции Бесселя в математике — семейство функций, являющихся
каноническими решениями дифференциального уравнения Бесселя:
Функции Бесселя впервые были определены швейцарским математиком
Даниилом Бернулли, а названы в честь Фридриха Бесселя.
В настоящей брошюре, возможно временно, будем полагать обозначения:
v – порядок функции является целым числом,
x - аргумент уравнения является вещественным числом.
6
https://ru.wikipedia.org/wiki/Бессель,_Фридрих_Вильгельм
Фридрих Вильгельм Бессель (1784—1846)
Фридрих Вильгельм Бессель (Friedrich Wilhelm
Bessel; 1784—1846) — немецкий математик,
астроном и геодезист, ученик Карла Фридриха
Гаусса. Не обучавшись в гимназии и
университете, получил докторскую степень
Гёттингенского университета.
Оставил
заметный
след
в
России:
мемориальная
мраморная
плита
в
Калининграде (бывш. Кёнигсберг) на холме
близ пересечения ул. Бесселя и ул. Генерала
Галицкого. Заслуги ученого были высоко
оценены избранием в члены многих академий,
в том числе Берлинскую (1812), и в
иностранные почетные члены Петербургской
академии наук (1814)...
Хорошо разработанная теория рассматриваемых функций, наличие подробных
таблиц и широкая область применений служат достаточным основанием для того,
чтобы отнести цилиндрические функции к числу наиболее важных специальных
функций. Приложения цилиндрических функций носят весьма разнообразный характер
— от абстрактных проблем теории чисел и теоретической астрономии до конкретных
прикладных проблем математической физики и техники. Для рассмотрения многих
проблем, связанных с применением цилиндрических функций, достаточно
ограничиться изучением специального класса этих функций, который соответствует
случаю, когда порядок функций v в уравнении Бесселя равен нулю или целому
положительному числу. Исследование данного класса носит более скромный характер,
чем теория, относящаяся к произвольным значениям v, и может служить хорошим
введением в эту общую теорию. Простейшими функциями рассматриваемого класса
являются функции Бесселя порядка нуль и единица: J0(z) b J1(z). При этом можно
показать, что функции Бесселя других порядков выражаются через эти две функции.
Они вполне могут использоваться для вычисления, например, определенных
интегралов. Кроме того, функции Бесселя первого и второго рода образуют полную
ортогональную систему функций. Это значит, что такие функции могут использоваться
для построения разложений прочих функций в ряды.
Для детального знакомства с теорией бесселевых функций можно посоветовать
две замечательные книги:
1. Трактат английского математика Джорджа Невилла Ватсона (George Neville
Watson; 1886-1965) Теория Бесселевых функций [02], в котором стоит обратить
внимание на главу II "Бесселевы функции с целыми индексами".
2. Фундаментальный Справочник по специальным функциям с формулами,
графиками и математическими таблицами под редакцией Абрамовица М. и
Стигана И. [03], который является переводом справочника, изданного в 1964
7
году в США. Сведения по обсуждаемой теме приведены в главе 9 "Функции
Бесселя целого порядка".
Материал, подаренный авторами в этих книгах, используется, особенно в
последние, не слишком плодотворные годы, для написания "новых" книг и, к
сожалению, не всегда с упоминанием первоисточников... Из "свежих" книг мне
понравилась отечественное пособие сотрудника кафедры высшей математики МФТИ
Владимира Ивановича Зубова [04], в котором на доступном уровне излагаются основ
теории функций Бесселя первого рода...
Приведенная ниже программа на языке FreeBASIC – ординарная переработка
программы (Алгоритм 5б) на языке Algol из нестареющей Библиотеки алгоритмов [05].
Там же была приписка авторов: значения функций можно вычислять с помощью более
современного алгоритма 236а.
'
P R O G R A M "CA_BesselJv"
'
10.05.2012
'------------------------------------------------------------'
Функция Бесселя первого рода целого порядка Jv(x)
'------------------------------------------------------------' Библиотека алгоритмов 1б-50б. Справочное пособие. Вып. 1.
' cтр. 19: Алгоритм 5б
' Вычисление функции Бесселя первого рода разложением в ряд.
'
v - порядок функции (целое число),
'
x - аргумент уравнения (вещественное число).
'
#Lang "FB"
' режим совместимости с FreeBASIC
'
Declare Function BesselJv(v As Integer, x as Double, Eps As
Double) As Double
'
' Основной модуль программы
Dim As Integer v
Dim As Double x, J0, J1, Eps
Open Cons For Output As #2
' вывод на консоль
'Open "CA_BesselJv.res" For Output As #2 ' вывод в файл
Eps = 0.000000001
'абсолютная погрешность
Print #2,: Print #2, " Bessel functions of the first kind:"
Print #2, "
x
J0
J1"
'
For x = -4 To 4 Step 0.25
v = 0
J0 = BesselJv(v, x, Eps)
v = 1
J1 = BesselJv(v, x, Eps)
Print #2, Using "
###.### ###.###########
###.###########"; x; J0; J1
Next x
Close #2
Print: Print " Press to complete!!!"
Sleep
8
'
' Подпрограмма вычисления функций Бесселя первого рода
Function BesselJv(v As Integer, x as Double, Eps As Double) As
Double
Dim As Integer s, t
Dim As Double k, term, sum, snfac
s = 0: x = x/2
sum = 1: term = 1: snfac = 1: k = 1
If x = 0 And v = 0 Then GoTo fin ' выход из подпрограммы
For t = 2 To v: snfac = snfac*t: Next t
k = x^v/snfac: x = x*x
iter: s = s + 1: term = -term*x/(s*(s + v))
sum = sum + term
If Abs(term) > eps Then GoTo iter ' большая пограшность
fin: BesselJv = k*sum
End Function
' Результат выполнения программы:
' Bessel functions of the first kind:
'
x
J0
J1
'
-4.000
-0.39714980986
0.06604332803
'
-3.750
-0.40140605495
-0.03322934913
'
-3.500
-0.38012773999
-0.13737752738
'
-3.250
-0.33275080217
-0.24111968802
'
-3.000
-0.26005195490
-0.33905895853
'
-2.750
-0.16414142781
-0.42597230295
'
-2.500
-0.04838377647
-0.49709410246
'
-2.250
0.08274985129
-0.54837835665
'
-2.000
0.22389077915
-0.57672480776
'
-1.750
0.36903253019
-0.58015619763
'
-1.500
0.51182767173
-0.55793650791
'
-1.250
0.64590608527
-0.51062326032
'
-1.000
0.76519768656
-0.44005058575
'
-0.750
0.86424227517
-0.34924360217
'
-0.500
0.93846980724
-0.24226845767
'
-0.250
0.98443592930
-0.12402597732
'
0.000
1.00000000000
0.00000000000
'
0.250
0.98443592930
0.12402597732
'
0.500
0.93846980724
0.24226845767
'
0.750
0.86424227517
0.34924360217
'
1.000
0.76519768656
0.44005058575
'
1.250
0.64590608527
0.51062326032
'
1.500
0.51182767173
0.55793650791
'
1.750
0.36903253019
0.58015619763
'
2.000
0.22389077915
0.57672480776
'
2.250
0.08274985129
0.54837835665
'
2.500
-0.04838377647
0.49709410246
'
2.750
-0.16414142781
0.42597230295
'
3.000
-0.26005195490
0.33905895853
'
3.250
-0.33275080217
0.24111968802
'
3.500
-0.38012773999
0.13737752738
'
3.750
-0.40140605495
0.03322934913
'
4.000
-0.39714980986
-0.06604332803
9
Ниже приведены графики функций Бесселя первого рода Jv, построенные по
значениям, полученным по программе:
На представленных графиках цветами линий обозначено:
v = 0 – красная линия;
v = 1 – синяя линия.
Думаю, для пытливых читателей, это уже хороший повод экспериментально
проверить некоторые теоретические утверждения...
Пожалуйста, сообщите о результатах!
Для любителей "прецизионных" вычислений ниже приведена программа
подсчета значений функции Бесселя первого рода нулевого порядка:
10
'
P R O G R A M "CA_BesselJ0"
'
10.05.2012
'------------------------------------------------------------'
Функция Бесселя первого рода нулевого порядка J0(x)
'------------------------------------------------------------' Справочник по специальным функциям, с формулами, графиками и
' математическими таблицами. Абрамовиц М., Стиган И. М.: 1979
' Стр. 177: Глава 9. Функции Бесселя целого порядка.
#Lang "FB"
' режим совместимости с FreeBASIC
#Include "crt/math.bi" ' математические функции
' объявление используемой подпрограммы-функции
Declare Function BesselJ0(X as Double) As Double
'
' Основной модуль программы
Dim As Double X, J0
Open Cons For Output As #2
' вывод на консоль
'Open "CA_BesselJ0.res" For Output As #2 ' вывод в файл
Print #2,: Print #2, " Bessel functions of the first kind:"
Print #2, "
x
J0"
For X = 0 To 4 Step 0.25
J0 = BesselJ0(X)
Print #2, Using "
###.### ###.################"; X; J0
Next X
Close #2
Print: Print " Press to complete!!!"
Sleep
'
' Вычисление значения функции Бесселя первого рода
Function BesselJ0(X as Double) As Double
Dim As Integer I, S
Dim As Double Y, Z, Q1, Q2, Q3, FX, X1, X2, X3, X4, XLG
Dim As Double A(0 To 17), C(0 To 17), D(0 To 17)
' Константа M_Pi_4 = Pi/4 = 0.78539816339744830962
' Коэффициенты формул разложения функции:
A(0) = -0.17e-18
A(1) = 0.1222e-16
A(2) = -0.75885e-15
A(3) = 0.4125321e-13
A(4) = -0.194383469e-11
A(5) = 0.7848696314e-10
A(6) = -0.267925353056e-8
A(7) = 0.7608163592419e-7
A(8) = -0.176194690776215e-5
A(9) = 0.3246032882100508e-4
A(10) = -0.4606261662062751e-3
A(11) = 0.4819180069467605e-2
A(12) = -0.3489376941140889e-1
A(13) = 0.1580671023320973e0
A(14) = -0.3700949938726498e0
A(15) = 0.2651786132033368e0
A(16) = -0.8723442352852221e-2
A(17) = 0.3154559429497802e0
11
'
C(0) = -0.1e-19
C(1) = 0.2e-19
C(2) = -0.11e-18
C(3) = 0.55e-18
C(4) = -0.288e-17
C(5) = 0.1631e-16
C(6) = -0.10012e-15
C(7) = 0.67481e-15
C(8) = -0.506903e-14
C(9) = 0.4326596e-13
C(10) = -0.4045789e-12
C(11) = 0.516826239e-11
C(12) = -0.7864091377e-10
C(13) = 0.163064646352e-8
C(14) = -0.5170594537606e-7
C(15) = 0.307518478751947e-5
C(16) = -0.5365220468132117e-3
C(17) = 0.1998920698695037e+1
'
D(0) = 0.1e-19
D(1) = -0.3e-19
D(2) = 0.13e-18
D(3) = -0.62e-18
D(4) = 0.311e-17
D(5) = -0.1669e-16
D(6) = 0.9662e-16
D(7) = -0.60999e-15
D(8) = 0.425523e-14
D(9) = -0.3336328e-13
D(10) = 0.30061451e-12
D(11) = -0.320674742e-11
D(12) = 0.4220121905e-10
D(13) = -0.72719159369e-9
D(14) = 0.1797245724797e-7
D(15) = -0.74144984110606e-6
D(16) = 0.683851994261165e-4
D(17) = -0.3111170921067402e-1
'
Y = Abs(X) ' Четная функция
S = 1: If X < 0 Then S = -1
Z = Y*0.125
If Z > 1 Then GoTo L3
If Z = 0 Then GoTo L2
' если 0 < Z <= 1
X2 = 4*Z*Z - 2
Q3 = 0: Q2 = 0
For I = 0 To 17
Q1 = Q2: Q2 = Q3
Q3 = X2*Q2-Q1+A(I)
Next I
BesselJ0 = (Q3-Q1)*0.5
12
GoTo L6
'
L2: ' если Z = 0
BesselJ0 = 1
GoTo L6
'
L3: ' если Z > 1
Z = 1/Z
X2 = 4*Z*Z - 2
Q3 = 0: Q2 = 0
For I = 0 To 17
Q1 = Q2: Q2 = Q3
Q3 = X2*Q2-Q1+C(I)
Next I
X1 = (Q3-Q1)*0.5
Q3 = 0: Q2 = 0
For I = 0 to 17
Q1 = Q2: Q2 = Q3
Q3 = X2*Q2-Q1+D(I)
Next I
FX = (Q3-Q1)*0.5
X2 = Cos(Y-M_Pi_4)
X3 = Sin(Y-M_Pi_4)
X4 = 1/(Sqr(M_Pi_2)*Sqr(Y)) ' M_Pi_2 = Pi/2
FX = FX*Z
BesselJ0 = S*X4*(X1*X2-FX*X3)
L6:
End Function
' Результат выполнения программы:
' Bessel functions of the first kind:
'
x
J0
Jc
' 0.000
1.0000000000000000
1.0000000000000000
' 0.250
0.9844359292958528
0.9844359292958527
' 0.500
0.9384698072408129
0.9384698072408129
' 0.750
0.8642422751666488
0.8642422751666486
' 1.000
0.7651976865579667
0.7651976865579666
' 1.250
0.6459060852712850
0.6459060852712853
' 1.500
0.5118276717359181
0.5118276717359181
' 1.750
0.3690325301851508
0.3690325301851508
' 2.000
0.2238907791412357
0.2238907791412357
' 2.250
0.0827498512887339
0.0827498512887340
' 2.500 -0.0483837764681980
0.0483837764681979
' 2.750 -0.1641414278085138 -0.1641414278085137
' 3.000 -0.2600519549019336 -0.2600519549019334
' 3.250 -0.3327508021706116 -0.3327508021706115
' 3.500 -0.3801277399872635 -0.3801277399872634
' 3.750 -0.4014060549361745 -0.4014060549361743
' 4.000 -0.3971498098638475 -0.3971498098638474
Примечание: третий столбец (Jc) заполнен контрольными значениями.
Результат ожидаемо отличный! Пользуйтесь!
13
1.2. Решение дифференциального уравнения Бесселя.
В пункте 1.1. указывалось, что значения функций Бесселя могут быть получены,
как решение задачи Коши для дифференциального уравнения Бесселя. Уравнение
представляет собой линейное обыкновенное дифференциальное уравнение второго
порядка, которое может быть приведено к системе двух дифференциальных уравнений
первого порядка.
Ниже приведена программа для решения поставленной задачи на языке
FreeBASIC, которая уже использовалась в брошюрах этой серии. В программе
используется метод Рунге—Кутты (распространено неправильное название Методы
Рунге—Кутта) четвертого порядка из семейства численных алгоритмов решения
обыкновенных дифференциальных уравнений и их систем. Данные итеративные
методы явного и неявного приближённого вычисления были разработаны около 1900
года немецкими математиками К. Рунге и М. Куттой. Программа получена ординарной
переработкой программы из книги по численным методам [06]. В книге много
маленьких программ по многим разделам вычислительной математике. Для
дальнейшего ознакомления с проблематикой решения дифференциальных уравнений
можно рекомендовать книгу Александра Ивановича Егорова [07], в которой
рассмотрены основные направления теории обыкновенных дифференциальных
уравнений и практические методы решения таких уравнений. Значительная часть книги
содержит стандартный учебный материал по курсу обыкновенных дифференциальных
уравнений. Кроме того, рассматриваются матричные дифференциальные уравнения,
основы теории устойчивости по Ляпунову, основы теории периодических решений
нелинейных уравнений, теория уравнений с разрывной правой частью
(дифференциальные включения) и применение теории групп Ли к решению
обыкновенных дифференциальных уравнений. Будет полезен и двухтомник "Решение
обыкновенных дифференциальных уравнений" [08], посвященный "нежестким
задачам" и [09], описывающий решение жестких и дифференциально-алгебраических
задач.
'
P R O G R A M "CA_RnKt4_method"
'
21.02.2013
'------------------------------------------------------------'
Решение системы обыкновенных дифференциальных уравнений
'------------------------------------------------------------'
Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик,
'
Фортран и Паскаль. Томск: МП "РАСКО", 1991. 272 с: ил.
'
Программа 6.4В: Метод Рунге-Кутты (Runge and Kutta
methods)
#Lang "FB"
' Режим совместимости с FreeBASIC
' В программе используются следующие подпрограммы:
' подпрограмма-функция вычисления правых частей
Declare Sub RPart(N As Double, X As Double, P As Double, Y()
As Double, F() As Double)
' подпрограмма шага Рунге-Кутты 4-го порядка
Declare Sub RnKt4(N As Integer, X As Double, P As Double, H As
Double, Y() As Double, F() As Double)
'
' Основной модуль программы
14
Dim As Integer N, M, I
N = 2
' порядок системы (размер матрицы)
Dim As Double P, X, X0, X9, H, Y(N), Y0(N), F(N)
P = 0
' порядок функции
X0 = 0.0: X9 = 2.5: M = 10 ' исходные данные
Y0(1) = 1.0: Y0(2) = 0.0
' значения при X0 = 0
' контрольная точка x = 2.5 y1 = -0.048383776468 y2 = 0.497094102464
Open Cons For Output As #2
'Open "CA_RnKt4_method.res" For Output As #2
Print #2, " Solution of differential equations of RungeKutta-4 method"
Do
Y(1) = Y0(1): Y(2) = Y0(2)
H = (X9-X0)/M
Print #2,: Print #2, Using " M = ######## Step =
##.#######"; M; H
Print #2, "
i
x
J0(x)
J0'(x)"
X = X0
For I = 1 To M
X = X + H
RnKt4(N, X, P, H, Y(), F()) ' Подпрограмма
интегрирования
If I Mod M/10 = 0 Then
Print #2, Using " ######### ##.###### ##.########
##.########"; I; X; Y(1); Y(2)
EndIf
Next I
M = 10*M
Loop While M <= 1000000
Close #2
Print: Print " Press to complete!!!"
Sleep
'
' Подпрограмма-функция вычисления правых частей
Sub RPart(N As Double, X As Double, P As Double, Y() As
Double, F() As Double)
F(1) = Y(2)
' Правые части уравнений
F(2) = ((P/X)^2-1)*Y(1) - Y(2)/X
End Sub
' Подпрограмма шага Рунге-Кутты 4-го порядка
Sub RnKt4(N As Integer, X As Double, P As Double, H As Double,
Y() As Double, F() As Double)
Dim As Integer I, J
Dim As Double H1, H2, Y0(N), Y1(N), Q
H1 = 0: H2 = H/2
For I = 1 To N: Y0(I) = Y(I): Y1(I) = Y(I): Next I
For J = 1 To 4
RPart(N, X+H1, P, Y(), F()) ' Вычисление правых частей
If J = 3 Then
H1 = H
Else
15
H1 = H2
EndIf
For I = 1 To N
Q = H1*F(I): Y(I) = Y0(I) + Q
If J = 2 Then Q = 2*Q
Y1(I) = Y1(I) + Q/3
Next I
Next J
For I = 1 To N: Y(I) = Y1(I): Next I
End Sub
' Результат выполнения программы:
' Solution of differential equations
'
' M =
10 Step = 0.2500000
'
i
x
J0(x)
'
1
0.250000
0.97469980
'
2
0.500000
0.91083368
'
3
0.750000
0.81664540
'
4
1.000000
0.69826623
'
5
1.250000
0.56196572
'
6
1.500000
0.41440842
'
7
1.750000
0.26250398
'
8
2.000000
0.11313506
'
9
2.250000 -0.02715116
'
10
2.500000 -0.15244778
'
' M =
100 Step = 0.0250000
'
i
x
J0(x)
'
10
0.250000
0.98207463
'
20
0.500000
0.93332274
'
30
0.750000
0.85648083
'
40
1.000000
0.75516077
'
50
1.250000
0.63405094
'
60
1.500000
0.49869759
'
70
1.750000
0.35522539
'
80
2.000000
0.21002681
'
90
2.250000
0.06943925
'
100
2.500000 -0.06057197
'
' M =
1000 Step = 0.0025000
'
i
x
J0(x)
'
100
0.250000
0.98414009
'
200
0.500000
0.93787989
'
300
0.750000
0.86338506
'
400
1.000000
0.76411274
'
500
1.250000
0.64464335
'
600
1.500000
0.51044480
'
700
1.750000
0.36759189
'
800
2.000000
0.22245626
'
900
2.250000
0.08138364
'
1000
2.500000 -0.04962432
'
of Runge-Kutta-4 method
J0'(x)
-0.18565538
-0.32032589
-0.42919497
-0.51368828
-0.57228391
-0.60357301
-0.60707658
-0.58351154
-0.53484145
-0.46420533
J0'(x)
-0.13517003
-0.25320677
-0.35911402
-0.44830568
-0.51685657
-0.56186624
-0.58162801
-0.57571440
-0.54498944
-0.49155045
J0'(x)
-0.12523595
-0.24340055
-0.35024139
-0.44086770
-0.51122329
-0.55829442
-0.58025942
-0.57657375
-0.54798621
-0.49648601
16
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
M =
10000 Step =
i
x
1000
0.250000
2000
0.500000
3000
0.750000
4000
1.000000
5000
1.250000
6000
1.500000
7000
1.750000
8000
2.000000
9000
2.250000
10000
2.500000
0.0002500
J0(x)
0.98440514
0.93840947
0.86415519
0.76508788
0.64577861
0.51168834
0.36888762
0.22374669
0.08261281
-0.04850803
J0'(x)
-0.12414796
-0.24238193
-0.34934329
-0.44013195
-0.51068271
-0.55797160
-0.58016571
-0.57670883
-0.54833826
-0.49703243
100000 Step =
i
x
10000
0.250000
20000
0.500000
30000
0.750000
40000
1.000000
50000
1.250000
60000
1.500000
70000
1.750000
80000
2.000000
90000
2.250000
100000
2.500000
0.0000250
J0(x)
0.98443283
0.93846375
0.86423355
0.76518669
0.64589332
0.51181373
0.36901803
0.22387636
0.08273614
-0.04839620
J0'(x)
-0.12403818
-0.24227981
-0.34925357
-0.44005872
-0.51062920
-0.55794001
-0.58015714
-0.57672320
-0.54837433
-0.49708792
0.0000025
J0(x)
0.98443562
0.93846920
0.86424140
0.76519659
0.64590481
0.51182628
0.36903108
0.22388934
0.08274848
-0.04838502
J0'(x)
-0.12402720
-0.24226959
-0.34924460
-0.44005140
-0.51062385
-0.55793686
-0.58015629
-0.57672465
-0.54837795
-0.49709348
M =
M =
1000000 Step =
i
x
100000
0.250000
200000
0.500000
300000
0.750000
400000
1.000000
500000
1.250000
600000
1.500000
700000
1.750000
800000
2.000000
900000
2.250000
1000000
2.500000
Здесь продемонстрировано влияние величины шага вычислений на результат
работы. Проблема необходимой и достаточной точности вычислений остается за
читателями. Одно замечу: при нынешних средствах вычислительно техники не следует
экономить на форматах представления числовых данных!
На рисунке ниже приведены фрагменты функции и ее производной:
17
Фрагмент функции Бесселя первого рода нулевого порядка J0(x) и ее производной.
18
2. Приведение матрицы к верхней форме Хессенберга.
Невозможно объяснить, что такое Матрица...
Ты должен увидеть это сам. Не поздно отказаться.
Потом пути назад не будет. Примешь синюю таблетку —
и сказке конец. Ты проснешься в своей постели и
поверишь, что это был сон. Примешь красную таблетку —
войдешь в страну чудес. Я покажу тебе, насколько
глубока кроличья нора.
"Матрица" (The Matrix)
культовый научно-фантастический боевик
братьев Вачовски, вышел на экраны 31 марта 1999.
...не смотрел и смотреть не буду! Просмотрел весь институтский учебник
по матричному исчислению, но так и не нашел таких операций над матрицами,
как "Революция" и "Перезагрузка"...
Молодежь ошибочно считает "максимум" - радиостанцией, "матрицу" фильмом, а "ориентацию" - половой. Но что интересно, на Мехмате в эти
понятия вкладывают точно такой же смысл...
19
Мы тоже не пальцем деланы – будем творить квадратные матрицы – чтобы
цифирки не разбегались!
Замечание читателя касается брошюры "Фрагмент 17. Вычислительные задачи
алгебры" в п.п. 4.2.2. Подпрограмма "MHQR.bas" которой сказано, что процедура
может быть использована для отыскания всех собственных значений произвольной
действительной матрицы ее следует преобразовать к верхней форме Хессенберга одной
из процедур elmHes, dirHes или ortHes, но ни одна из них не была представлена.
Досадно, конечно...
Замечу, что поднятая в этом пункте тема имеет непосредственное отношение к
теории автоматического управления (ТАУ), но поскольку брошюры предназначаются
для старшеклассников и младшекурсников речь пойдет исключительно о линейных
стационарных непрерывных системах. Здесь уместно вспомнить слова специалиста в
области теории автоматического регулирования и систем автоматического управления
Александра Аркадьевича Красовского (1921—2003): "Периодизация развития теории
автоматического управления, как и многих других, относительно новых областей
пауки, не является установившейся и общепринятой" [10]. До пятидесятых годов
сложилась классическая "Теория автоматического регулирования" (ТАР),
рассматривающая процессы в системе "объект — регулятор". Основы этой теории были
заложены И.А.Вышнеградским, А.Б.Стодолой, Д.К.Максвеллом... Под влиянием
потребностей автоматизации управления технологическими процессами в пятидесятых
годах начала формироваться "Теория автоматического управления" (ТАУ), появились
известные работы Л.С. Поптрягипа, Р. Беллмана, Р. Калмана. Пока не существует
общепринятого определения понятия "Современная теории автоматического
управления" (CTAУ) - авторы не могут сформулировать ее характерные признаки!
Многие из них называют описание процессов в пространствах состояний, хотя и в
классической теории широко применялось описание движения в фазовом пространстве
(особенно — на фазовой плоскости).
В упомянутой книге [10] Красовского отмечает: Классическую ТАР в основном
создавали инженеры для инженеров и лишь частично математики для инженеров.
СТАУ создают в основном математики для инженеров и во все большей мере
математики для математиков. Последнее с точки зрения практики вызывает
определенное беспокойство... Главное негативное влияние на практическое внедрение
методов СТАУ оказывает масса оторванных от практических потребностей и
возможностей работ и даже направлений, интересных в математическом отношении, но
бесплодных в отношении современных приложений. Нельзя отрицать право на
существование математической СТАУ как раздела математики, развивающегося по
собственным законам и находящего применение по мере возникновения
соответствующих потребностей. Однако такая математическая СТАУ должна быть
достаточно четко выделена по отношению к прикладной СТАУ".
Но, оставим эти проблемы академикам...
Посмотрим на некоторые особенности "математической СТАУ":
https://ru.wikipedia.org/wiki/Пространство_состояний_(теория_управления)
Пространство состояний в теории управления один из основных методов
описания поведения динамической системы. Движение системы в пространстве
20
состояний отражает изменение её состояний, а модель системы - набор переменных
входа, выхода и состояния, связанных между собой дифференциальными уравнениями
первого порядка, которые записываются в матричной форме. Иными словами, модель
системы представляется в виде связки: "вход - состояние - выход".
Как вяжутся между собой понятия "матрица", "пространство состояний"
применительно к "Теории автоматического управления" можно узнать из книги [11],
посвященной современной теории управления и разнообразным приложениям ее
результатов. Она состоит из четырех частей, в которых излагаются вопросы
устойчивости, управления в детерминированных и стохастических системах, метода
расчета систем управления. Вскользь упоминается, что для решения многих задач ТАУ
используются алгоритмы отыскания собственных значений матрицы, которая
представляет собой модель системы. О решении полной и частичной проблем
собственных значений можно узнать в книге [12] или в книге [13] по вычислительным
методам линейной алгебры.
К слову, следует заметить, что в линейной алгебре "центр тяжести"
рассматриваемых проблем мигрирует – достаточно сравнить перечень задач в 1960
году, представленный в книге [12]:
 решение систем линейных алгебраических уравнений,
 вычисление определителей матриц,
 нахождение обратной матрицы,
 нахождение собственных значений и собственных векторов матриц,
и перечень задач в 1997 году, представленный в книге [14]:
 решение систем линейных алгебраических уравнений;
 линейная задача о наименьших квадратах;
 нахождение собственных значений и собственных векторов матриц;
 нахождение сингулярных чисел и сингулярных векторов матриц.
Указанную тенденцию следует учитывать при планировании перечня необходимых
программ для нужд ТАУ!
Когда говорят о вычислении собственных значений матрицы, то различают
полную и частичную проблему собственных чисел. В полной проблеме вычисляют
все собственные числа и соответствующие им собственные векторы матриц, а под
частичной проблемой обычно понимают задачу нахождения одного или нескольких
собственных чисел и соответствующих им собственных векторов. Методы решения
указанных задач по типу примененной вычислительной схемы разделяют на прямые
(точные, например, метод Крылова) и итерационные (приближенные, например, метод
Якоби). В прямых методах обычно строится правило вычисления коэффициентов
характеристического многочлена матрицы, затем уже собственные значения находятся
как корни характеристического многочлена по какому-либо численному методу. После
этого вычисляются собственные векторы, что обычно считается простой задачей. В
итерационных методах коэффициенты характеристического уравнения не вычисляют, а
составляют некоторые итерационные последовательности, позволяющие найти одно
или несколько, а иногда и все собственные значения исходной матрицы А.
Итерационные методы почти всегда более трудоемкие, однако они надежнее прямых
методов, так как менее чувствительны к ошибкам округления.
Интерес к нахождению собственных значений матрицы связан с тем, что
нахождение собственных значений матрицы сводится к решению алгебраического
21
характеристического уравнения и любой алгебраический полином, является
характеристическим полиномом для некоторой матрицы. А если вспомнить
теорему Абеля-Руффини о том, что для алгебраических полиномов степени выше 4 не
существует прямых (конечных) методов нахождения корней, в которых выражения для
этих корней даются в виде композиций четырёх арифметических действий и операции
взятия корня, то интерес к нахождению собственных значений еще более
увеличивается. Появляется естественная идея предварительно привести матрицу, для
которой решается проблема собственных значений, к некоторой, по возможности,
простейшей форме, для которой собственные значения и/или собственные векторы
могут быть найдены проще, чем для исходной матрицы.
Есть три специальных класса матриц, которые имеют большое значение в
решении проблемы собственных значений и для которых выбор ведущего элемента
осуществляется много проще.
Класс I. Матрицы Хессенберга.
Возможно использование двух вариантов матрицы Хессенберга:
1) верхняя матрица Хессенберга, если
aij = 0 при i >= j + 2;
2) нижняя матрица Хессенберга, если
aij = 0 при j >= i + 2.
Класс II. Трехдиагональные матрицы.
Матрица является трехдиагональной, если
aij = 0 при |j – i| > 1.
Класс III. Положительно определенные симметричные матрицы.
Симметричная матрица называется положительно определенной, если все ее
главные миноры положительны.
На этом "введение" заканчивается – переходим к матрицам Хессенберга:
http://ru.wikipedia.org/wiki/Матрица_Хессенберга
Рассматриваемые матрицы названы в честь немецкого математика и инженера
Карла Адольфа Хессенберга (Karl Adolf Hessenberg; 1904-1959). В линейной алгебре,
матрицами Хессенберга называют "почти" треугольные матрицы.
Верхняя матрица Хессенберга — это квадратная матрица, у которой все элементы
ниже первой поддиагонали равны нулю.
Нижняя матрица Хессенберга — это квадратная матрица, у которой все члены выше
первой наддиагонали равны нулю.
Матрицы, которые являются одновременно как нижними, так и верхними матрицами
Хессенберга — это трёхдиагональные матрицы.
Диагональные элементы указанных матриц являются собственными значениями...
Матрицы Хессенберга служат отправными формами для многих алгоритмов,
например, вычисления собственных значений, поскольку нулевые значения уменьшают
сложность задачи. Существует несколько методов сведения матриц к матрицам
Хессенберга с теми же собственными значениями.
В замечательном справочнике Уилкинсона и Райнша [15] приведены алгоритмы
решения всех основных задач линейной алгебры, реализованные в виде процедур на
22
языке Алгол-60, в том числе несколько процедур приведения матрицы к верхней
форме Хессенберга:




ElmHes - предназначена для приведения произвольной действительной
матрицы к верхней форме Хессенберга с помощью действительных устойчивых
элементарных преобразований подобия.
DirHes - выполняет такое же приведение, что и процедура ElmHes, но с
меньшими ошибками за счет накопления скалярного произведения с удвоенной
точностью.
ComHes - позволяет осуществить приведение комплексной матрицы к верхней
форме Хессенберга с помощью элементарных устойчивых преобразований
подобия.
OrtHes - предназначена для приведения действительной матрицы к верхней
форме Хессенберга с помощью элементарных ортогональных преобразований.
Ниже приведена программа на языке FreeBASIC, реализующая алгоритм
процедуры ElmHes (is used to reduce a real matrix to upper-Hessenberg form by real
stabilized elementary similarity transformations). Пример взят из книги [15, стр. 309].
'
P R O G R A M "CA_ElmHes"
'
12.10.2012
'------------------------------------------------------------'
Приведение квадратной матрицы к верхней форме Хессенберга
'------------------------------------------------------------' Алгоритм и пример взяты из Справочника Уилкинсона и Райнша
' Обозначено:
'
А0()
- квадратная матрица (исходная, сохраняется);
'
H()
- квадратная матрица Хессенберга (результирующая);
'
NM
- порядок матрицы;
'
I
– номер строки;
'
J
– номер столбца;
'
VPrm
- вектор перестановок.
'
#Lang "FB"
' режим совместимости с FreeBASIC
' Объявление используемой подпрограммы
Declare Sub MHES(NM As Integer, A1() As Double, H() As Double,
VPrm() As Integer)
' Блок исходных данных
Data 4
Data 1.0, 2.0, 3.0, 5.0
Data 2.0, 4.0, 1.0, 6.0
Data 1.0, 2.0, -1.0, 3.0
Data 2.0, 0.0, 1.0, 3.0
'
' Основной модуль программы
Dim As Integer NM, I, J
Read NM
' Ввод порядка матрицы
Dim As Integer VPrm(NM)
Dim As Double A1(NM,NM), H(NM,NM)
' Ввод исходной матрицы
23
For I = 1 To NM
For J = 1 To NM: Read A1(I,J): Next J
Next I
'
Open Cons For Output As #2
'Open "CA_ElmHes.res" For Output As #2
Print #2,: Print #2, " Bringing the matrix to form
Hessenberg"
'
MHES(NM, A1(), H(), VPrm())
'
Print #2,: Print #2, " Resulting transformed matrix:"
Print #2,
For I = 1 To NM
For J = 1 To NM: Print #2, Using " ###.##### "; H(I,J);: Next
J
Print #2,
Next I
Print #2,
Print #2, " Vector permutations:"
Print #2,
For J = 1 To NM: Print #2, Using "
### "; VPrm(J);:
Next J
Print #2,
Close #2
Print: Print " Press to complete!!!"
Sleep
'
' Подпрограмма приведения к верхней форме Хессенберга
Sub MHES(NM As Integer, A() As Double, H() As Double, VPrm()
As Integer)
' В подпрограмме обозначено:
'
А0()
- квадратная матрица (исходная, сохраняется);
'
H()
- квадратная матрица Хессенберга (результирующая);
'
NM
- порядок матрицы;
'
I
– номер строки;
'
J
– номер столбца;
'
VPrm
- вектор перестановок.
'
Dim As Integer I, J, M
Dim As Double X, Y
' Сохранение исходной матрицы
For I = 1 To NM
For J = 1 To NM: H(I,J) = A(I,J): Next J
Next I
' Обнуление вектора перестановок
For I = 1 To NM: VPrm(I) = 0: Next I
'
For M = 2 To NM-1
'
' Определение максимального по модулю элемента
I = M: X = 0
24
For J = M TO NM
If Abs(H(J,M-1)) > Abs(X) Then X = H(J,M-1): I = J
Next J
VPrm(M) = I
'
If I <> M Then
' Перестановка строк и столбцов матрицы
If M-1 <= NM Then
For J = M-1 To NM
Y = H(I,J)
H(I,J) = H(M,J)
H(M,J) = Y
Next J
EndIf
For J = 1 To NM
Y = H(J,I)
H(J,I) = H(J,M)
H(J,M) = Y
Next J
EndIf
'
If X <> 0 Then
For I = M+1 To NM
Y = H(I,M-1)
If Y <> 0 Then
Y = Y/X
H(I,M-1) = Y
For J = M To NM
H(I,J) = H(I,J) - Y*H(M,J)
Next J
For J = 1 To NM
H(J,M) = H(J,M) + Y*H(J,I)
Next J
EndIf
Next I
EndIf
'
Next M
' Обнуление "лишних" элементов матрицы
For I = 2 To NM
For J = 1 To I-2
H(I,J) = 0
Next J
Next I
'
End Sub
' Результат выполнения программы:
25
'
'
'
'
'
'
'
'
Bringing the matrix to form Hessenberg
Resulting transformed matrix:
1.00000
8.50000
5.32143
3.00000
2.00000
10.50000
6.10714
1.00000
0.00000
-7.00000
-3.00000
0.00000
0.00000
0.00000
0.16071
-1.50000
Vector permutations:
0
2
4
0
Результат, конечно, ожидаемый...
Однако на этом "поле чудес" еще пахать и пахать!
...не следует близко к сердцу принимать слова комментария в последних строках
подпрограммы – элементы матрицы еще пригодятся...
Удачи!
26
3. Идентификация объектов управления.
Следующую претензию счел обоснованной – она касается сразу двух брошюрок:
Фрагмент 11. Обработка цифровых данных.
На рисунке приведена упрощенная схема замкнутой системы управления – это на тот
случай, если, Бог даст, доберемся до программирования алгоритмов управления
промышленного назначения...
Фрагмент 13. Решение дифференциальных уравнений.
Материал требует более детального рассмотрения. Например, представляют интерес
следующие вопросы:
- идентификация параметров моделей реальных технологических процессов;
- реакция системы на всевозможные возмущающие сигналы;
- разработка алгоритмов расчета параметров стандартных регуляторов.
3.1. Трогательное введение.
Обсуждение начнем, пожалуй, с цитирования "неиссякаемого источника
знаний", с ВикипЕдии - пусть она послужит нам Нитью Ариадны. Напомню: когда
Тесей (сын афинского царя Эгея) решился убить Минотавра (плод неестественной
любви Пасифаи, жены царя Миноса, к посланному Посейдоном быку), он получил от
пламенно любившей его Ариадны клубок ниток, выведший его из лабиринта, где
обитал Минотавр (применять нить научил её Дедал).
Тесей и Ариадна. Гравюра флорентийского мастера 15 века.
27
https://ru.wikipedia.org/wiki/Автоматизация
Автоматизация — одно из направлений научно-технического прогресса,
использующее саморегулирующие технические средства и математические методы с
целью освобождения человека от участия в процессах получения, преобразования,
передачи и использования энергии, материалов, изделий или информации, либо
существенного уменьшения степени этого участия или трудоёмкости выполняемых
операций. Объектами автоматизации являются следующие области человеческой
деятельности:
 производственные процессы;
 проектирование;
 организация, планирование и управление;
 научные исследования;
 обучение;
 бизнес-процессы.
Договоримся, что будем рассматривать лишь одну область, а именно
автоматизацию производственных процессов, наиболее интересную для творческих
натур. В ней автоматизация предназначена для частичного или полного исключения
человека из непосредственного участия в производственном процессе. Автоматизация
производства направлена на то, чтобы процессы получения, преобразования, передачи
и использования энергии, материалов и информации выполнялись автоматически.
Необходимость автоматизации обусловлена факторами прогресса: ростом
производительности и эффективности труда, увеличением мощности и быстродействия
машин, повышением точности и качества выполнения операций, сложности и
опасности технологических процессов, ростом ответственности за безопасность
человека и его среды обитания... При этом под автоматикой понимают отрасль науки и
техники, охватывающую теорию и принципы построения систем управления,
действующих без непосредственного участия человека. В узком смысле, автоматика —
совокупность методов и технических средств, исключающих участие человека при
выполнении операций конкретного процесса. Т.е. автоматизация это инженерный
процесс, запускаемый для реализации автоматики. Этот процесс (условно) можно
разделить на два главных направления:
 Анализ систем управления.
 Синтез систем управления.
Дальше Нить Ариадны тянется к статье:
https://ru.wikipedia.org/wiki/Теория_автоматического_управления
Теория автоматического управления (ТАУ) — научная дисциплина, которая
изучает процессы автоматического управления объектами разной физической природы.
При этом при помощи математических средств выявляются свойства систем
автоматического управления и разрабатываются рекомендации по их проектированию.
Конечно, здесь узнаем много полезного, но... продолжим;
https://ru.wikipedia.org/wiki/Идентификация_систем
Идентификация систем — совокупность методов для построения
математических моделей динамической системы по данным наблюдений.
Математическая модель в данном контексте означает математическое описание
поведения какой-либо системы или процесса в частотной или временной области, к
примеру, физических процессов (движение механической системы под действием силы
тяжести), экономического процесса (реакция биржевых котировок на внешние
28
возмущения) и т.п. В настоящее время эта область теории управления хорошо изучена
и находит широкое применение на практике.
Ну, не может все эти премудрости человек удержать в голове – нужно выбрать
что-то самое главное, а остальное "оставить на потом"... Для этой вот цели заглянул на
сайт "Alma Mater" (МИСиС) и с ужасом обнаружил, что направлением
"Автоматизация технологических процессов" занимается только одна кафедра!
Постучался на E-mail, попросил помочь хотя бы актуальными ссылками на
доступную в сети литературу по этому профилю (признаюсь - отстал от жизни)... В
тайне подумал: было бы прекрасно соблазнить кого-нибудь с кафедры пописать
пособие... Или их у кафедры уже перебор методичек? Но поддержки, к сожалению, не
нашел...
Крымский вал, дом 3: 4-ый корпус МИСиС
Слева виднеется главный корпус МИСиС, а направо дорожка в Парк Культуры...
Теперь уже трудно сказать, что главнее в становлении инженера...
Здесь прошли мои, наверное, лучшие годы!
А какие великолепные примеры объектов для автоматизации предлагает металлургия!
Ниже упомяну только два из них:
3.1.1. АО "Северсталь".
АО "Северсталь" — российская вертикально-интегрированная сталелитейная и
горнодобывающая компания, владеющая Череповецким металлургическим комбинатом
(Вологодская область), вторым по величине сталелитейным комбинатом России.
29
Надо быть настоящим металлургом, чтобы увидеть и оценить красоту этого
сооружения.
Пятая доменная печь Череповецкого металлургического комбината —
знаменитая "Северянка", самая крупная домна в бывшем СССР и в мире... была
задута 12 апреля 1986 года (в День космонавтики) в 5 часов 55 минут, а 13 апреля в 20
часов 15 минут выдала первый чугун. На сегодняшний день ДП №5 – крупнейшая
доменная печь Европы: Полезный объем — 5580 м3, высота центрального узла печи —
105 метров. Для обеспечения работы домны требуется ежегодно 1 млн. 700 тыс. тонн
кокса, около 2 млн. тонн окатышей, и свыше 4 млн. тонн агломерата. Проектная
производительность печи составляет до 4 миллионов тонн в год.
3.1.2. Институт атомных реакторов.
В 80 километрах от Ульяновска, на реке Черемшан, находится город
Димитровград с населением около 100 000 человек. Его главное предприятие —
Научно-исследовательский институт атомных реакторов (НИИАР), который был создан
в 1956 году по инициативе Курчатова. Изначально он был опытной станцией для
испытаний ядерных реакторов, но в настоящее время спектр направлений деятельности
значительно расширился. Сейчас в НИИАР испытывают различные материалы, чтобы
определить, как они себя ведут в условиях продолжительного радиоактивного
излучения, создают радионуклидные источники и препараты, которые применяют в
медицине и исследованиях, решают технические вопросы экологически чистых
технологий и просто ведут научную деятельность. В НИИАР работает около 3 500
сотрудников и 6 реакторов.
30
Дмитровград. Реакторный зал
и это тоже металлургия...
Вы можете представить себе цену в 27 000 000 долларов США за один грамм
вещества? Именно столько стоит радиоактивный элемент Калифорний-252. Дороже
только антиматерия, которая является самой дорогой субстанцией в мире (около 60
триллионов долларов за грамм антиводорода). На сегодняшний день в мире накоплено
всего 8 грамм Калифорния-252, а ежегодно производится не более 40 миллиграмм.
3.2. Предложение для заинтересовавшихся.
Если вдруг у кого-нибудь из читателей, прочитанное выше вызвало интерес, без
всякого стеснения обращайтесь! В качестве тестовой задачки картинка ниже:
31
Читателя, приславшего правильный ответ (зашифрованное слово) и несколько
слов о себе нижайше буду просить стать соавтором следующих брошюрок.
Как бы я хотел видеть в нашей компании древних изобретателей, о которых до
обидного мало написано в наших учебниках:
3.2.1. Ктесибий (годы деятельности 285–222 год до
н.э.) — древнегреческий изобретатель, математик и
механик,
живший
в
Александрии
в
Эллинистическом Египте. Ктесибия считают «отцом
пневматики». Он написал первые научные трактаты
об упругой силе сжатого воздуха и её использовании
в воздушных насосах и других механизмах (даже в
пневматическом
оружии),
заложил
основы
пневматики, гидравлики и теории упругости
воздуха. Ни одна из его письменных работ не
сохранилась, в том числе его Воспоминания.
Наиболее ясное описание пожарной машины
Ктесибия содержится в трактате его ученика Герона
Александрийского "Опыты с воздухом". Герону
приписывается очень важное усовершенствование, благодаря которому пожарная труба
стала приносить действительную пользу при тушении пожаров — изобретение
поворотной трубы (шейки).
3.2.2. Герон Александрийский, подробности его
жизни неизвестны, время жизни отнесено ко
второй половине I века н.э. Герон древнегреческий математик и инженер, создатель
множества
машин:
различных
сифонов,
механического театра марионеток, устройства для
открывания дверей и даже торговых автоматов,
продававших за мелкие монеты воду в египетских
храмах. Однако самое известное его изобретение –
"эолипил" или "шар Эола", представлявший собой
первую паровую турбину — шар, вращаемый
силой струй водяного пара...
32
3.3. Определение параметров методом Симою.
Начало идентификации систем, как этапа построения математических моделей
на основе наблюдений, связывают с работой Карла Фридриха Гаусса «Theoria motus
corporum coelestium in sectionibus conicis solem ambientium», в которой он использовал
разработанный им метод наименьших квадратов для предсказания траектории
движения планет.
Приблизительно до 50-х годов XX века, большая часть процедур идентификации
в автоматике, основывалась на наблюдении реакций управляемых объектов при
наличии некоторых управляющих воздействий... Проблема заключалась в том, что
область приложений этих методов была ограничена чаще всего скалярными
системами(SISO, Single-input,single-output). В 1960 году Рудольф Калман представил
описание управляемой системы в виде пространства состояний, что позволяло работать
и с многомерными (MIMO, Many-input,many-output) системами, и заложил основы для
оптимальной фильтрации и оптимального управления, основывавшихся на данном типе
описания. Конкретно для задач управления, методы идентификации систем были
разработаны в 1965 году в работах Хо и Калмана, Острёма и Болина. Эти работы
открыли путь разработке двух методов идентификации, популярных до сих пор: методу
подпространства и методу ошибки предсказания. Первый основан на использовании
проекций в евклидовом пространстве, а второй на минимизации критерия, зависящего
от параметров модели.
Общий подход к построению математической модели объекта требует три
базовые составляющие:
1. набор данных, полученных в результате нормальной работы изучаемого объекта
либо при целенаправленном эксперименте;
2. множество моделей-кандидатов для использования;
3. правило, по которому каждая модель-кандидат может быть принята или
отброшена.
С достаточно общих позиций математическое моделирование в книге [16]
рассматривается как один из методов познания реального мира в период
формирования так называемого информационного общества, как интеллектуальное
ядро быстро развивающихся информационных технологий. Этот метод не
противоречит хорошо известной формуле: "от живого созерцания к абстрактному
мышлению и от него к практике – таков диалектический путь познания истины,
познания объективной реальности" (В.И. Ленин).
Многообразие исследуемых объектов управления порождает многообразие
методов идентификации. По виду оператора систем математические модели систем
управления разделяют на следующие классы:
1. линейные и нелинейные;
2. непрерывные, дискретные, непрерывно-дискретные;
3. нестационарные и стационарные;
4. детерминированные и стохастические;
5. одномерные и многомерные;
6. системы с сосредоточенными и с распределенными параметрами.
33
Внешние воздействия подразделяют следующим образом:
1. непрерывные воздействия, которые описываются функциями непрерывного
аргумента, и дискретные воздействия, которые описываются соответственно,
функциями дискретной переменной;
2. детерминированные и случайные воздействия;
3. одномерные внешние воздействия, когда система имеет один входной сигнал, и
многомерные воздействия.
Таким образом, для определения типа конкретной системы (модели) надо указать
шесть классов, к которым относится оператор системы и три признака,
характеризующие внешнее воздействие. Например, характеристика некоторой системы
по классификационным признакам может быть такой:
линейная непрерывная нестационарная детерминированная одномерная система
с сосредоточенными параметрами при непрерывных детерминированных
многомерных внешних воздействиях.
При таком разнообразии типов систем считается хорошим тоном вспомнить
фундаментальную, правда, немного устаревшую книгу [17] профессора из Голландии
Питера Эйкхоффа, известного специалиста в области теории управления и одного из
создателей теории идентификации. В монографии содержатся необходимые сведения
из теории вероятностей, математической статистики, теории сигналов и динамических
систем. Систематически излагаются основные методы идентификации, оценивания
параметров и состояния систем управления. Теоретические положения
иллюстрируются многочисленными примерами. Отдельная глава посвящена
практическим применениям.
Построение модели по результатам наблюдений представляет собой
формализацию, необходимую для определения основных признаков, связей,
закономерностей, присущих объекту-оригиналу, и отсеивания второстепенных
признаков. Для одного и того же объекта в зависимости от конкретных требований
практики и типа решаемой задачи может быть построен ряд моделей, осуществлена
формализация различных функций этого объекта или внешних воздействий на этот
объект...
Кроме того различают два подхода к методам идентификации:
 методы идентификации вне контура регулирования;
 методы идентификации в замкнутом контуре;
и два способа реализации процедур:
 реализация, использующая явные математические выражения;
 реализация, использующая настраиваемые модели.
Метод площадей Симою М.П. позволяет определить передаточную функцию
модели объекта по кривой разгона. Кривая разгона – реакция динамического звена
(объекта регулирования) на скачкообразное воздействие произвольной амплитуды.
Кривая разгона может быть получена как экспериментально, так и расчетным путем.
Этот метод мгновенно получил широкое распространение (НИИ и ВУЗы) после
журнальной публикации [18], автором которой был Макс Паульсович Симою (1922 1973), сотрудник ЦНИИКА (Институт Комплексной автоматизации, Москва).
ЦНИИКА сотрудничает с многими зарубежными и российскими фирмами - среди
партнеров фирмы производители технических и программных средств автоматизации:
Siemens, Compressor Controls Corporation, Yokogawa, Motorola, Техноконт, RT-Soft и
другие... В мае 2012 года Федеральная Антимонопольная служба (ФАС) России
34
разрешила ООО "Балтнефтепровод" приобретение контрольного пакета (53%) акций
ЗАО ЦНИИКА... Необъяснимо!
Дальнейшее развитие метод получил в работах Евгения Павловича Стефани
(1915-1982),
например,
по
основам
расчета
настройки
регуляторов
теплоэнергетических процессов [19]. В книге автор преследовал цель дать в руки
инженера-практика, работающего в области автоматизации, систематизированное
изложение инженерных основ расчета систем автоматического регулирования. Пример
удачного использование методики, предложенной в публикациях [18] и [19], при
выполнении расчетов можно найти в курсовой работе студентов Уфимского
Государственного авиационного технического Университета Бескаравайного В.А. и
Хадиуллина С.Х., которая используется далее.
Для успешного синтеза систем управления сложными техническими объектами
необходимо иметь информацию об объекте управления. При синтезе линейных САУ
чаще всего используется описание в виде передаточных функций, для чего получают
переходную характеристику объекта и по ней строю аппроксимирующую
передаточную функцию, которую используют для разработки САУ. В данной работе
представлен алгоритм идентификации на основе метода Симою.
Для проверки алгоритма и программы был выбран объект с передаточной
функцией вида:
W(s ) 
K (b 1 s  1)
a 3 s 3  a 2 s  a1s  1
с параметрами:
K2
b1  1
a 1  1, a 2  0.33, a 3  0.036
т.е. моделируемый объект имеет передаточную функцию:
W (s ) 
2(s  1)
2(s  1)

,
2
0.036s  0.33s  s  1 (0.4s  1)(0.3s  1) 2
3
а полученная аналитически реакция на единичный скачек:
h0 (t )  48e

1
t
0.4
 50e

1
t
0.3
 46.66(6)  t  e

1
t
0.3
 2  1(t )
'
P R O G R A M "CA_SimoyusMethod"
'
14.08.2015
'------------------------------------------------------------'
Определение передаточной функции методом Симою
'------------------------------------------------------------' Курсовая работа: "Идентификация и синтез САУ"
' Алгоритм идентификации вариант 0
' Обозначено:
'
t - "физическое" время;
'
D - additive disturbance (аддитивноея возмущение);
'
DK - ключ аддитивного возмущения.
35
#Lang "FB"
' режим совместимости с FreeBASIC
' Общая область памяти
Common Shared h() As Double
' Объявление подпрограмма вывода графика на экран
Declare Sub ScrnOut(N As Integer, dt As Double)
'
' Основной модуль программы
Dim As Integer N, i, j, DK, hmax
N = 1600 ' переходный процесс dt*N = 8 секунд
Dim As Double h(N+1), D, sigm(N+1), sigm1, t, dt
Dim As Double I1(N+1), I2(N+1), I3(N+1), I4(N+1)
Dim As Double Fj(4), Ij(4)
Dim As Double K, b1, a1, a2, a3
Open Cons For Output As #2
'Open "CA_SimoyusMethod.res" For Output As #2
' параметры физического времени
For DK = 0 To 1
t = 0: dt = 0.005: hmax = 2
Print #2,: Print #2, "
Simoyus Method"
Print #2, Using " N = #### dt =##.#### DK =## "; N ; dt; DK
' данные тестового примера:
K = 2.0: b1 = 1.0: a1 = 1.0: a2 = 0.33: a3 = 0.036
' вывод точек графика теоретической функции
Print #2, "
t
h"
For i = 0 To N
t = i*dt
' D - возмущение (помеха), действующее на объект
If i > 0 And DK <> 0 Then D = 0.05*(Rnd() - 0.5) Else D = 0
h(i+1) = 48*exp(-1/0.4*t)-50*exp(-1/0.3*t)46.666666667*t*exp(-1/0.3*t)+2+D
If i Mod 100 = 0 Then Print #2, Using " ##.## ##.#######";
t; h(i+1)
Next i
' вычисление подинтергальных значений I1,I2,I3,I4
For i = 1 To N: sigm(i) = h(i)/hmax: Next i
t = 0
For i = 1 To N
t = (i-1)*dt
sigm1 =(1-sigm(i))
I1(i) =sigm1
I2(i) =sigm1*t
I3(i) =sigm1*t*t
I4(i) =sigm1*t*t*t
Next i
' вычисление непосредственно самих I1,I2,I3,I4
For j = 1 To 4: Ij(j) = 0: Next j
For i = 1 to N-1 ' интегрирование (середина прямоугольников)
Ij(1) = Ij(1)+((I1(i)+I1(i+1))/2)*dt
Ij(2) = Ij(2)+((I2(i)+I2(i+1))/2)*dt
Ij(3) = Ij(3)+((I3(i)+I3(i+1))/2)*dt
Ij(4) = Ij(4)+((I4(i)+I4(i+1))/2)*dt
Next i
36
' вычислене значений Fj
Fj(1) = Ij(1)
Fj(2) = Ij(1)*Ij(1)-Ij(2)
Fj(3) = Ij(1)*Ij(1)*Ij(1)-2*Ij(1)*Ij(2)+0.5*Ij(3)
Fj(4) =-Ij(4)/6+0.5*Ij(1)*Ij(3)+Fj(3)*Ij(1)-Fj(2)*Ij(2)
' вычисление параметров модели
b1 = -Fj(4)/Fj(3)
a1 = b1+Fj(1)
a2 = Fj(2)+b1*Fj(1)
a3 = Fj(3)+b1*Fj(2)
K = h(N+1)
Print #2, "
K
b1
a1
a2
a3"
Print #2, Using " ###.##### ###.##### ###.##### ###.#####
###.#####"; K; b1; a1; a2; a3
Next DK
Close #2
ScrnOut(N, dt)
Print: Print " Press to complete!!!"
Sleep
'
' Подпрограмма вывода графика на экран
Sub ScrnOut(N As Integer, dt As Double)
Dim As Integer i
Dim As Double t
Dim As Integer Wdt, Hgt
Wdt = 840: Hgt = 540 ' значения ширины и высоты
ScreenRes Wdt, Hgt
' вытянутый по ширине прямоугольник
' выделить область отсечения
View (20, 20) - (820, 520), 1, 15 ' часть экрана с синим фоном
' установить видимые координаты
Window (-1, -1) - (8, 3)
' "физическое" окно
' вычерчивание оси X
Line (-1,0)-(8,0), 7
Draw String(7.5, -0.1), "X"
' вычерчивание оси Y
Line (0,-1)-(0,3), 7
Draw String(-0.2, 2.8), "Y"
For i = 0 To N
' перебрать все точки
t = i*dt
PSet(t, h(i)), 13
' вывод розовой точки графика
Next i
End Sub
' Результат выполнения программы:
'
'
Simoyus Method
' N = 1600 dt = 0.0050 DK = 0
'
t
h
' 0.00
0.0000000
' 0.50
1.9013527
' 1.00
2.4915939
' 1.50
2.3202982
37
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
'
2.00
2.50
3.00
3.50
4.00
4.50
5.00
5.50
6.00
6.50
7.00
7.50
8.00
2.1410106
2.0526002
2.0179221
2.0057768
2.0017959
2.0005448
2.0001625
2.0000479
2.0000140
2.0000041
2.0000012
2.0000003
2.0000001
K
2.00000
b1
1.00000
Simoyus Method
N = 1600 dt = 0.0050
t
h
0.00
0.0000000
0.50
1.9227797
1.00
2.4964339
1.50
2.3348922
2.00
2.1579559
2.50
2.0373975
3.00
2.0032391
3.50
1.9935189
4.00
2.0057459
4.50
1.9912256
5.00
1.9857328
5.50
2.0250444
6.00
2.0163569
6.50
1.9828801
7.00
1.9982803
7.50
1.9762452
8.00
2.0111176
K
b1
2.01112
1.00947
a1
1.00000
a2
0.33000
a3
0.03600
a2
0.32879
a3
0.03789
DK = 1
a1
1.00933
Здесь используется подпрограмма ScreenRes, которая инициализирует
графический режим с указанными горизонтальным и вертикальным разрешениями. Он
уже была упомянута в брошюре "Фрагмент 18. Опустите ваши кисти..." Конечно,
программа должна быть существенно облагорожена и комментирована, но главное –
результат вполне обнадеживающий.
Дерзайте!!!
Как говорил древнегреческий философ Сократ (469 до н.э. — 399 до н.э.): "Нет
предела совершенству"! Ему вторит испанский художник и писатель Сальвадор Дали
(1904 —1989): "Не бойтесь совершенства. Вам его не достичь"!
А я буду закругляться...
38
Чтобы всегда под рукой была подходящая литература, завершая "лапидарий"
могу указать еще трехтомный "свежий" учебник под редакцией Н.Д.Егупова (МГТУ
им.Н.Э.Баумана), посвященный методы классической и современной теории
автоматического управления:






Том 1 [20] Анализ и статистическая динамика систем автоматического
управления.
В первом томе учебника изложены основные положения классической теории
автоматического управления: основные понятия и принципы управления,
методы математического описания стационарных, нестационарных и
нелинейных непрерывных систем и исследования их устойчивости и качества
процессов управления; подробно рассмотрен метод пространства состояний.
Том 2 [21] Синтез
регуляторов
и
теория
оптимизации
систем
автоматического управления.
Во втором томе учебника изложены методы синтеза регуляторов,
обеспечивающих заданное качество процессов управления и позволяющих
определить состав, структуру САУ и параметры всех ее устройств из условия
удовлетворения заданному комплексу технических требований в классе
линейных (стационарных и нестационарных), нелинейных, дискретных и
многомерных систем Отражены основные положения модального управления.
Изложены основы теории оптимальных систем.
Том 3 [22] Методы современной теории автоматического управления.
В третьем томе учебника изложены основные теоретические положения
некоторых направлений теории автоматического управления, интенсивно
развиваемых в последние десятилетия. Рассмотрены методы расчета и
проектирования систем, использующие аппарат дифференциальной геометрии.
Отражены центральные положения теории катастроф, теории хаоса; приведены
понятия, связанные с теорией фракталов и их использованием. Рассмотрена
теория робастного управления.
39
Список литературы.
[01] Даль В.И. Пословицы русского народа. Сборник. CПб.-М.: Издание
книгопродавца-типографа М.О.Вольфа. 1879. 750+642 с.
[02] Ватсон Г.Н. Теория Бесселевых функций. Часть 1. М.: ИнЛит, 1949. 800 с.
[03] Справочник по специальным функциям, с формулами, графиками и
математическими таблицами. Абрамовиц М., Стиган И. М.: 1979. 832 с.
[04] Зубов И.В. Функции Бесселя. М.: МФТИ, 2007. 51 с.
[05] Библиотека алгоритмов 1б-50б. Справочное пособие. Вып. 1. / М.И.Агеев,
В.П.Алик, Ю.И.Марков, Р.М.Галис. М.: Советское радио, 1975. 175 с.
[06] Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль.
Томск: МП "РАСКО", 1991. 272 с.
[07] Егоров А.И. Обыкновенные дифференциальные уравнения с приложениями. М.:
ФизМатЛит, 2005. 384 с.
[08] Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных
уравнений. Нежесткие задачи. М.: Мир, 1990. 512 с.
[09] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений.
Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.
[10] Красовский А.А. Справочник по теории автоматического управления. М.: Наука,
1987. 712 с.
[11] Афанасьев В.Н., Колмановский В.Б., Носов В.Р. Математическая теория
конструирования систем управления. М.: Высшая школа, 2003. 615 с.
[12] Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.:
ФизМатГиз, 1960. 655 с.
[13] Крылов В.И., Бобков В.В., Монастырный П.И. Начала теории вычислительных
методов. Линейная алгебра и нелинейные уравнения. Мн.: Наука и техника, 1985. 280 с.
[14] Деммель Дж. Вычислительная линейная алгебра. Теория и приложения. М.: Мир.
2001. 430 с.
[15] Уилкинсон Дж.Х., Райнш С. Справочник алгоритмов на языке Алгол. Линейная
алгебра. М.: Машиностроение, 1976. 390 с.
[16] Математическое моделирование в технике. Учеб. для вузов / Под ред. B.C.
Зарубина, А.П. Крищенко. М.: Изд-во МГТУ им. Н.Э. Баумана, 2003. 496 с.
40
[17] Эйкхофф П. Основы идентификации систем управления. М.: Мир, 1975. 680 с.
[18] Симою М.П. Определение коэффициентов передаточных функций
линеаризованных звеньев систем регулирования. Автоматика и телемеханика, 1957. №
6, с.514–527.
[19] Стефани Е.П. Основы расчета настройки регуляторов теплоэнергетических
процессов. М.: Энергия, 1972. 376 с.
[20] Методы классической и современной теории автоматического управления. Том 1:
Анализ и статистическая динамика систем автоматического управления/Под ред.
Н.Д.Егупова. М.: Изд-во МГТУ им.Н.Э.Баумана, 2000, 748 с.
[21] Методы классической и современной теории автоматического управления. Том 2:
Синтез регуляторов и теория оптимизации систем автоматического управления/Под
ред. Н.Д.Егупова. М.: Изд-во МГТУ им.Н.Э.Баумана, 2000, 736 с.
[22] Методы классической и современной теории автоматического управления. Том 3:
Методы современной теории автоматического управления/Под ред. Н.Д.Егупова. М.:
Изд-во МГТУ им.Н.Э.Баумана, 2000, 748 с.
Простите, что изложение получилось такое сумбурное –
был жестко ограничен по времени.
Кроме того были затруднения с размещением файла в хранилище,
что привело к скачку порядкового номера. Простите!
До новых встреч!
Жду соавтора!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
47
Размер файла
1 405 Кб
Теги
Симою, алгоритм, хессенберг, бессель, freebasic
1/--страниц
Пожаловаться на содержимое документа