Статистический анализ и оптимизация САР Привод сопла ракеты носителя
Московский
Государственный Авиационный Институт
(Технический
университет)
Кафедра 704
Информационно-управляющие комплексы
Êóðñîâàÿ
ðàáîòà
Òåìà: “Ñòàòèñòè÷åñêèé
àíàëèç ÑÀД
Âûïîëíèë: ñò. ãð. 07-403
Êîðíèëîâ
Ä.Ì.
Ðóêîâîäèòåëü: Êóäðÿøîâ
Ñ. Â.
Ìîñêâà
1998ã.
1.
Задание
Ïðèâîä
ãèáêîãî
ñîïëà
ðàêåòû-íîñèòåëÿ:
Äàííûå:
Ïàðàìåòðû
âîçäåéñòâèé
íà âõîäàõ
ñèñòåìû
çàäàíû â
âèäå
êîððåëÿöèîííûõ
ôóíêöèé:
ãäå:
l=1.2
Óñèëèòåëüíîå
çâåíî:
Kp=13 1/c
Ïàðàìåòðû
ïåðâîé
íåëèíåéíîñòè:
S=0.5;
K1=2.5
K2=1.9
Ïàðàìåòðû
âòîðîé
íåëèíåéíîñòè:
S=27 °/c
Ïàðàìåòðû
òðåòüåé
íåëèíåéíîñòè:
S=3 °
1. Ïðîìîäåëèðîâàòü
ñîñòîÿíèå
ñèñòåìû â
íåëèíåéíîì
âèäå.
2.
Ëèíåàðèçîâàòü
íåëèíåéíûå
ýëåìåíòû è
ïðîìîäåëèðîâàòü
ñîñòîÿíèå
ñèñòåìû â
ëèíåàðèçîâàííîì
âèäå.
3.
Ïîñòðîèòü
ýâîëþöèþ
ìàòðèöû
êîâàðèàöèé
Содержание
1.
Задание.............................................................................................................................................................
2. Теоретическая часть............................................................................................................................
2.1
Случайные процессы и их математическое описание............................................
2.2
Прохождение стационарного процесса через линейную динамическую систему
2.3
Формирующий фильтр......................................................................................................................
2.4
Априорный статистический анализ......................................................................................
2.5
Статистическая линеаризация................................................................................................
3. Реализация....................................................................................................................................................
3.1
Система дифференциальных уравнений............................................................................
3.2
Расчет системы в нелинейной форме....................................................................................
3.3
Расчет линеаризованной системы.........................................................................................
4. Заключение..................................................................................................................................................
5. Список литературы................................................................................................................................
2.
Теоретическая
часть
2.1 Случайные процессы и их
математическое описание
Ïóñòü
t ïðèíàäëåæèò
T (äîïóñòèìîìó
ìíîæåñòâó).
Åñëè t
ïðîáåãàåò
íåïðåðûâíûå
çíà÷åíèÿ íà
ìíîæåñòâå T, òî x(t) ïðèíÿòî
íàçûâàòü ñëó÷àéíûì
ïðîöåññîì.
Ïðè
êàæäîì
ôèêñèðîâàííîì
t=t* âîçíèêàåò
ñëó÷àéíàÿ
âåëè÷èíà x(t*)
êîòîðóþ
ïðèíÿòî
íàçûâàòü çíà÷åíèåì
ñëó÷àéíîãî
ïðîöåññà.
Ñëó÷àéíûé
ïðîöåññ
õàðàêòåðèçóåòñÿ
ñîâîêóïíîñòüþ
ïëîòíîñòåé
ðàñïðåäåëåíèÿ
âåðîÿòíîñòåé
ñ
âîçðàñòàþùåé
ðàçìåðíîñòüþ
k=1,2,...,n.
Äåéñòâèòåëüíî
âåëè÷èíà
ðàâíà
âåðîÿòíîñòè
òîãî, ÷òî
Ïîýòîìó
÷åì áîëüøå n, òåì
áîëåå ïîëíîé
èíôîðìàöèåé
î ïîâåäåíèè x(t) â
èíòåðåñóþùåì
íàñ
èíòåðâàëå
âðåìåíè ìû
ðàñïîëàãàåì.
Ïðàêòè÷åñêè
îãðàíè÷èâàþòñÿ
ðàññìîòðåíèåì
òîëüêî
îäíîìåðíûõ
è äâóìåðíûõ
ïëîòíîñòåé
ðàñïðåäåëåíèÿ
ëèáî èíûõ õàðàêòåðèñòèê
ñëó÷àéíûõ
ïðîöåññîâ
(ãëàâíûì
îáðàçîì
ìîìåíòîâ
ïåðâîãî è
âòîðîãî ïîðÿäêîâ),
êîòîðûå
îïðåäåëÿþòñÿ
äàííûìè ïëîòíîñòÿìè.
Ïðèìåðîì
ñëó÷àéíîãî
ïðîöåññà,
ïîëíîñòüþ õàðàêòåðèçóåìîãî
îäíîìåðíîé
è äâóìåðíîé
ïëîòíîñòÿìè,
ÿâëÿåòñÿ ìàðêîâñêèé
ñëó÷àéíûé
ïðîöåññ.
Çàâèñèìîñòü
ìåæäó
çíà÷åíèÿì x(ti)
ÿâëÿåòñÿ
ïðîñòåéøåé,
òàê êàê
ðàñïðîñòðàíÿåòñÿ
ëèøü íà
ñîñåäíèå
çíà÷åíèÿ x(ti-1) è x(ti).
Íàëè÷èå
ïîäîáíîé
çàâèñèìîñòè
ïðèâîäèò ê
òîìó, ÷òî
âåðîÿòíîñòü
íàõîæäåíèÿ x(ti) â
èíòåðâàëå [xi, xi+dxi]
â ìîìåíò
âðåìåíè t=ti ÿâëÿåòñÿ
óñëîâíîé è
çàâèñèò îò
çíà÷åíèÿ
ñëó÷àéíîãî
ïðîöåññà â
ïðåäûäóùèé
ìîìåíò
âðåìåíè ti-1.Çàâèñèìîñòü
x(ti)
îò áîëåå
ðàííèõ
ìîìåíòîâ
âðåìåíè t1, t2,
ti-2, (ò. å. îò
áîëåå
ãëóáîêîé
ïðåäûñòîðèè
ïðîöåññà)
îòñóòñòâóåò.
Ýòî
îçíà÷àåò,
÷òî äëÿ
ìàðêîâñêîãî
ïðîöåññà
óñëîâíàÿ
(èëè
ïåðåõîäíàÿ)
ïëîòíîñòü:
Îòñþäà:
Òàêèì
îáðàçîì,
íà÷àëüíàÿ
áåçóñëîâíàÿ
îäíîìåðíàÿ
ïëîòíîñòü è
ñîâîêóïíîñòü
óñëîâíûõ
(ïåðåõîäíûõ)
ïëîòíîñòåé
ïîëíîñòüþ
îïèñûâàþò
ìàðêîâñêèé
ñëó÷àéíûé
ïðîöåññ.
Àáñîëþòíî
ñëó÷àéíûì
ïðîöåññîì ïðèíÿòî
íàçûâàòü
òàêîé
ïðîöåññ,
ëþáûå äâà
çíà÷åíèÿ
êîòîðîãî
ñóòü
íåçàâèñèìûå
ñëó÷àéíûå
âåëè÷èíû. Â
ýòîì ñëó÷àå
ïëîòíîñòü
âåðîÿòíîñòè
èìååò
ñëåäóþùèé
âèä:
Ñëó÷àéíûé
ïðîöåññ
íàçûâàåòñÿ ñòàöèîíàðíûì,
åñëè âñå åãî
ïëîòíîñòè
âåðîÿòíîñòåé
íå çàâèñÿò
îò âûáîðà
íà÷àëà
îòñ÷åòà
âðåìåíè, ò. å.
èíâàðèàíòíû
ê
âðåìåííîìó
ñäâèãó t:
Èç
ýòîãî
ñëåäóåò, ÷òî
îäíîìåðíàÿ
ïëîòíîñòü
ðàñïðåäåëåíèÿ
ñòàöèîíàðíîãî
ïðîöåññà
âîîáùå íå
çàâèñèò îò
âðåìåíè.
Ãàóññîâñêèé
ïðîöåññ -ýòî
òàêîé
ñëó÷àéíûé
ïðîöåññ
ñêîëü óãîäíî
ìåðíàÿ
ïëîòíîñòü
âåðîÿòíîñòè
êîòîðîãî ãàóññîâñêàÿ.
n-ðàçìåðíîñòü
âåêòîðà X,
Kx-ìàòðèöà
êîâàðèàöèè
mx-ìàòåìàòè÷åñêîå
îæèëàíèå.
Ãàóññîâñêèé
ñëó÷àéíûé
ïðîöåññ
ÿâëÿåòñÿ
ñòàöèîíàðíûì
è
ìàðêîâñêèì.
Ê
íàèáîëåå
âàæíûì
ìîìåíòíûì
õàðàêòåðèñòèêàì
ñòàöèîíàðíîãî
ñëó÷àéíîãî
ïðîöåññà
îòíîñÿòñÿ
ìàòåìàòè÷åñêîå
îæèäàíèå,
äèñïåðñèÿ,
êîððåëÿöèîííàÿ
ôóíêöèÿ.
Ìàòåìàòè÷åñêîå
îæèäàíèå
õàðàêòåðèçóåò
ñðåäíåå
òå÷åíèå
ïðîöåññà x(t) ïî
âðåìåíè.
Äèñïåðñèÿ
ñëó÷àéíîãî
ïðîöåññà
Êîððåëÿöèîííàÿ
ôóíêöèÿ
ãäå , ;ýòà
ôóíêöèÿ
ïðåäñòàâëÿåò
ñîáîé
ñðåäíåå ïðîèçâåäåíèå
öåíòðèðîâàííûõ
çíà÷åíèé ñëó÷àéíîãî
ïðîöåññà â
ìîìåíòû
âðåìåíè t è t+t.
Êîððåëÿöèîííàÿ
ôóíêöèÿ
õàðàêòåðèçóåò
ñòåïåíü
ëèíåéíîé
ñâÿçè
(êîððåëÿöèè)
ìåæäó
çíà÷åíèÿìè
ïðîöåññà,
îòñòîÿùèìè
äðóã îò äðóãà
íà âðåìÿ t.
Ïðè t=0, êîððåëÿöèîííàÿ
ôóíêöèÿ
ðàâíà
äèñïåðñèè.
Ïîíÿòèå
êîððåëÿöèîííîé
ôóíêöèè
ìîæåò áûòü
èñïîëüçîâàíî
è äëÿ õàðàêòåðèñòèêè
ñòåïåíè
ñâÿçè äâóõ
ñëó÷àéíûõ
ïðîöåññîâ x(t) è y(t). Â
ýòîì ñëó÷àå
îíà
íàçûâàåòñÿ
âçàèìíîé êîððåëÿöèîííîé
ôóíêöèåé:
Â
òåîðèè
àâòîìàòè÷åñêîãî
óïðàâëåíèÿ
øèðîêî
èñïîëüçóþòñÿ
îïèñàíèÿ
ñëó÷àéíûõ ïðîöåññîâ
â ÷àñòîòíîé
îáëàñòè èëè,
ïî èíîìó,
ñïåêòðàëüíîå
ïðåäñòàâëåíèå
ñëó÷àéíûõ
ïðîöåññîâ.
Ðàññìîòðèì
ïðåîáðàçîâàíèå
Ôóðüå îò êîððåëÿöèîííîé
ôóíêöèè:
Ïîëó÷åííàÿ
ôóíêöèÿ
åñòü ÷åòíàÿ
âåùåñòâåííàÿ
ôóíêöèÿ
íàçûâàåìàÿ
ñïåêòðàëüíîé
ïëîòíîñòüþ
ñòàöèîíàðíîãî
ñëó÷àéíîãî
ïðîöåññà.
Ñïðàâåäëèâî
îáðàòíîå:
Äëÿ
ñòàöèîíàðíîãî
ñëó÷àéíîãî
ïðîöåññà n(t) (ñ
íóëåâûì
ìàòåìàòè÷åñêèì
îæèäàíèåì) òèïà
áåëîãî øóìà,
êîððåëÿöèîííàÿ
ôóíêöèÿ
èìååò âèä:
ãäå
d(t)-äåëüòà-ôóíêöèÿ
Äèðàêà, à N
-èíòåíñèâíîñòü
øóìà.
Ñïåêòðàëüíàÿ
ïëîòíîñòü
ýòîãî
ïðîöåññà áóäåò:
÷òî
ìîæåò áûòü
ïðèíÿòî â
êà÷åñòâå
îïðåäåëåíèÿ
áåëîãî øóìà.
Âûðàæåíèå
îçíà÷àåò,
÷òî ìîùíîñòü
ïàðöèàëüíûõ
ñîñòàâëÿþùèõ
ñëó÷àéíîãî
ïðîöåññà n(t) äëÿ ëþáûõ
÷àñòîò
îäíà è òà æå.
Ïîýòîìó
áåëûé øóì
ÿâëÿåòñÿ
íàèáîëåå
èíòåíñèâíûì
âèäîì ïîìåõè.
2.2 Прохождение стационарного
процесса через линейную динамическую систему
Ðàññìîòðèì
ëèíåéíóþ
äèíàìè÷åñêóþ
ñèñòåìó ñ
ïîñòîÿííûìè
êîýôôèöèåíòàìè.
Íà åå âõîä
ïîñòóïàåò
ñòàöèîíàðíûé
ñëó÷àéíûé
ïðîöåññ x(t), íà
âûõîäå
èìååò ìåñòî
ïðîöåññ y(t).
Òåîðåòè÷åñêè
âûõîäíîé ñëó÷àéíûé
ïðîöåññ y(t) ÿâëÿåòñÿ
ñòàöèîíàðíûì
òîëüêî ïîñëå
çàòóõàíèÿ
ñâîáîäíûõ
êîëåáàíèé â
ñèñòåìå, òî åñòü
ïðè
t®¥.
Îäíàêî, â
èíæåíåðíûõ
ïðèëîæåíèÿõ
ìû áóäåì
ñ÷èòàòü, ÷òî
ïåðåõîäíûé
ïðîöåññ â
ñèñòåìå
çàêàí÷èâàåòñÿ
çà âðåìÿ,
îïðåäåëÿåìîå
â
ñîîòâåòñòâèå
ñ ïðàâèëàìè
òåîðèè óïðàâëåíèÿ.
Èíûìè
ñëîâàìè ìû
áóäåì
ñ÷èòàòü
ïðîöåññ y(t)
ñòàöèîíàðíûì
ïî
èñòå÷åíèè
âðåìåíè çàòóõàíèÿ.
Ñïåêòðàëüíàÿ
ïëîòíîñòü
âûõîäíîãî
ïðîöåññà
èìååò âèä:
Äèñïåðñèÿ
âûõîäíîãî
ïðîöåññà:
Êîððåëÿöèîííàÿ
ôóíêöèÿ
âûõîäíîãî
ïðîöåññà:
2.3 Формирующий фильтр
Êàê
áûëî
ïîêàçàíî
âûøå áåëûé
øóì èìååò ïîñòîÿííóþ
ñïåêòðàëüíóþ
ïëîòíîñòü âî
âñåì
äèàïàçîíå
÷àñòîò.
Ñïåêòðàëüíàÿ
ïëîòíîñòü
âñåõ
ôèçè÷åñêè ñóùåñòâóþùèõ
ñòàöèîíàðíûõ
ñëó÷àéíûõ ïðîöåññîâ
ïðåäñòàâëÿþò
ñîáîé
äðîáíî-ðàöèîíàëüíóþ
ôóíêöèþ
÷àñòîòû:
Ïðè÷åì
ñòåïåíü
ïîëèíîìà â
çíàìåíàòåëå
âûøå
ñòåïåíè
ïîëèíîìà â
÷èñëèòåëå.
Òàêèå
ôóíêöèè
äîïóñêàþò
ôàêòîðèçàöèþ:
,
ãäå
êâàäðàò
ìîäóëÿ
àìïëèòóäíîé
õàðàêòåðèñòèêè
íåêîåé
ôèêòèâíîé
ìèíèìàëüíî-ôàçîâîé
äèíàìè÷åñêîé
ñèñòåìû,
êîòîðóþ â
äàëüíåéøåì
ìû áóäåì íàçûâàòü
ôîðìèðóþùåì
ôèëüòðîì ñîîòâåòñòâóþùèì
ñïåêòðàëüíîé
ïëîòíîñòè íåêîåãî
ñëó÷àéíîãî
ïðîöåññà.
2.4 Априорный статистический
анализ
Ïîä
àïðèîðíûì
ñòàòèñòè÷åñêèì
àíàëèçîì
(èëè
àíàëèçîì
òî÷íîñòè)
ïîíèìàåòñÿ
îïðåäåëåíèå
ñòàòèñòè÷åñêèõ
õàðàêòåðèñòèê
(ìàòåìàòè÷åñêèõ
îæèäàíèé,
äèñïåðñèé,
ñïåêòðàëüíûõ
ïëîòíîñòåé,
ðàñïðåäåëåíèå
âåðîÿòíîñòåé
è ò. ï.)
êîîðäèíàò
óïðàâëÿåìîãî
äèíàìè÷åñêîãî
îáúåêòà ïî
èçâåñòíîìó
åãî äèôôåðåíöèàëüíîìó
óðàâíåíèþ
äâèæåíèÿ è
ñòàòèñòè÷åñêèì
õàðàêòåðèñòèêàì
ñëó÷àéíûõ
ôàêòîðîâ.
Ïóñòü
ëèíåàðèçîâàííûå
óðàâíåíèÿ
âîçìóùåííîãî
äâèæåíèÿ
óïðàâëÿåìîãî
îáúåêòà èìåþò
âèä:
ãäå
x(t)-âåêòîð
ñîñòîÿíèÿ
(ôàçîâûé âåêòîð),
ðàçìåðíîñòè
nx1,
A(t)-ìàòðèöà
êîýôôèöèåíòîâ,
ðàçìåðíîñòè
nxn.
B(t)-
ìàòðèöà
êîýôôèöèåíòîâ
,áåëûõ øóìîâ,
ðàçìåðíîñòè nxm.
n-âåêòîð
áåëûõ øóìîâ,
ðàçìåðíîñòè
mx1.
Òîãäà
äèôôåðåíöèàëüíûå
óðàâíåíèÿ
äëÿ âåêòîðà
ìàòåìàòè÷åñêèõ
îæèäàíèé è
ìàòðèöû
êîâàðèàöèé
èìåþò ñëåäóþùèé
âèä:
Ðàçìåðíîñòü
ìàòðèöû
êîâàðèàöèè nxn.
N-äèàãîíàëüíàÿ
ìàòðèöà
èíòåíñèâíîñòåé áåëûõ
øóìîâ.
Äèôôåðåíöèàëüíûå
óðàâíåíèÿ (1)-(3)
ðåøàþòñÿ
îäíèì èç
÷èñëåííûõ
ìåòîäîâ
èíòåãðèðîâàíèÿ.
Òàêèì
îáðàçîì, ìû
îïðåäåëÿåì
âåêòîð
ñîñòîÿíèÿ è
ñòàòèñòè÷åñêèå
õàðàêòåðèñòèêè
ñèñòåìû â
ëþáîé
ìîìåíò
âðåìåíè.
Ïåðåä
íà÷àëîì
èíòåãðèðîâàíèÿ,
äîëæíû áûòü
èçâåñòíû
àïðèîðíûå
çíà÷åíèÿ
âåêòîðà
ñîñòîÿíèÿ,
âåêòîðà ìàòåìàòè÷åñêèõ
îæèäàíèé è
ìàòðèöû
êîâàðèàöèé
â íà÷àëüíûé
ìîìåíò
âðåìåíè.
2.5 Статистическая линеаризация
Ëåãêî
âèäåòü, ÷òî
äëÿ ðåøåíèÿ
óðàâíåíèé èç
ïóíêòà 2.4
íåîáõîäèìû
ëèíåéíûå
ñèñòåìû óðàâíåíèé.
Îäíàêî íà
ïðàêòèêå
ñèñòåìû
óïðàâëåíèÿ
ìîãóò ñîäåðæàòü
(è ÷àùå âñåãî
ñîäåðæàò) íåëèíåéíûå
ýëåìåíòû, è
óðàâíåíèå
äëÿ âåêòîðà ñîñòîÿíèé
ïðèíèìàåò
âèä:
Â
ýòîì ñëó÷àå
ïðèìåíÿåòñÿ
ìåòîä
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè,
êîãäà
íåëèíåéíûé ýëåìåíò
çàìåíÿåòñÿ
ëèíåéíûì â
íåêîòîðîì
ñìûñëå
ýêâèâàëåíòíûì.
Ïóñòü
íåëèíåéíûé
ýëåìåíò
èìååò
ñëåäóþùèé
âèä:
Ââåäåì
ëèíåéíûé
ýëåìåíò
ñëåäóþùåãî
âèäà:
,
ãäå
Íåîáõîäèìî
÷òîáû
âåëè÷èíà íà
âûõîäå ëèíåéíîãî
ýëåìåíòà
áûëà
ýêâèâàëåíòíà,
â íåêîòîðîì
ñìûñëå,
âåëè÷èíå íà
âûõîäå íåëèíåéíîãî
ýëåìåíòà.
Ñóùåñòâóþò
äâà ïîäõîäà:
1.
Êðèòåðèé
âèäà:
M{z}=M{h}
D{z}=D{h}
Ôîðìóëû
äëÿ
êîýôôèöèåíòîâ
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè:
2.
Âòîðîé
ñïîñîá
çàêëþ÷àåòñÿ
â
âûïîëíåíèè êðèòåðèÿ
âèäà:
M{z}=M{h}
D{h-z}®min
Êîýôôèöèåíò
b
âû÷èñëÿåòñÿ
ïî ôîðìóëå
3.
Реализация
Äëÿ
ðåøåíèÿ
ïîñòàâëåííîé
çàäà÷è áûëî
íàïèñàíî
ïðîãðàììíîå
îáåñïå÷åíèå
ñ ïîìîùüþ
ñðåäû Microsoft Visual C++ 4.0 äëÿ
ìàòðè÷íûõ
îïåðàöèé,
÷èñëåííûõ
ìåòîäîâ
èíòåãðèðîâàíèÿ.
Îñíîâíàÿ
çàäà÷à ðåøàåòñÿ
â äâóõ
ïðîãðàììàõ,
äëÿ ðàñ÷åòà
íåëèíåéíîé
ñèñòåìû è
ëèíåàðèçîâàííîé.
3.1 Система дифференциальных
уравнений
Äëÿ
òîãî ÷òîáû
ââåñòè â
ñèñòåìó
ñëó÷àéíûå
âîçìóùåíèÿ
ñ
òðåáóåìûìè
êîððåëÿöèîííûìè
ôóíêöèÿìè
âîñïîëüçóåìñÿ
ïîíÿòèåì
ôîðìèðóþùåãî
ôèëüòðà,
äèíàìè÷åñêîãî
çâåíà íà
âõîä êîòîðîãî
ïîñòóïàåò
áåëûé øóì, à
íà âûõîäå ïðîöåññ
ñ
òðåáóåìûìè
ïàðàìåòðàìè.
Èòàê,
ñïåêòðàëüíàÿ
ïëîòíîñòü
òðåáóåìîãî
ïðîöåññà
èìååò âèä:
Ñîãëàñíî
ôîðìóëå ïåðåäàòî÷íàÿ
ôóíêöèÿ
ôîðìèðóþùåãî
ôèëüòðà èìååò
âèä:
Ïîñëå
íåñëîæíûõ
ïðåîáðàõîâàíèé
ïîëó÷àåì
äèôôåðåíöèàëüíûå
óðàâíåíèÿ
äëÿ ôîðìèðóþùèõ
ôèëüòðîâ
äâóõ âõîäîâ.
Òðåòüåé
ôàçîâîé
êîîðäèíàòîé
áóäåò x,
çíà÷åíèå íà
âõîäå â
òðåòüþ
íåëèíåéíîñòü.
Åãî
óðàâíåíèå
èìååò âèä:
Ýòè
óðàâíåíèÿ è
ñîñòàâÿò
ñèñòåìó
äèôôåðåíöèàëüíûõ
óðàâíåíèé:
3.2
Расчет
системы в нелинейной форме
Äëÿ
ðàñ÷åòà
äàííîé
ñèñòåìû â
íåëèíåéíîì
âèäå áûëà
ðàçðàáîòàíà
ñîîòâåòñòâóþùàÿ
ïðîãðàììà
êîòîðàÿ
èíòåãðèðîâàëà
ñèñòåìó
óðàâíåíèÿ
ìåòîäîì Ðóíãå-Êóòòà
÷åòâåðòîãî
ïîðÿäêà.
Äëÿ
ðåøåíèå
äèôôåðåíöèàëüíîãî
óðàâíåíèÿ
ìåòîäîì
Ðóíãå-Êóòòà
íåîáõîäèìî
çíàòü íà÷àëüíûå
óñëîâèÿ, òî
åñòü
çíà÷åíèÿ
âåêòîðà
ñîñòîÿíèé â
íóëåâîé
ìîìåíò
âðåìåíè, â
äàííîì
ñëó÷àå îíè
íóëåâûå.
Ïîñëå
îñóùåñòâëåíèÿ
èíòåãðèðîâàíèÿ
ðåçóëüòàòû
áûëè
çàïèñàíû â
ôàéë, è çàòåì
ïîñòðîåíû.
Ãðàôèêè
èìåþò
ñëåäóþùèé
âèä:
Ðèñóíîê
1.
3.3
Расчет
линеаризованной системы
Äëÿ
ðàñ÷åòà
ñèñòåìû â
ëèíåàðèçîâàííûì
âèäå ïî
ôîðìóëàì,
íåîáõîäèìî
ëèíåàðèçîâàòü
ñèñòåìó
äèôôåðåíöèàëüíûõ
óðàâíåíèé, ò.
å. çàìåíèòü
òðè
íåëèíåéíûõ ýëåìåíòà
ëèíåéíûìè
êàê ýòî
ïîêàçàíî â 2.5. Äëÿ
òîãî ÷òîáû
ïîëó÷èòü
êîýôôèöèåíòû
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè
íàì íåîáõîäèìî
çíàòü
ïàðàìåòðû
ñëó÷àéíîãî
ïðîöåññà íà
âõîäå â
íåëèíåéíîñòü
(ìàòåìàòè÷åñêîå
îæèäàíèå è
äèñïåðñèþ).
Ïàðàìåòð
âõîäà â
òðåòüþ
íåëèíåéíîñòü
èìååòñÿ â
âåêòîðå
ñîñòîÿíèé.
Ïîýòîìó åãî
äèñïåðñèÿ
ïðèñóòñòâóåò
â ìàòðèöå
êîâàðèàöèè.
Äëÿ
íàõîæäåíèÿ
äèñïåðñèé
íà âõîäå â
ïåðâóþ è
âòîðóþ
íåëèíåéíîñòè
ïðîèçâåäåì
ñëåäóþùèå
äåéñòâèÿ.
Ïðåîáðàçóåì
ñèñòåìó
òàêèì îáðàçîì
÷òîáû
âõîäîì
ñèñòåìû
îñòàâàëñÿ
îäèí èç
âõîäîâ â
íàøó
ñèñòåìó, à
âûõîäîì âõîä
â
íåëèíåéíîñòü
(äâà âõîäà
ðàññìàòðèâàþòñÿ
îòäåëüíî, à
çàòåì èõ
äèñïåðñèè
ñêëàäûâàþòñÿ).
Èòàê,
ïîñ÷èòàåì
äèñïåðñèè
îò
âîçäåéñòâèé
a
è h,
íà âõîäàõ
ïåðâîé è
âòîðîé
íåëèíåéíîñòè:
Íàçîâåì
ñèãíàë íà
âõîäå â
ïåðâóþ
íåëèíåéíîñòü
õ,
âî âòîðóþ y òîãäà
óðàâíåíèå
ñâÿçûâàþùèå
x è a
áóäóò èìåòü
âèä:
Çäåñü
êîýôôèöèåíòû
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè
ñîîòâåòñòâåííî
ïåðâîé è âòîðîé
íåëèíåéíîñòåé, äàëåå:
Ñëåäîâàòåëüíî
ïåðåäàòî÷íàÿ
ôóíêöèÿ
áóäåò èìåòü
âèä:
Ñïåêòðàëüíàÿ
ïëîòíîñòü
ñëó÷àéíîãî
ïðîöåññà íà
âõîäå a
èìååò âèä:
Òîãäà
äèñïåðñèÿ íà
âõîäå ïåðâîé
íåëèíåéíîñòè
áóäåò èìåòü
âèä:
Â
òåîðèè
èçâåñòíî,
÷òî
èíòåãðàëû
âèäà:
ãäå:
Â
àíàëèòè÷åñêîé
ôîðìå èìåþò
âèä:
, ãäå
Ðàññìîòðèì
îòäåëüíî
çíàìåíàòåëü
íàøåãî
èíòåãðàëà,
ïðèâåäåííûé
ê âèäó:
Òîãäà:
Îòñþäà
ñëåäóåò:
Îòêóäà:
Ðàññóæäàÿ
àíàëîãè÷íûì
îáðàçîì
ïîëó÷èì îñòàëüíûå
äèñïåðñèè:
Îòñþäà
ñëåäóåò, ÷òî
äèñïåðñèÿ íà
âõîäå â ïåðâóþ
íåëèíåéíîñòü
èìååò âèä:
Íà
âõîäå âòîðîé
íåëèíåéíîñòè:
Äàëåå
íåîáõîäèìî
ïîëó÷èòü
âûðàæåíèÿ
äëÿ ñàìèõ
êîýôôèöèåíòîâ
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè,
âîñïîëüçîâàâøèñü
âûâåäåííûìè
ðàíüøå ñîîòíîøåíèÿìè.
Ëåãêî âèäåòü,
÷òî âñå èíòåãðàëû
â ýòèõ
ôîðìóëàõ
áóäóò èìåòü
îäèí èç òðåõ,
íèæå
ïåðå÷èñëåííûõ,
âèäîâ:
Ýòè
èíòåãðàëû,
ñ÷èòàþòñÿ â
÷èñëåííîì
âèäå è ïîëó÷àþòñÿ
ñ ïîìîùüþ
ôóíêöèè
îøèáîê è
ãàóññîâñêîé
ïëîòíîñòè
âåðîÿòíîñòè.
Îíè ðåàëèçîâàíû
â ïðîãðàììå â
âèäå
ôóíêöèé,
òîãäà êîýôôèöèåíòû
ñòàòèñòè÷åñêîé
ëèíåàðèçàöèè
äëÿ ïåðâîãî
íåëèíåéíîãî
ýëåìåíòà áóäóò
èìåòü âèä:
K0=k2*J1(0,-s,m,D,-1)+l1*J0(0,-s,m,D,-1)+
k1*J1(-s,s,m,D,0)+k2*J1(s,0,m,D,1)+l2*J0(s,0,m,D,1);
K1=(k2*J2(0,-s,m,D,-1)+l1*J1(0,-s,m,D,-1)+k1*J1(-s,s,m,D,0)+
k2*J2(s,0,m,D,1)+l2*J1(s,0,m,D,1)-m*K0)/D;
Äëÿ
âòîðîãî:
K0=-s*J0(0,-s,m,D,-1)+J1(-s,s,m,D,0)+s*J0(s,0,m,D,1);
K1=(-s*J1(0,-s,m,D,-1)+J2(-s,s,m,D,0)+s*J1(s,0,m,D,1)+
s*m*J0(0,-s,m,D,-1)-m*J1(-s,s,m,D,0)-s*m*J0(s,0,m,D,1))/D;
Äëÿ
òðåòüåãî:
K0=-s*J0(0,-s,m,D,-1)+J1(-s,s,m,D,0)+s*J0(s,0,m,D,1);
K1=(-s*J1(0,-s,m,D,-1)+J2(-s,s,m,D,0)+s*J1(s,0,m,D,1)+
s*m*J0(0,-s,m,D,-1)-m*J1(-s,s,m,D,0)-s*m*J0(s,0,m,D,1))/D;
Ëèíåàðèçîâàííàÿ
ñèñòåìà
äîëæíà
èìåòü âèä:
Âîñïîëüçîâàâøèñü
óðàâíåíèÿìè
äëÿ íåëèíåéíîé
ñèñòåìû
ìàòðèöû
ìîæíî
çàïèñàòü:
Ëåãêî
âèäåòü, ÷òî
ìàòðèöà A íà
êàæäîì øàãå
èíòåãðèðîâàíèÿ
áóäåò èçìåíÿòüñÿ,
â
çàâèñèìîñòè
îò
êîýôôèöèåíòîâ
ëèíåàðèçàöèè,
êîòîðûå â
ñâîþ î÷åðåäü
çàâèñÿò îò
äèñïåðñèé
íà âõîäå
íåëèíåéíûõ
ýëåìåíòîâ,
êîòîðûå
çàâèñÿò îò
äèñïåðñèé
íà âõîäå.
Ñëåäîâàòåëüíî
íà êàæäîì
øàãå
èíòåãðèðîâàíèÿ
íàøåé ñèñòåìû
ìû äîëæíû
èíòåãðèðîâàòü
äèôôåðåíöèàëüíîå
óðàâíåíèå
äëÿ ìàòðèöû
êîâàðèàöèé K:
Ðåàëèçóÿ
âñå
âûøåñêàçàííîå
ïîëó÷èì ãðàôèê
èçìåíåíèÿ
âåêòîðà
ñîñòîÿíèé
ëèíåàðèçîâàííîé
ñèñòåìû:
Ðèñóíîê 2.
Ýâîëþöèÿ
ìàòðèöû
êîâàðèàöèè
ïî âõîäíûì âîçäåéñòâèÿì
èìååò âèä:
Ðèñóíîê 3
Äèñïåðñèÿ
âûõîäíîé
êîîðäèíàòû:
Ðèñóíîê 4
4.
Заключение
Êàê
âèäíî èç
ðèñóíêîâ 1 è 2
ñòàòèñòè÷åñêàÿ
ëèíåàðèçàöèÿ
ïðîâåäåíà
ïðàâèëüíî.
Ïîâåäåíèå
ñèñòåìû â
íåëèíåéíîé
è
ëèíåàðèçîâàííîé
ôîðìå
ïðèìåðíî îäèíàêîâî.
Èç ðèñóíêîâ 3
è 4 ñëåäóåò,
÷òî äèñïåðñèè
ñõîäÿòñÿ, ò. å.
÷åðåç
íåêîòîðûé
ïðîìåæóòîê
âðåìåíè
ïåðåõîäíûå
ïðîöåññû â
ñèñòåìå
çàêàí÷èâàþòñÿ
è íà âûõîäå
ñèñòåìû
èìååòñÿ ñòàöèîíàðíûé
ñëó÷àéíûé
ïðîöåññ.
5.
Список
литературы
1.
Êîíñïåêò
ëåêöèé
2.
Ëåáåäåâ
è äð.
Ñòàòèñòè÷åñêàÿ
äèíàìèêà è îïòèìèçàöèÿ
óïðàâëåíèÿ
ËÀ.
3.
Êîðí Ã. è
Ò.
Ñïðàâî÷íèê
ïî
ìàòåìàòèêå
4.
Ïîëÿêîâ Þ.
Â, Êðóãëîâ È. Þ.
Ïðîãðàììèðîâàíèå
â ñðåäå Òóðáî Ïàñêàëü
5.5
5.
Áàáàê Ñ.
Â., Âàñèëüåâ Â.
È. è äð. Îñíîâû
òåîðèè ìíîãîñâÿçíûõ
ñèñòåì
àâòîìàòè÷åñêîãî
óïðàâëåíèÿ
ËÀ.
//КуЁёютря
Ёрсютр яю ётртшётшчхёъющ фшэрьшъх ЛА
//Выяюыэшы
туфхэт уЁ. 07-403
//КюЁэшыют
Д. М.
//Ррёчхт
ышэхрЁшчютрээющ ёшётхьы
//ЛшэхрЁшчртрээря
шётхьр
#include
<math.h>
#include
<iostream.h>
#include
<fstream.h>
#include
<stdlib.h>
#include
<time.h>
#include
"matrix.h"
#include
"runge.h"
#include
"random.h"
//Зрфрэшх
const
double Da=0.1;
const
double Dn=0.03;
const
double K1=2.5;
const
double K2=1.9;
const
double l=1.2;
const
double Kp=13;
const
double S1=0.5;
const
double S2=27;
const
double S3=3;
const
double dh=0.01;
//Кырёё
тхътюЁр ёюётюяэшя схч ъютрЁшрцшш
template<class
T> class CMyRunge:public CRunge4<T>{
public:
CMyRunge(TMatrix<T> t,T h, T
x):CRunge4<T>(t,h,x){};
TMatrix<T> f(T,TMatrix<T>);
};
//Кырёё
ьртЁшцы ъютрЁшрцшш
template<class
T> class CKovariation:public CRunge4<T>{
public:
CKovariation(TMatrix<T> t,T h, T
x):CRunge4<T>(t,h,x){};
TMatrix<T> f(T,TMatrix<T>);
};
//МртЁшцы
уЁртэхэшя
TMatrix<double>
A(3,3);
TMatrix<double>
B(3,1);
//ВхътюЁ
юётюяэшя
TMatrix<double>
X(3,1);
//МртЁшцр
ъютрЁшрцшш
TMatrix<double>
K(3,3);
//Сыучрщэыщ
тхътюЁ
TMatrix<double>
s(3,1);
//Кюэффшцшхэты
трт.ышэхрЁшчрцшш
double
K11,K21,K31;
//ВючтЁрщрхт
чэръ x
double
Sign(double x)
{
if(x<0)
return -1; else return 1;
}
//ИэтхуЁры
юшшсъш
double
erf(double x)
{
double x1=fabs(x);
const
double p=0.3275911;
const
double a1=0.254829592;
const
double a2=-0.284496736;
const
double a3=1.421413741;
const
double a4=-1.453152027;
const
double a5=1.061405429;
double
t=1/(1+p*x1);
return
Sign(x)*(1-(a1*t+a2*pow(t,2)+a3*pow(t,3)+a4*pow(t,4)+a5*pow(t,5))*exp(-x*x));
}
//Пыютэюёть
тхЁюятэюётш
double
KGauss(double x,double m,double D)
{
return
1/sqrt(2*3.14159*D)*exp(-1/2*(x-m)*(x-m)/D);
}
//ИэтхуЁрыы
//ИэтхуЁры
ют p(x)
double
J0(double Low, double Up, double m, double D, int Inf)
{
switch( Inf) {
case -1:return
(1+erf((Up-m)/sqrt(2*D)))/2;
case 0: return
(erf((Up-m)/sqrt(2*D))-erf((Low-m)/sqrt(2*D)) )/2;
case +1:return
(1-erf((Low-m)/sqrt(2*D)))/2;
}
return 0;
}
//ИэтхуЁры
ют x*p(x)
double
J1(double Low, double Up, double m, double D, int Inf)
{
switch(Inf){
case
-1:return (-D*KGauss(Up,m,D)+m*J0(0,Up,m,D,Inf));
case 0:
return(-D*(KGauss(Up,m,D)-KGauss(Low,m,D))+m*J0(Low,Up,m,D,Inf) );
case 1:
return(D*KGauss(Low,m,D)+m*J0(Low,0,m,D,Inf));
}
return
0;
}
//ИэтхуЁры
ют x*x*p(x)
double
J2(double Low, double Up, double m, double D, int Inf)
{
switch(Inf){
case
-1:return (-D*(Up*KGauss(Up,m,D)-J0(0,Up,m,D,Inf))+m*J1(0,Up,m,D,Inf));
case 0:
{
double d;
d=-D*(Up*KGauss(Up,m,D)-Low*KGauss(Low,m,D)-
J0(Low,Up,m,D,Inf))+m*J1(Low,Up,m,D,Inf);
return d;}
case 1:
return(-D*(-Low*KGauss(Low,m,D)-J0(Low,0,m,D,Inf))+m*J1(Low,0,m,D,Inf));
}
return
0;
}
//Стрт.
ышэхрЁшчрцшя
//Кюэффшцшхэты
яхЁтющ эхышэхщэюётш
void
Linear1(double m, double D, double& K0,double& K1)
{
const
double s=0.5;
const
double k1=2.5;
const
double k2=1.9;
double
l1=s*(k2-k1);
double
l2=-l1;
K0=k2*J1(0,-s,m,D,-1)+l1*J0(0,-s,m,D,-1)+
k1*J1(-s,s,m,D,0)+k2*J1(s,0,m,D,1)+l2*J0(s,0,m,D,1);
K1=(k2*J2(0,-s,m,D,-1)+l1*J1(0,-s,m,D,-1)+k1*J1(-s,s,m,D,0)+
k2*J2(s,0,m,D,1)+l2*J1(s,0,m,D,1)-m*K0)/D;
}
//Кюэффшцшхэты
ттюЁющ эхышэхщэюётш
void
Linear2(double m, double D, double& K0,double& K1)
{
const
double s=27;
K0=-s*J0(0,-s,m,D,-1)+J1(-s,s,m,D,0)+s*J0(s,0,m,D,1);
K1=(-s*J1(0,-s,m,D,-1)+J2(-s,s,m,D,0)+s*J1(s,0,m,D,1)+
s*m*J0(0,-s,m,D,-1)-m*J1(-s,s,m,D,0)-s*m*J0(s,0,m,D,1))/D;
}
//Кюэффшцшхэты
тЁхтьхщ эхышэхщэюётш
void
Linear3(double m, double D, double& K0,double& K1)
{
const
double s=3;
K0=-s*J0(0,-s,m,D,-1)+J1(-s,s,m,D,0)+s*J0(s,0,m,D,1);
K1=(-s*J1(0,-s,m,D,-1)+J2(-s,s,m,D,0)+s*J1(s,0,m,D,1)+
s*m*J0(0,-s,m,D,-1)-m*J1(-s,s,m,D,0)-s*m*J0(s,0,m,D,1))/D;
}
//ДшёяхЁёшш
эр тхюфх эхышэхщэюётхщ
//ПхЁтря
//ПрЁрьхтЁы:
ъюэффшцшхэты ышэхрЁшчрцшш яхЁтющ ш ттюЁющ эхышэхщэюётхщ
double
Dx(double k1,double k2)
{
double
d1,d2;
d1=1/(2*3.14159)*Kp*Kp*K.Element(1,1)*l/(l+Kp*k1*k2);
d2=1/(2*3.14159)*Kp*Kp*K.Element(2,2)/(k1*k2-l)/k1/k2;
return
(d1+d2);
}
//ВтюЁря
double
Dy(double k1,double k2)
{
double
d1,d2;
d1=1/(2*3.14159)*k1*k1*Kp*Kp*K.Element(1,1)*l/(l+Kp*k1*k2);
d2=1/(3.14159)*K1*Kp*K.Element(2,2)/(l-k1*k2*Kp)/k1/k2;
return
(d1+d2);
}
//Рхрышчрцшя
ъырёёр тхътюЁр ёюётюяэшя
template
<class T> TMatrix<T> CMyRunge<T>::f(T h,TMatrix<T> m)
{
double d=1;
TMatrix<T>*
t=new TMatrix<T>;
s.Element(1,1,B.Element(1,1)*gauss(0,d));
s.Element(2,1,B.Element(2,1)*gauss(0,d));
s.Element(3,1,B.Element(3,1));
*t=A*m+s;
return
*t;
}
//Рхрышчрцшя
ъырёёр ьртЁшцы ъютрЁшрцшш
template
<class T> TMatrix<T> CKovariation<T>::f(T h,TMatrix<T>
m)
{
TMatrix<T>*
t=new TMatrix<T>(3,3);
*t=A*m+m*A.Transpose()+B*B.Transpose();
return
*t;
}
//Пюххрыш
main()
{
//Иэшцшрышчрцшя
ьртЁшцы A
for(int
i=1; i<=3;i++)
for(int
j=1; j<=3;j++)
A.Element(i,j,0);
A.Element(1,1,-l);
A.Element(2,2,-l);
A.Element(3,2,1);
//Иэшцшрышчрцшя
ьртЁшцы B
B.Element(1,1,sqrt(Da*l/3.14159));
B.Element(2,1,sqrt(Dn*l/3.14159));
B.Element(3,1,0);
//Иэшцшрышчрцшя
ьртЁшцы ъютрЁшрцшш
for(i=1;
i<=3;i++)
for(int
j=1; j<=3;j++)
if(i!=j)K.Element(i,j,0);
K.Element(1,1,Da);
K.Element(2,2,Dn);
K.Element(3,3,(Da*Kp*Kp+Dn)*10);
//Иэшцшрышчрцшя
тхътюЁр ёюётюяэшя X
for(i=1;
i<=3;i++)
X.Element(i,1,0);
//Кырёё
Руэух фыя тхътюЁр ёюётюяэшя
CMyRunge<double>
Vector(X,dh,0);
//Кырёё
Руэух фыя ъютрЁшрцшш
CKovariation<double>
Kov(K,dh,0);
//ЗртЁртър
double
k10,k11,k20,k21;
Linear1(0,Da*Kp*Kp,k10,k11);
Linear2(0,Da*Kp*Kp,k20,k21);
double
d1=Dx(k11,k21);
double
d2=Dy(k11,k21);
Linear1(0,d1,k10,k11);
Linear2(0,d2,k20,k21);
K11=k11;
K21=k21;
Linear3(0,d2+Dn,k10,k11);
K31=k11;
A.Element(3,1,K11*K21*Kp);
A.Element(3,3,-K11*K21*Kp);
int
first=1;
filebuf
f;
f.open("x.dat",ios::out
);
ostream
o(&f);
filebuf
f1;
f1.open("k.dat",ios::out
);
ostream
o1(&f1);
double
timec=0;
double
dt=2.0/100;
//Оёэютэющ
цшъы
for(i=1;i<=200;i++)
{
timec+=dt;
//Зряюыэяхь ьртЁшцу A
if
(!first){
d1=Dx(K11,K21);
d2=Dy(K11,K21);
Linear1(0,d1,k10,K11);
Linear2(0,d2,k20,K21);
Linear3(0,K.Element(3,3),k10,K31);
A.Element(3,1,K11*K21*Kp);
A.Element(3,3,-K11*K21*Kp);
first=0;
};
X=Vector.Next();
K=Kov.Next();
o<<X.Element(1,1)<<"\t"<<X.Element(2,1)<<"\t"<<
K31*X.Element(3,1)<<"\t"<<timec<<endl;
o1<<K.Element(1,1)<<"\t"<<K.Element(2,2)<<"\t"<<
K31*K31*K.Element(3,3)<<"\t"<<timec<<endl;
}
cout<<K;
cin>>i;
return
0;
}
//НхышэхрЁшчютрээря
шётхьр
#include
<iostream.h>
#include
<fstream.h>
#include
<stdlib.h>
#include
<time.h>
#include
"matrix.h"
#include
"runge.h"
#include
"random.h"
//Зрфрэшх
const
double Da=0.1;
const
double Dn=0.03;
const
double K1=2.5;
const
double K2=1.9;
const
double l=1.2;
const
double Kp=13;
const
double S1=0.5;
const
double S2=27;
const
double S3=3;
const
double dh=0.001;
template<class
T> class CMyRunge:public CRunge4<T>{
public:
CMyRunge(TMatrix<T> t,T h, T
x):CRunge4<T>(t,h,x){};
TMatrix<T> f(T,TMatrix<T>);
static T F1(T);
static T F2(T);
static T F3(T);
protected:
T k1,k2;
};
template
<class T> TMatrix<T> CMyRunge<T>::f(T h,TMatrix<T> m)
{
double d=10;
k1=sqrt(Da*l/3.14159);
k2=sqrt(Dn*l/3.14159);
TMatrix<T>*
t=new TMatrix<T>(3,1);
t->Element(1,1,-l*m.Element(1,1)+k1*gauss(0,d));
t->Element(2,1,-l*m.Element(2,1)+k2*gauss(0,d));
t->Element(3,1,m.Element(2,1)+F2(F1(Kp*(m.Element(1,1)-m.Element(3,1)))));
return
*t;
}
template
<class T> T CMyRunge<T>::F1(T x)
{
if(fabs(x)>=S1)
return
(K2*x+x/fabs(x)*S1*(K1-K2)); else
return(K1*x);
return
0;
}
template
<class T> T CMyRunge<T>::F2(T x)
{
if(x<=-S2)
return
(-S2);
if(fabs(x)<S2)
return(x);
if(x>=S2)
return
(S2);
return
0;
}
template
<class T> T CMyRunge<T>::F3(T x)
{
if(x<=-S3)
return
(-S3);
if(fabs(x)<S3)
return(x);
if(x>=S3)
return
(S3);
return
0;
}
main()
{
int i;
TMatrix<double>
m(3,1);
m.Element(1,1,0);
m.Element(2,1,0);
m.Element(3,1,0);
CMyRunge<double>
r(m,dh,0);
filebuf
f;
f.open("a.dat",ios::out
);
ostream
o(&f);
double
t=0;
double
dt=2.0/100;
for(i=1;i<=200;i++){
t+=dt;
m=r.Next();
o<<m.Element(1,1)<<"\t"<<m.Element(2,1)<<
"\t"<<CMyRunge<double>::F3(m.Element(3,1))<<"\t"<<t<<endl;
//cout<<m.Element(1,1)<<"\t"<<m.Element(2,1)<<
//"\t"<<CMyRunge<double>::F3(m.Element(3,1))<<endl;
}
cin>>i;
return
0;
}
//Matrix
Template Library, version 1.00
//(c)
by Dmitry Kornilov, 1997
//Designed
for Microsoft Visual C++ and Borland C++
#ifndef
_IOSTREAM_H
#include
<iostream.h>
#endif
#ifndef
_MATRIX_H
#define
_MATRIX_H
#endif
#define
INDEX(i,j,col) (i-1)*col+j-1
template
<class T> class TMatrix
{public:
TMatrix(int ,int );
TMatrix(){col=row=0;buffer=NULL;};
TMatrix(TMatrix<T>&);
~TMatrix();
void
Element(int,int,T);
T
Element(int,int);
TMatrix<T>
Minor(int, int);
TMatrix<T>
Invert();
TMatrix<T>
Transpose();
T
Determinant();
int
GetRows(){return row;}
int
GetCols(){return col;}
friend
TMatrix operator+(TMatrix<T>&,TMatrix<T>&);
friend
TMatrix operator+(TMatrix<T>&,T);
friend
TMatrix operator-(TMatrix<T>&,TMatrix<T>&);
friend
TMatrix operator*(T,TMatrix<T>&);
friend
TMatrix operator*(TMatrix<T>&,T);
friend
TMatrix operator*(const TMatrix<T>& ,const TMatrix<T>& );
friend
ostream& operator<< (ostream&, const TMatrix<T>&);
void
operator=(TMatrix&);
protected:
T
*buffer;
int
row,col;
};
class
TMatrixError
{public:
int
Error;
char*
Message;
TMatrixError(int
n) {Error=n; };
char*
GetString();
};
char*
TMatrixError::GetString()
{
if(Error==5)
Message="Index Error";
return
Message;
}
template
<class T> TMatrix<T>::TMatrix(int rows,int cols)
{
row=rows;
col=cols;
buffer=new
T[row*col];
}
template
<class T> TMatrix<T>::~TMatrix()
{
delete
[] buffer;
}
template
<class T> TMatrix<T>::TMatrix(TMatrix<T> &m)
{
row=m.row;
col=m.col;
buffer=new
T [row*col];
for(int
i=1;i<=row;i++)
for(int
j=1;j<=col;j++)
buffer[(i-1)*col+j-1]=m.buffer[(i-1)*m.col+j-1];
}
template
<class T>void TMatrix<T>::Element(int i,int j, T f)
{
buffer[(i-1)*col+j-1]=f;
}
template
<class T>T TMatrix<T>::Element(int i,int j)
buffer[(i-1)*col+j-1];
template
<class T> TMatrix<T>
operator+(TMatrix<T>& m1, TMatrix<T>& m2)
{
TMatrix<T>&
m= *new TMatrix<T>(m1.row,m1.col);
for(int
i=1;i<=m1.row;i++)
for(int
j=1;j<=m1.col;j++)
m.buffer[(i-1)*m1.col+j-1]=m1.buffer[(i-1)*m1.col+j-1]+
m2.buffer[(i-1)*m1.col+j-1];
//m.Element(i,j,m1.Element(i,j)+m2.Element(i,j));
return
m;
}
template
<class T> TMatrix<T>
operator+(TMatrix<T>& m1, T f)
{
TMatrix<T>&
m= *new TMatrix<T>(m1.row,m1.col);
for(int
i=1;i<=m1.row;i++)
for(int
j=1;j<=m1.col;j++)
m.buffer[(i-1)*m1.col+j-1]=m1.buffer[(i-1)*m1.col+j-1]+f;
//m.Element(i,j,m1.Element(i,j)+m2.Element(i,j));
return
m;
}
template
<class T> TMatrix<T>
operator-(TMatrix<T>& m1, TMatrix<T>& m2)
{
TMatrix<T>
*m;
m=new
TMatrix<T>(m1.row,m1.col);
for(int
i=1;i<=m1.row;i++)
for(int
j=1;j<=m1.col;j++)
m->buffer[(i-1)*m1.col+j-1]=m1.buffer[(i-1)*m1.col+j-1]-
m2.buffer[(i-1)*m1.col+j-1];
return
*m;
}
template
<class T> TMatrix<T>
operator*(T f, TMatrix<T>& m1)
{
TMatrix<T>
*m;
m=new
TMatrix<T>(m1.row,m1.col);
for(int
i=1;i<=m1.row;i++)
for(int
j=1;j<=m1.col;j++)
m->buffer[(i-1)*m1.col+j-1]=m1.buffer[(i-1)*m1.col+j-1]*f;
return
*m;
}
template
<class T> TMatrix<T>
operator*(TMatrix<T>& m1, T f)
{
TMatrix<T>&
m =*new TMatrix<T>(m1.row,m1.col);
for(int
i=1;i<=m1.row;i++)
for(int
j=1;j<=m1.col;j++)
m.buffer[(i-1)*m1.col+j-1]=m1.buffer[(i-1)*m1.col+j-1]*f;
return
m;
}
template
<class T> void TMatrix<T>::operator=(TMatrix<T>& m)
{
row=m.row;
col=m.col;
if(!buffer)
buffer=new T [row*col];
for(int
i=1;i<=row;i++)
for(int
j=1;j<=col;j++)
buffer[(i-1)*col+j-1]=m.Element(i,j);
}
template<class
T> TMatrix<T> operator*(const TMatrix<T>&
m1,const TMatrix<T>& m2)
{
int
i,j,k;
T
f3,f2,f1;
TMatrix<T>*
m =new TMatrix<T>(m1.row,m2.col);
for(i=1;i<=m1.row;i++)
for(j=1;j<=m2.col;j++)
{
f3=0;
for(k=1;k<=m1.col;k++){
f1=m1.buffer[INDEX(i,k,m1.col)];
f2=m2.buffer[INDEX(k,j,m1.col)];
f3=f3+f1*f2;}
m->Element(i,j,f3);
}
return *m;
}
template
<class T> ostream& operator<< (ostream& o, const
TMatrix<T>& m)
{
for(int
i=1;i<=m.row;i++){
for(int
j=1;j<=m.col;j++)
o<<m.buffer[INDEX(i,j,m.col)]<<"\t";
o<<endl;
}
o<<endl;
return
o;
}
template
<class T> TMatrix<T> TMatrix<T>::Minor(int i, int j)
{
TMatrix<T>
*m=new TMatrix<T>(row-1,col-1);
int
count=0;
for
(int i1=1;i1<=row;i1++)
for
(int j1=1;j1<=col;j1++)
if
((i1!=i) && (j1!=j))
{
m->buffer[count]=buffer[INDEX(i1,j1,col)];
count++;
}
return
*m;
}
template
<class T> T TMatrix<T>::Determinant()
{
if(
col==2) {
return
buffer[INDEX(1,1,2)]*buffer[INDEX(2,2,2)]-
buffer[INDEX(1,2,2)]*buffer[INDEX(2,1,2)];}
T f=0;
TMatrix<T>
p;
for(int
i=1;i<=col;i++){
p=Minor(i,1);
f=f+buffer[INDEX(i,1,col)]*p.Determinant()*sign(i+1);
}
return
f;
}
int
sign(int i)
{
if( (i
% 2)==0) return 1; else
return
-1;
}
template
<class T> TMatrix<T> TMatrix<T>::Invert()
{
TMatrix<T>
*m=new TMatrix<T>(row,col);
T
f=Determinant();
for
(int i=1;i<=row;i++)
for
(int j=1;j<=col;j++)
m->buffer[INDEX(i,j,col)]=sign(i+j)*(Minor(j,i).Determinant())/f;
return
*m;
}
template
<class T> TMatrix<T> TMatrix<T>::Transpose()
{
TMatrix<T>
*m=new TMatrix<T>(col,row);
for
(int i=1;i<=row;i++)
for
(int j=1;j<=col;j++)
m->buffer[INDEX(j,i,row)]=buffer[INDEX(i,j,col)];
return
*m;
}
//Random
functions
//randomize-чрфрэшх
эютющ чртЁртъш
//random-ёыучрщэюх
чшёыю т фшрярчюэх 0-1
//gauss-эюЁьрыьэюх
ЁрёяЁхфхыхэшх ё ьрт юцшфрэшхь
// m ш
Ё ътрфЁртшчхёъшь ютъыюэхэшхь sigma
#ifndef
_INC_MATH
#include
<math.h>
#endif
#include
<stdlib.h>
#include
<TIME.h>
void
randomize()
{
srand(
(unsigned)time( NULL ) );
}
double
random()
{
double*
d=new double;
*d=(double)rand()/RAND_MAX;
return
*d;
}
double
gauss(double m,double sigma)
{
double
a,b,r,sq;
do{
a=2*random()-1;
b=2*random()-1;
r=pow(a,2)+pow(b,2);
}while
(r>1);
sq=sqrt(-2*log(r)/r);
return(m+sigma*a*sq);
}
//Runge-Kutt
template
//version
1.00
//By
Dmitry Kornilov
#ifndef
_MATRIX_H
#include
"matrix.h"
#endif
template
<class T> class CRunge4 {
public:
CRunge4(){};
CRunge4(TMatrix<T>,T,T);
virtual TMatrix<T>
f(T,TMatrix<T>)=0;
TMatrix<T> Next();
private:
TMatrix<T>
prev;
T h;
T
xprev;
};
template
<class T> CRunge4<T>::CRunge4(TMatrix<T> m,T f,T x)
{
h=f;
prev=m;
xprev=x;
}
template
<class T> TMatrix<T> CRunge4<T>::Next()
{
TMatrix<T>&
m= *new TMatrix<T>(prev.GetRows(),prev.GetCols());
TMatrix<T>
k1,k2,k3,k4;
k1=f(xprev,prev)*h;
k2=f(xprev+h/2,prev+k1*0.5)*h;
k3=f(xprev+h/2,prev+k2*0.5)*h;
k4=f(xprev+h,prev+k3)*h;
m=prev+((T)1/6)*(k1+2.0*k2+2.0*k3+k4);
prev=m;
xprev+=h;
return
m;
}
|
