M.Luft, Z.A
ukasik
340
10. Symulacja układów dynamicznych
10.1. Zasady symulacji analogowej
Badanie dynamiki układów z zakresu automatyki, budowy maszyn,
techniki lotniczej i rakietowej, procesów chemicznych i innych jest
zagadnieniem bardzo złożonym. Nawet badanie stabilności układów
dynamicznych opisanych liniowymi równaniami różniczkowymi niż-
szych rzędów jest stosunkowo uciążliwe. Badanie jakości działania
układów opisywanych równaniami różniczkowymi wyższych rzędów
jest jeszcze bardziej złożone. Trudności te wzrastają wielokrotnie przy
badaniu układów nieliniowych, gdzie przeprowadzenie teoretycznych
obliczeń jest często niemożliwe. Nowoczesne badania powyższych
problemów są prowadzone innymi metodami, które skojarzone z teo-
rią zapewniają rozwiązywanie tych zagadnień w sposób szybki i bar-
dzo wygodny w praktycznej realizacji. SÄ… to metody modelowania.
Zasada modelowania polega na wykrywaniu zależności charakteryzu-
jących badany układ za pomocą analizy przeprowadzonej na innym
układzie, zwanym modelem.
Procesy zachodzące w modelu można obserwować, rejestrować
oraz badać z różnych interesujących punktów widzenia. Tak więc mo-
delowanie pozwala rozwiązywać eksperymentalnie zadania dotyczące
przebiegów różnych procesów fizycznych, bez użycia obiektów rze-
czywistych.
IstniejÄ… dwie zasadnicze metody modelowania, a mianowicie:
- modelowanie fizyczne,
- modelowanie matematyczne.
W modelach fizycznych istnieje zgodność równań matematycznych
opisujących układ rzeczywisty i model. Odpowiednie części modelu
oraz układu rzeczywistego zachowują ren sam charakter fizyczny ba-
danego eksperymentu.
W modelu matematycznym, podobnie jak i fizycznym, istnieje
zgodność równań matematycznych opisujących układ rzeczywisty
i model. Różnica polega natomiast na możliwości zmiany charakteru
M.Luft, Z.A
ukasik
341
fizycznego modelu w stosunku do modelu rzeczywistego. Korzysta
się z faktu, że różne pod względem fizycznym zjawiska, są opisywane
analogicznymi równaniami matematycznymi. Jest to tzw. analogia
matematyczna.
Jeżeli układ rzeczywisty jest opisany np. równaniem różniczko-
wym:
n
d x ( n )
&ð x
f (t, x, x, &ð&ð,..., x ) (10.1)
n
dt
przy czym:
k
( k ) d x
x
k
dt
lub też układem n równań pierwszego rzędu o n niewiadomych typu:
dx
i
f (t, x , x ,..., xn ) (i , , ,..., n) (10.2)
i
dt
wystarcza również, aby zmienne Xi modelu były odwzorowaniem
zmiennych rzeczywistych x1 w odpowiedniej skali, tzn. model byłby
opisany równaniem lub układem równań różniczkowych:
dX
i
f ( , X , X , ... , X ) (i 1,2 ,3, ... , n) (10.3)
i 1 2 n
d
gdzie:
X xi ,
i x
i
= t ,
- współczynnik skali amplitudy,
- współczynnik skali czasu.
W modelach elektronicznych wielkości wyrażenia (10.3) X1, X2, . .
.,Xn są przedstawione w postaci odpowiednich napięć. W modelach
M.Luft, Z.A
ukasik
342
może być również zastosowana różna skala czasu. W ten sposób po-
wolne procesy zachodzące w układach rzeczywistych, trwające kilka
czy kilkanaście godzin, mogą w modelach matematycznych trwać
ułamki sekund i odwrotnie. Przy modelowaniu całych układów stosuje
się takie współczynniki skali czasu, przy których można uzyskać mi-
nimalne błędy, w przypadku problemów ilościowych. Przy badaniach
jakościowych natomiast przyjmuje się możliwie najmniejsze współ-
czynniki skali czasu, zapewniając krótki czas oczekiwania na rozwią-
zanie oraz zmniejszenie czasu zużycia maszyny.
Modelowanie matematyczne za pomocÄ… maszyn analogowych mo-
że być rozwiązane dwiema metodami:
- metodÄ… modelowania strukturalnego,
- metodą modelowania równań.
Pierwsza z nich jest znacznie prostsza i bardziej przejrzysta. Polega
ona na modelowaniu problemów, które są podane w postaci schematu
strukturalnego i równań opisujących jego poszczególne bloki. Metoda
modelowania równań jest mniej przejrzysta i stosuje się ją wówczas,
gdy dynamika badanego układu jest podana tylko w postaci równania
różniczkowego lub układu równań.
10.1.1. Wzmacniacz operacyjny
Podstawowym elementem maszyn analogowych jest wzmacniacz
operacyjny. Służy on do wykonywania większości matematycznych
operacji liniowych i stanowi element pomocniczy urządzeń i układów
nieliniowych. Wzmacniacz operacyjny jest to wzmacniacz napięcia
stałego, charakteryzujący się następującymi cechami:
- bardzo dużym współczynnikiem wzmocnienia napięciowego
(rzędu 104 - 106 V/V),
- bardzo dużą impedancją wejściową i małą impedancją wyj-
ściową,
- dobrÄ… stabilizacjÄ… poziomu zerowego,
- dostatecznie szerokim pasmem przenoszenia częstotliwości,
- dwupolarnÄ… charakterystykÄ… statycznÄ…,
- odwracalną fazą sygnału wyjściowego o 180o (zmiana znaku na
przeciwny).
M.Luft, Z.A
ukasik
343
Duże zastosowanie znajdują ostatnio elektroniczne wzmacniacze ope-
racyjne zbudowane na elementach scalonych, charakteryzujÄ…ce siÄ™
lepszymi parametrami, niż produkowane wcześniej wzmacniacze
tranzystorowe lub lampowe.
Właściwości wzmacniacza operacyjnego wynikają z zastosowania
silnego ujemnego sprzężenia zwrotnego. Rozpatrzmy układ przedsta-
wiony na rysunku 10.1 przyjmując, że wzmacniacz ma nieskończenie
duże wzmocnienie, nieskończenie dużą impedancję wejściową i zero-
wą impedancję wyjściową.
Ze względu na dużą wartość impedancji wejściowej wzmacniacza
i0 = i1 - i2 0 (10.4)
Z2(s)
i2
Z1(s)
U (s) U (s) U (s)
x 0 y
-Ku
i1 i0
Rys. 10.1. Uproszczony schemat ideowy jednowejściowego wzmac-
niacza operacyjnego: Ux - napięcie wejściowe, Uy - napięcie
wyjściowe, Z1 - impedancja obwodu wejściowego, Z2 - im-
pedancja sprzężenia zwrotnego, Ku - współczynnik wzmoc-
nienia
Ponieważ
U ( s)
y
U ( s) (10.5)
0
Ku
zatem na podstawie wyrażenia (10.4) można napisać:
M.Luft, Z.A
ukasik
344
U ( s) U ( s)
y y
U ( s) U ( s)
x y
K K
u u
(10.6)
Z1 ( s) Z2 ( s)
skąd po przekształceniu otrzymujemy wyrażenie transmitancji opera-
torowej układu z rysunku 10.1:
U (s)
Z (s) 1
y
2
(10.7)
U (s) Z1 (s)
1 Z (s)
x
2
1 1
K Z1(s)
u
Wyrażenie (10.7) przy założeniu Ku , upraszcza się do postaci:
U ( s)
Z2 ( s)
y
(10.8)
U ( s) Z1 ( s)
x
Z wyrażenia (10.8) wynika, że w zależności od wyboru impedancji Z1
i Z2 wzmacniacz operacyjny może realizować różne funkcje.
Z1n(s)
U (s)
1n
:
:
:
:
Z12(s) Z2(s)
U (s)
12
Z11(s)
U (s) U (s)
11 y
-Ku
U (s)
0
Rys. 10.2. Schemat ideowy sumatora
Wyrażenie opisujące wzmacniacz z wieloma wejściami (rys. 10.2)
można wyprowadzić w taki sam sposób jak wyrażenie (10.6). W tym
przypadku
M.Luft, Z.A
ukasik
345
n
i0 iw i2 iij i2 0 (10.9)
j 1
oraz
n
U ( s)
x 1
U ( s) Z2 ( s) (10.10)
y
Z1 j ( s)
j 1
Ze wzoru (10.10) wynika, że wzmacniacz z wieloma wejściami speł-
nia rolÄ™ sumatora.
Szeregowe połączenie układów podstawowych z rysunków 10.1 i
10.2 dzięki dużej impedancji wejściowej i małej impedancji wyjścio-
wej wzmacniaczy operacyjnych, pozwala zrealizować układ, którego
transmitancja jest równa iloczynowi transmitancji połączonych w ten
sposób wzmacniaczy (rys. 10.3);
n
U ( s) Z2 j ( s)
y
( 1)n (10.11)
U ( s) Z1 j ( s)
j 1
x
Z21(s) Z22(s) Z2n(s)
Z (s)
Z (s) Z (s) 1n
11 12 U (s)
U (s) y
x
Rys. 10.3. Szeregowe połączenie układów zrealizowanych na wzmac-
niaczach operacyjnych
10.1.2. Modelowanie strukturalne układów dynamicznych
na wzmacniaczach operacyjnych
Odpowiedni dobór impedancji sprzęgających we wzmacniaczu
operacyjnym umożliwia realizację różnych liniowych operacji mate-
matycznych, takich jak: sumowanie, odejmowanie, mnożenie przez
stały współczynnik, całkowanie, różniczkowanie, itp. Włączając do
M.Luft, Z.A
ukasik
346
pętli sprzężenia zwrotnego impedancję nieliniowe, można realizować
różne operacje nieliniowe. Dla podstawowego wzmacniacza odwraca-
jącego fazę (rys 10.1) zestawiono w tablicy 10.1 najczęściej spotykane
impedancje liniowe Z1 i Z2. Zgodnie ze wzorem (10.8), stosunek
dwóch dowolnie wybranych impedancji określa transmitancję
wzmacniacza. Rozpatrzmy kilka najczęściej spotykanych przypad-
ków.
Człon proporcjonalny
PrzyjmujÄ…c impedancje rzeczywiste Z1(s) = R1 i Z2(s) = R2 uzysku-
jemy transmitancjÄ™ wzmacniacza
G ( s) k (10.12)
przy czym:
R2
k .
R1
Z transmitancji (10.12) wynika, że jest to człon proporcjonalny
o współczynniku wzmocnienia k (rys. 10.4a). Analogicznie przyjmu-
jąc impedancje rzeczywiste w układzie z rysunku 10.2 otrzymamy
sumator, stanowiący niezbędny element wszelkich układów regulacji
(rys. 10.4b). PodstawiajÄ…c do wzoru (10.10) Z2(s) = Rij gdzie (i, j = 1,
2, . . . ,n) otrzymamy:
n
U k U (10.13)
y j xj
j 1
R R
2 2
a) b)
R R
1 11
U (s) U (s) U (s) U (s)
x y x1 y
:
:
:
R
1n
U (s
xn
)
Rys. 10.4. Człon proporcjonalny: a) schemat blokowy, b) sumator
M.Luft, Z.A
ukasik
347
Człon inercyjny
C
R
2
R
1
U (s) U (s)
x y
Rys. 10.5. Człon inercyjny pierwszego rzędu
Przyjmując impedancję Z2 jako równoległe połączenie pojemności C i
rezystancji R2 oraz jako impedancjÄ™ Z1 rezystancjÄ™ R1, otrzymamy
człon inercyjny pierwszego rzędu (rys. 10.5) o transmitancji:
k
G ( s) (10.14)
1 sT
przy czym:
R2
k , T R2C .
R1
A
ącząc szeregowo n układów z rysunku 10.5, uzyskuje się człon iner-
cyjny n-tego rzędu o transmitancji:
n
k
n j
G ( s) ( 1) (10.15)
1 sTj
j 1
M.Luft, Z.A
ukasik
348
Tablica 10.1
Impedancje liniowe Z(s) stosowane na wejściach i w pętlach sprzężeń zwrot-
nych wzmacniaczy operacyjnych
Impedancje
Z (s)
Lp. Obwód Zależności
Z (s)
Z (s)
1
R
k k = R
2
C
Ä… C
Ä… s
R
3
k R
C k
T s
T RC
4
C R
Ä… C
( sT )
Ä… s
T RC
5
k R
R
1
Ä… Ts
k R
C R
2
Ä…
T s
R R
T ( R R )C
6
Ä… C C
C R
1 1
T R C
1 (1 sT1 )(1 sT2 )
C R
2 2
C C
s
s 1 sT3
T ( R R )
C C
T R C
(T T T )
M.Luft, Z.A
ukasik
349
Integrator
Całkowanie jest ważną funkcją wzmacniacza operacyjnego. Układ
całkujący (integrator) otrzymuje się poprzez włączenie do pętli
wzmacniacza kondensatora, a na wejściu rezystora (rys. 10.6a).
Transmitancja takiego członu ma postać:
k
G ( s) (10.16)
s
przy czym:
1
k .
RC
W praktyce w integratorze należy zapewnić możliwość rozpoczęcia
całkowania z określonym warunkiem początkowym, czyli działanie
zgodnie z zależnością:
t
U (t ) k U ( )dx u0 (10.17)
y x
0
gdzie: u0 - warunek poczÄ…tkowy.
a)
C
R
U (s) U (s)
x y
r1 r2
b) u0
1
2
C
K
1
R
1
U (s) 1 U (s)
x y
K
2
2
Rys.10.6. Integrator: a) wzmacniacz w układzie całkującym,
b) wprowadzenie warunku poczÄ…tkowego
Istnieje wiele sposobów wprowadzenia niezerowego warunku począt-
kowego do integratora. Jeden z nich, stosowany w maszynach analo-
gowych, pokazano na rysunku 10.6b. Polega on na Å‚adowaniu konden-
M.Luft, Z.A
ukasik
350
satora w układzie inercyjnym, bezpośrednio do wartości napięcia wa-
runku początkowego. Jeżeli klucze K1, K2 znajdą się w położeniu (1),
to układ całkuje. Zmiana położenia kluczy do pozycji (2) powoduje
odłączenie napięcia wejściowego Ux i dołączenie napięcia warunków
początkowych u0. Kondensator jest ładowany ze stałą czasową r2C. Po
przełączeniu kluczy K1, K2 ponownie do pozycji (1), układ całkuje
napięcie poczynając od napięcia warunku początkowego.
Schematy i transmitancje podstawowych układów liniowych zrealizo-
wanych na wzmacniaczu operacyjnym zestawiono w tablicy 10.2.
Przedstawione proste układy dynamiczne, realizowane za pomocą
wzmacniaczy operacyjnych, są wykorzystywane najczęściej do mode-
lowania lub budowy układów dynamicznych danych w postaci schematu
blokowego. Jest to metoda modelowania strukturalnego, w którym każ-
demu elementowi schematu blokowego przyporzÄ…dkowuje siÄ™ jego model
analogowy zrealizowany na wzmacniaczu operacyjnym poprzez odpo-
wiedni dobór impedancji Z1 i Z2. Prosty przykład modelowania struktural-
nego pokazano na rysunku 10.7. W układzie tym występują elementy iner-
cyjne, całkujące, regulator PI oraz dwa węzły sumacyjne. Współczynniki
proporcjonalności oraz stałe czasowe poszczególnych elementów dobiera
się zgodnie z następującymi zależnościami:
Ri R3
k , k , k2 ,
p
R1 R4
R C
Ti RiCi , T R3C2 .
Jeżeli w omawianym modelu użyta zostanie naturalna skala czasu, to
przykładową wartość Ti = 100s można ustalić, dobierając Ri = 100M i
Ci = 1 F. Metoda modelowania strukturalnego jest powszechnie stosowa-
ną do analizy właściwości układów regulacji automatycznej, gdyż mode-
lowany układ nie traci swojej struktury blokowej. Gałąz
sprzężenia zwrot-
nego w układzie modelującym jest bezpośrednio widoczna, istnieje więc
możliwość rozwierania pętli sprzężenia zwrotnego i badania modelu ukła-
du otwartego. Bezpośrednia analogia między schematem blokowym a
analogowym, pozwalająca na oddzielne badania części składowych bez
zmiany jego struktury blokowej, uzasadnia ewentualne użycie więk-
szej ilości członów operacyjnych.
M.Luft, Z.A
ukasik
351
Tablica 10.2.
Schematy i transmitancje podstawowych układów liniowych realizowanych na
wzmacniaczach operacyjnych
Transmitancja
Lp. Schemat połączeń Rodzaj
U (s)
elementów
G (s)
U (s)
R
2
1 P
R
R
1
U U
1 2
R
C
2 I
R
U U
1 2
RCs
R
2
3 D
R
1
R R Cs
U C U
1 2
rzeczywisty
R R Cs
C
4 inercyjny
R
2
R
I rzędu
R
1
U U
1 2
R R Cs
R C
2
5 PI
R R
1
U U
1 2
R R Cs
C R
2
6 PD
R
R
1
U U
1 2
R Cs
R
C
2
C R
1 2
7
R
1
U U
1 2
R R C s
korekcyjny
R R C s
M.Luft, Z.A
ukasik
352
z
1 2 3
a)
-
1 x
k1 + k
2
+
k (t )
p
2
Ti s
s (Ts 1)
x0
-
C C
2 2
b)
R R C C R R
0 i i 1 3 3
R R R R R
+x0
0 1 2 4 3
-x
R R
0 4
+z
Rys.10.7. Modelowanie strukturalne układu regulacji: a) schemat blo-
kowy, b) schemat analogowy; 1 - regulator, 2 - człon wyko-
nawczy, 3 - obiekt regulacji
10.2. Zasady modelowania układów dynamicznych na
maszynie analogowej
Modelowanie symulowanego urzÄ…dzenia na podstawie jego opisu
matematycznego, a więc rozwiązywanie równań różniczkowych na
maszynie analogowej, jest zwykle trudniejsze od modelowania struk-
turalnego, zwłaszcza w przypadku równań nieliniowych.
Nowoczesna maszyna analogowa ma trzy podstawowe rodzaje
elementów, wykorzystywanych przede wszystkim przy modelowaniu
problemów liniowych (rys.10.8) :
- integrator, realizujÄ…cy funkcjÄ™ (rys.10.8a),
- sumator, realizujÄ…cy operacjÄ™ (rys.10.8b),
- potencjometr (rys.10.8c), realizujący operację mnożenia przez
dodatni współczynnik mniejszy od jedności.
M.Luft, Z.A
ukasik
353
a)
k1
x1
b)
k2
y(0) x2
: y
:
x(t) y(t)
:
k
kn
xn
c)
x
y = x
Rys. 10.8. Podstawowe liniowe elementy operacyjne maszyny analo-
gowej: a) integrator, b) sumator, c) potencjometr
Oprócz tych elementów w maszynie analogowej występują jeszcze
inne przyrządy przeliczające (bloki mnożące, dzielące, pierwiastkują-
ce) oraz bloki pozwalające na kształtowanie podstawowych nielinio-
wości.
Liniowe układy dynamiczne są opisane liniowymi równaniami
różniczkowymi o stałych lub zmiennych współczynnikach. Liniowe
równanie różniczkowe ma postać:
n n 1
d x d x dx
an n an 1 n 1 ... a1 a0 x f (t ) 0 (10.18)
dt dt dt
z warunkami poczÄ…tkowymi:
i
d x
xi (i 0,1,2,..., n 1) (10.19)
dti
t 0
gdzie:
f(t) - wymuszenie działające na układ (funkcja wymuszająca).
M.Luft, Z.A
ukasik
354
W równaniu o stałych współczynnikach ai = const (i = 0, 1, 2, ..., n).
W równaniu o zmiennych współczynnikach wszystkie lub część
współczynników ai zależą od czasu (tj. od zmiennej niezależnej).
W celu zbudowania modelu równania różniczkowego (10.18) na ma-
szynie analogowej, wygodnie jest wyznaczyć z tego równania naj-
wyższą pochodną i przedstawić go w postaci:
n n 1
d x dx d x
A0 x A1 ... An 1 n 1 An f (t ) (10.19)
dtn dt dt
Realizację modelu powyższego równania na maszynie analogowej
(przy użyciu tylko elementów z rysunku 10.8) przedstawiono na ry-
sunku 10.9. Równanie (10.19) można również przedstawić w postaci
układu równań pierwszego rzędu, uzyskując w ten sposób równania
integratorów. Podstawiając nowe zmienne.1
( n )
&ð &ð&ð
x1 x, x x, x3 x, x x
2 n 1
otrzymamy na podstawie (10.19) układ równań:
&ð
x1 x
2
&ð
x x3
2
Mð (10.20)
&ð
x x
n n 1
x A0 x1 A1 x ... An 1 x An f (t )
n 1 2 n
z których ostatnie jest tzw. równaniem sumatora a pozostałe równa-
niami integratorów. Układ równań (10.20) prowadzi również do mo-
delu analogowego z rysunku 10.9.
1
Należy zauważyć, ze jest to po prostu wybór zmiennych stanu (patrz rozdział 9).
M.Luft, Z.A
ukasik
355
Rys. 10.9. Ogólny model analogowy równania (10.19)
Modelowanie na maszynie analogowej układów dynamicznych
opisanych transmitancjÄ… operatorowÄ…:
Y ( s) bm sm bm 1sm 1 ... b1s b0
G ( s) (10.21)
X ( s) an sn an 1sn 1 ... a1s a0
sprowadza się do przekształcenia danej transmitancji w równanie róż-
niczkowe jednej zmiennej lub układ równań różniczkowych pierw-
szego rzędu. Mając daną transmitancję (10.21), można napisać odpo-
wiadające jej równanie różniczkowe:
n n
d y d y dy
a a ... a a y (t )
n n
n n
dt
dt dt
(10.22)
m m
d y d y dy
bm m bm m ... b b x(t )
dt
dt dt
gdzie:
x(t) - wielkość wejściowa (wymuszająca),
y(t) - wielkość wyjściowa danego układu dynamicznego.
M.Luft, Z.A
ukasik
356
Równanie (10.22) jest niewygodne do programowania na maszynie
analogowej, dlatego też należy je odpowiednio przekształcić. Jedną ze
stosowanych w tym przypadku metod jest sprowadzenie tego równa-
nia do dwóch równań równoważnych za pomocą zmiennej pomocni-
czej. Mnożąc licznik i mianownik transmitancji (10.21) przez nową
zmienną Z(s), możemy ją przedstawić w postaci dwóch równoważ-
nych równań operatorowych (dla an = 1):
n n
X (a ) (a s a s ... a s a ) Z (s) (10.23)
n n
oraz
m m
Y (a ) (bm s bm s ... b s b ) Z (s) (10.24)
Równania powyższe przedstawione w postaci różniczkowej mają po-
stać:
n n
d z d z dz
a ... a a z x (10.25)
n
n n
dt
dt dt
m m
d z d z dz
y bm bm ... b b z (10.26)
m m
dt
dt dt
Ogólną budowę modeli równań (10.25) i (10.26), czyli model analo-
gowy układu o transmitancji (10.21) pokazano na rysunku 10.10.
M.Luft, Z.A
ukasik
357
Rys. 10.10. Ogólny model analogowy układu o transmitancji (10.21)
10.3 Zasady symulacji cyfrowej układów dynamicznych
10.3.1. Numeryczne całkowanie równań różniczkowych
Symulacja cyfrowa układu dynamicznego opisanego równaniami
różniczkowymi jest równoznaczna z ich numerycznym całkowaniem
przy użyciu maszyny cyfrowej. Symulacja cyfrowa w porównaniu
z analogową ma wiele zalet. Do najważniejszych należy zaliczyć: du-
żą dokładność obliczeń, brak konieczności skalowania równań, mniej-
szą pracochłonność. Ma to szczególne znaczenie przy symulacji nieli-
niowych równań różniczkowych.
Omawiane poniżej proste metody numeryczne polegają na ogół na
przybliżeniu krzywej rozwiązania ścisłego za pomocą odcinka linii
prostej, jakkolwiek np. metoda Adamsa zastępuje krzywą odcinkami
paraboli. Sposób obliczania odcinków jest zmienny, zależnie od przy-
jętej metody. Dokładność obliczeń jest zwykle lepsza, jeżeli korzysta
M.Luft, Z.A
ukasik
358
się z wartości obliczonych w kilku poprzednich krokach, a nie jedynie
z nachylenia w ostatnim punkcie. Niektóre z tych ogólnych metod
umożliwiają ponowne przeliczenie każdego kroku w celu zmniejsze-
nia błędu początkowego.
Naszym zadaniem jest rozwiązanie równania różniczkowego
dx
f (t , x ) (10.27)
dt
z warunkiem początkowym x(t0) = x0. Zakładamy, przy tym, że funk-
cja f(t,x) spełnia warunek Lipschitza i istnieje jednoznaczne rozwiąza-
nie zagadnienia Cauchy ego. OznaczajÄ…c przyrost (krok) zmiennej
niezależnej t = h, mamy
tk t0 k h ( k 0,1,2 ,...) (10.28)
Obliczanie (ściślej mówiąc, przybliżone) wartości rozwiązania x(t) i
jego pochodnej x (t) wyznacza się ze wzorów:
x x(t ) x(t k h)
k k 0
(k 0,1,2,...) (10.29)
x f f (t , x ) x (t )
k k k k k
Różnicę xk 1 x(t ) między wartością obliczoną a dokładną, nazy-
k 1
wamy błędem aproksymacji. Większość metod numerycznych dla
obliczenia (k+1) kroku wymaga co najmniej informacji o krokach (k-
1) oraz k. Oznacza to, że na ogół pierwszy krok obliczany jest inaczej
niż pozostałe. W bardziej złożonych metodach rzadko wymaga się
znajomości więcej niż sześciu poprzedzających kroków. Z reguły nie
używa się więcej niż trzech.
M.Luft, Z.A
ukasik
359
10.3.2. Metody Eulera i Adamsa
Rozwinięcie w szereg Taylora w otoczeniu wybranego punktu wią-
że funkcję w tym punkcie z kolejnymi pochodnymi tej funkcji
w punkcie poprzednim:
dx d xk h d xk
h
k
xk xk h Kð (10.30)
dt ! !
dt dt
Zauważmy, że dla obliczenia (k+1) wartości, wymagana jest informa-
cja tylko o k-tym punkcie.
W obliczeniach praktycznych korzysta siÄ™ tylko z kilku pierwszych
wyrazów. Metoda Eulera uwzględnia jedynie dwa pierwsze wyrazy
szeregu Taylora.
dx
k
xk xk h (10.31)
dt
Biorąc pod uwagę, że
dx
k
xk f
k
dt
mamy
x x h f (10.32)
k k k
Metodę Eulera można stosować dla dostatecznie małych wartości h
i tylko dla kilku pierwszych kroków. Dla następnych punktów zaleca
się stosowanie metody dokładniejszej - np. Adamsa.
Metoda Adamsa opiera siÄ™ na ekstrapolacji parabolicznej krzywej,
a nie na ekstrapolacji liniowej. Sposób postępowania zapewnia, że
błąd pozostaje mały.
Równanie paraboli zastosowanej do ekstrapolacji ma postać:
M.Luft, Z.A
ukasik
360
x x a (t t ) b(t t ) (10.33)
k k k
przy czym współczynniki a, b należy wyznaczyć.
Różniczkując wyrażenie (10.33) mamy
dx
a b(t tk ) (10.34)
dt
Jeżeli t = tk , to fk = a, a jeżeli t = tk-1 , to
f a b (t t ) (10.35)
k k k
Biorąc powyższe pod uwagę i podstawiając do (10.35) t = tk+1 otrzy-
mujemy podstawowy wzór metody Adamsa:
f f
k k
xk xk h( f ) (10.36)
k
Można wykazać, że aby błąd był mały, drugie wyrażenie w nawiasie
równania (10.36) należy utrzymać małe w stosunku do fk. Ponieważ w
obliczeniach wymagane są różnice kolejnych wierszy, początkowo
trzeba obliczyć kilka wierszy za pomocą metody Eulera lub w inny
sposób. Metoda Adamsa jest właściwie najprostszą z metod odpo-
wiednich dla maszyn cyfrowych.
10.3.3. Metoda Runge-Kutta
Omawiane powyżej metody obliczeń numerycznych wymagają
użycia danych liczbowych z poprzednich kroków. Dokładność przy
tym wzrasta w miarę zwiększenia liczby poprzedzających kroków.
Grupa metod, której nadano ogólnie nazwę Runge-Kutta, różni się
pod tym względem od wyżej wymienionych.
Metoda Runge-Kutta charakteryzuje się następującymi zaletami:
M.Luft, Z.A
ukasik
361
1. Jest samostartująca, tzn. nie wymagane są żadne dane z po-
przednich punktów. Metoda może być stosowana nawet do obliczania
pierwszego kroku.
2. Potrzebna jest mniejsza pojemność pamięci operacyjnej maszyny
cyfrowej, gdyż nie są wymagane dane liczbowe z poprzednich kro-
ków.
3. Metoda jest samorzutnie stabilna i może być użyta do obliczeń,
które wymagają dużej liczby kroków.
4. Zmiana kroku może być przeprowadzona w dowolnym punkcie
obliczeń.
Do nielicznych wad tej metody zaliczyć należy nieco trudniejsze
niż w innych metodach oszacowania błędu.
Spośród metod Runge-Kutta różnych rzędów, najszersze zastoso-
wania znalazła metoda czteropunktowa. Rozwiązanie równania
(10.27) według tej metody wyraża się wzorem:
1
xk 1 xk ( K1 2 K 2 K K ) (10.37)
2 3 4
6
przy czym:
K h f (t , x )
1 k k
h K
1
K h f (t , x )
2 k k
2 2
(10.38)
h K
2
K h f (t , x )
3 k k
2 2
K h f (t h, x K )
4 k k 3
tk t0 k h ( k 0,1,2 ,...)
Zastosowanie powyższej metody rozpatrzymy na następującym przy-
kładzie:
Przykład 10.1.
M.Luft, Z.A
ukasik
362
Wyznaczamy funkcję x(t), będącą rozwiązaniem równania różnicz-
kowego
dx
1 x; x (0) 0 (10.39)
dt
W naszym przypadku f(t,x) = 1 - x i wzory (10.38) przyjmują postać:
K h (1 xk )
1
K
K h 1 xk 1
2
2
(10.40)
K
K h 1 xk 2
3
2
K h 1 xk K
4 3
10.3.4. Numeryczne rozwiązywanie układów równań
różniczkowych
Każde zwyczajne równanie różniczkowe n-tego rzędu jest równo-
ważne układowi n równań pierwszego rzędu, jak pokazano w p. 10.2.
Jeżeli każde z powyższych równań pierwszego rzędu zapisać można
w postaci (10.27), to można je rozwiązać numerycznie metodami
omówionymi w poprzednich punktach. Zalety metody Runge-Kutta
przemawiają za jej stosowaniem. Rozpatrzmy zatem sposób rozwią-
zywania układu równań różniczkowych tą właśnie metodą.
Nie tracąc nic z ogólności rozwiązań, rozpatrzmy układ dwóch
równań różniczkowych pierwszego rzędu:
x f (t, x, y )
(10.41)
y g (t, x, y )
których rozwiązaniami są x = x(t) oraz y = y(t).
M.Luft, Z.A
ukasik
363
StosujÄ…c znane wzory czteropunktowej metody Runge-Kutte mamy:
1
x x ( K 2 K 2 K K )
k 1 k 1 2 3 4
6
(10.42)
1
y y (M 2 M 2 M M )
k 1 k 1 2 3 4
6
przy czym
t t0 h k
k
K h f (t , x , y )
1 k k k
M h g (t , x , y )
1 k k k
h K M
1 1
K h f (t , x , y )
2 k k k
2 2 2
h K M
1 1
M h g (t , x , y )
2 k k k
2 2 2
(10.43)
h K M
2 2
K h f (t , x , y )
3 k k k
2 2 2
h K M
2 2
M h g (t , x , y )
3 k k k
2 2 2
K h f (t h, x K , y M )
4 k k 3 k 3
M h g (t h, x K , y M )
4 k k 3 k 3
W każdym kroku obliczania wartości xk i yk należy najpierw kolej-
no wyznaczyć wartości współczynników (10.43), a następnie wstawić
je do wzorów na xk+1 i yk+1.
Prześledz
my to na poniższym przykładzie.
M.Luft, Z.A
ukasik
364
Przykład 10.2.
Należy znalez
ć rozwiązanie równania różniczkowego
x x x (10.44)
przy czym x(0) = 1, x (0) = 0.
Podstawiając y = x mamy zamiast równania różniczkowego drugiego
rzędu, układ dwóch równań rzędu pierwszego.
x y
(10.45)
y y x
gdzie:
x(0) = 1 oraz y(0) = 0.
Stosując wzory (10.37) i (10.38) otrzymujemy wartości
K h y
1 k
M h ( y x )
1 k k
M
1
K h ( y )
2 k
2
M K
1 1
M h ( y x )
2 k k
2 2
(10.46)
M
2
K h ( y )
3 k
2
M K
2 2
M h ( y x )
3 k k
2 2
K h ( y M )
4 k 3
M h ( y M x K )
4 k 3 k 3
M.Luft, Z.A
ukasik
365
które podstawiamy do wzorów:
1
x x ( K 2 K 2 K K )
k 1 k 1 2 3 4
6
(10.47)
1
y y (M 2 M 2 M M )
k 1 k 1 2 3 4
6
Uzyskamy w ten sposób kolejne wartości rozwiązania x(t), jak i rów-
nież jego pochodnej x (t).
10.4. Pytania kontrolne
1. Omówić właściwości idealnego wzmacniacza operacyjnego.
2. Wyprowadzić wzór na transmitancję operatorową wzmacniacza ze
sprzężeniem zwrotnym.
3. Podać modele analogowe podstawowych członów dynamicznych.
4. Zasady modelowania liniowych równań różniczkowych ze skalo-
waniem amplitudy i czasu.
5. Zasady modelowania układów dynamicznych na podstawie danej
transmitancji
6. operatorowej.
7. Zbudować schemat maszynowy układu modelującego równanie
różniczkowe:
x 2 x 0,5 x 0,2 x 2 1(t )
dla x(0) = 0, x (0) = 4, x  (0) = 0.
x 10 x 1000 x 100 1(t )
dla x(0) = x (0) = 0.
8. Podać schemat maszynowy układu modelującego układu o trans-
mitancji operatorowej:
M.Luft, Z.A
ukasik
366
2
s 5 s 0,6 s 0,025
a) G (s) ; b) G (s)
2 2
s 0,5s 0,01 s 0,9 s 0,02
9. Napisać algorytm numerycznego rozwiązywania równania róż-
niczkowego (metodÄ… Runge-Kutta):
y 2 y 100 y 1
a)
y (0) y (0) 0
x 5 x 20 x 2 x 0
b)
x(0) 1, x (0) 6, x (0) 0
3
x 0.2 x x 0
c)
x(0) 1, x (0) 0