DOI 10.21008/j.1897-0737.2017.89.0004
__________________________________________
* Politechnika Opolska.
Bernard BARON*
Joanna KOLAŃSKA-PŁUSKA*
METODY CAŁKOWANIA NUMERYCZNEGO ZAGADNIEŃ SZTYWNYCH W ELEKTROTECHNICE
W pracy przedstawiono testy własnej biblioteki obliczeniowej zawierającej szereg metod rozwiązywania równań różniczkowych zwyczajnych. W zastosowaniu do bada- nia dynamiki układów sztywnych zaprezentowano metody niejawne Rungego-Kutty oraz metody wielokrokowe. Jako przykład testujący opracowanej biblioteki w tym za- kresie zaproponowano układ dwóch równań różniczkowych Kapsa [2]. Przetestowane metody zastosowano także do rozwiązywania równań na przykładzie prostego układu RLC z prostownikiem.
SŁOWA KLUCZOWE: równania różniczkowe, metody niejawne IRK, metoda wielokrokowa Geara
1.WPROWADZENIE
Tradycyjnie do badania dynamiki układów elektrycznych [2], w których wy- stępują równania sztywne stosuje się metody wielokrokowe Geara. W praktyce można stosować także metody niejawne Rungego-Kutty, które mimo większego kosztu numerycznego obliczeń niewiele odbiegają od metod Geara. Ponadto charakteryzują się wysokimi rzędami oraz większą stabilnością [3]. Spośród wielu metod niejawnych Rungego-Kutty zaprezentowano metodę RADAU II rzędu od 3–11. Wykonano testy umożliwiające analizę błędów estymacji tych algorytmów dla równań sztywnych zadanych w postaci normalnej
] t ), t ( dt [
) t (
dX F X
(1)
gdzie: X(t0) = X0 – wektor warunków początkowych, X(t) – N wymiarowy wek- tor w chwili t, F[X(t),t] – jest funkcją wektorową zmiennej wektorowej X(t) i parametru t.
2. METODY NIEJAWNE RUNGEGO-KUTTY (IRK) Metody niejawne typu Runge-Kuta [4] mają ogólnie postać:
) i ( i j m
1 j i 1
i X w K
X
(2)
gdzie wj są stałymi, natomiast wektory Kj(i)
wyrażają się wzorami:
] h c t , a [
h jl (li) i j i
m
1 l i i ) i (
i
K X
F
K dla j = 1,...,m (3)
przy czym hi = ti+1 – ti,
lj m
1 j
l a
c
(4)
Tabela współczynników Butchera dla metody niejawnej ma postać:
Tabela 1. Tabela Butchera [8] dla metody niejawnej
W przypadku metod niejawnych [4] wektory współczynników tworzą nieli- niowy układ równań algebraicznych. Przy rozwiązywaniu układu N równań różniczkowych wektory te mają N składowych więc dla metody m etapowej otrzymujemy układ Nm równań nieliniowych. Fakt rozwiązywania w każdym kroku całkowania układu równań nieliniowych oznacza stosunkowo duży koszt obliczeń numerycznych.
Do rozwiązywania układu równań (1) stosuje się metodę Newtona [5–6, 10].
Rozwiązywanie w procesie iteracyjnym Newtona układu równań liniowych spowalnia niewątpliwie proces całkowania w metodach niejawnych Rungego- Kutty. Dzieje się tak tym bardziej im większy jest układ równań różniczkowych.
W praktyce modelowania dynamiki różnych układów jest bardzo często tak, że im większy układ równań różniczkowych tym rzadsza jest macierz Jacobiego J(F)(X,t) prawej strony równaniu (1) tj. funkcji F(X,t). Żeby więc nie spowodo-
wać wzrostu kosztów obliczeń w tego typu przypadkach należy zastosować techniki macierzy rzadkich.
Dla niejawnych m–etapowych metod Rungego-Kutty istnieją pewne wybory węzłów c1, c2, ..., cm, dla których można otrzymać wysokie rzędy metody. Jak wiadomo, kwadraturą o maksymalnym rzędzie aproksymacji jest kwadratura Gaussa. Dlatego celowym jest taki wybór węzłów c1, c2, ..., cm, które stanowią zera kwadratowej formuły wysokiego rzędu. Kwadratury Gaussa-Legendrea, Radau i Lobatto pozwalają na otrzymanie rzędów metody odpowiednio 2m, 2m–1 i 2m–2.
W algorytmie RADAU IIA węzły c1, c2, ..., cm, stanowią zera wielomianu Legendre’a z podstawieniem 1-2x, co gwarantuje że ostatni węzeł jest końcem przedziału całkowania (cm = 1).
) x 2 1 ( P ) x 2 1 ( P ) x (
P m m 1
) def 1 (
m m = 1,2,...,s (5)
Węzły czasowe cj zadaje się równaniem 0 ) c (
Pm(1) j j = 1,2,...,m (6) co gwarantuje że cm = , a więc ostatnim węzłem czasowym jest węzeł odpowia- dający prawemu końcowi przedziału całkowania.
Rys. 1. Wielomiany Legendre’a w przedziele [0,1]
Macierz Vandermonde'a dla tych współczynników ma postać:
U =
m 1
m 1
m 3 1 m 2 1 m 1
2 m 2
3 2
2 2
1
m 3
2 1
) c ( )
c ( ) c ( ) c (
) c ( )
c ( ) c ( ) c (
c c
c c
1 1
1 1
(7)
Wyznaczając jej macierz odwrotną otrzymuje się macierz V w postaci:
V = U-1 =
mm 1
m 1 m
m 2 22
21
m 1 12
11
V V
V
V V
V
V V
V
(8)
Znajomość macierzy odwrotnej V-1 pozwala zgodnie z teorią metod niejaw- nych Rungego-Kutty na mocy warunku uproszczonego B(m) wyznaczyć współ- czynniki wi tablicy Butchera:
k v 1
w ik
m
1 k
i
i = 1,2,...,m (9)
Wykazuje się, że m–etapowej metody Radau IIA osiąga się rząd 2m–1.
Z warunku uproszczonego C(m) IRK [6] wyznacza się elementy tablicy Bu- tchera aij:
k ) c v ( a
k i jk m
1 k
ik
i,j = 1,2,...,m (10)
które są generowane przez konstruktory klas rozpatrywanej biblioteki.
W bibliotece metod całkowania do estymacji błędu całkowania zastosowano koncepcję de Swarta i Sӧderlinda. Zastąpiono w ten sposób formułę włożoną Hairera [4] formułą niejawną [de Swarta i Sӧderlinda]. Pozwoliło to na realiza- cję w programie automatycznego doboru kroku całkowania z większą gwarancją stabilności rozwiązania.
3. METODY WIELOKROKOWE GEARA
Metoda wielokrokowa Geara ma postać:
) t , ( h b
ai n i 1 0 n 1 n 1
k
1 i 1
n
X F XX (11)
gdzie ai (i = 1, 2, ..., k) oraz b0 są rozwiązaniem układu równań 1
ai
k
1 i
(12)
1 jb a ) 1 i
( j i 0
k
1 i
, j > 0 (13)
W algorytmie Geara jedyną nieznaną wielkością w równaniu (11) jest Xn+1, co można zapisać w postaci:
0 a
) t , ( h b )
( i n 1 1
k
1 i 1 n 1 n 0
1 n 1
n
X F X
XF X (14)
Do rozwiązania równania nieliniowego ze względu na Xn+1, stosuje się itera- cję prostą, jak dla algorytmu Adamsa-Multona. Jak wykazał Gear [1], metoda iteracji prostej nie jest właściwym członem korekcyjnym do wyprowadzenia sztywno stabilnych algorytmów. Ponieważ zasadniczym zyskiem z zastosowania algorytmów sztywno stabilnych jest możliwość przyjęcia dużego kroku po znik- nięciu szybkiej składowej przejściowej, za człon korekcyjny przyjmuje się me- todę Newtona.
W programie zastosowano macierz Nordsiecka [9] do konstrukcji algorytmu Geara, wygodnej ze względu na dobór kroku całkowania i rzędu metody.
3. PRZYKŁAD TESTUJĄCY
Zaprezentowane w punktach 2 i 3 algorytmy przetestowano na przykłado- wym układzie równań różniczkowych sztywnych, dla których istnieje analitycz- ne rozwiązanie. Pozwoli to na oszacowanie maksymalnych błędów rzeczywi- stych popełnianych przez algorytmy typu Geara oraz Radau IIA wysokich rzę- dów oraz porównania ich z błędem estymowanym. W pierwszej kolejności roz- patrywany będzie problem Kapsa stanowiący układ dwóch równań różniczko- wych z małym parametrem sztywności .
Zadany jest sztywny układ dwóch równań różniczkowych nieliniowych o po- staci (15) Kapsa P. (1981):
3 2 t 1
1 e (x (t))
) t ( x ) 1 2 dt (
) t ( dx
2 2 2 1 t 2
2 e (x (t)) (x (t)) dt
) t (
dx (15)
0 t 6, 0 < << 1
Układ powyższych równań różniczkowych przy warunkach początkowych x1 (0) = x2(0) = 1 ma następujące rozwiązanie analityczne:
t 2 1(t) e
x , x2(t)et (16) W konstrukcji układu równań mały parametr (np. = 1,0e-8) silnie wpływa na sztywność tego układu.
Znajomość rozwiązania analitycznego równań Kapsa pozwala na badanie błędów rzeczywistych lokalnych i globalnych testowanych algorytmów dowol- nych rzędów. Umożliwia także badanie skuteczności estymacji błędów lokal- nych na tle błędów rzeczywistych metody.
Błąd rzeczywisty lokalny dla dowolnej k-tej zmiennej stanu xk(t) definiuje się następująco:
dla x (t) 1
) t ( x
) t ( x ) h , t ( x
1 ) t ( x dla ) t ( x ) h , t ( x ) h , t ( e
k k
k k
k k
def k
k (17)
e (t,h)
max ) h , t (
e k
N k 1
gdzie xk(t,h) oraz xk(t) są rozwiązaniami odpowiednio numerycznymi i anali- tycznymi dla danych problemów. Globalny błąd rzeczywisty definiuje się nastę- pująco:
)}
h , t ( e { max ) h (
E i
i
ti t0 ih (18)
Interfejs użytkownika projektu realizującego badanie równań Kapsa przed- stawiono na rys. 2.
Rys. 2. Interfejs projektu testującego metody IRK
Na przykładzie układu równań Kapsa przeprowadzono również badania te- stowe opracowanych algorytmów w zależności od parametru sztywności tego układu przy automatycznym doborze kroku całkowania (tabela 2). Z badań tych wynika, że przy zmianie parametru równań + 1,0e-1 1,0e-8 ilość iteracji dla metod Radau IIA47 utrzymuje się mnie więcej na tym samym poziomie od 10 do 20 iteracji przy podobnych czasach obliczeń.
Rys. 3. Błędy rzeczywiste i estymowane metody niejawnej Rungego-Kutty z formułą Radau IIA dla różnych rzędów z parametrami estymacji w0 = 0,01 oraz w0/ = 0,067 funkcją kroku całkowa- nia h oraz funkcją czasu t przy automatycznym doborze kroku całkowania dla zadanego błędu
względnego i absolutnego = 1,0e-8 dla estymacji błędu wg. de Swarta-Sӧderlinda Tabela 2. Wyniki testu metod IRK całkowania w zależności od parametru sztywności równania Kapsa w postaci ilości kroków i czasu obliczeń przy automatycznym doborze kroków całkowania z estymacją błędu wg. de Swarta, Sӧderlinda dla zadanego błędu względnego i absolutnego = 1e-6
Wykonany test wykazał, że metoda wielokrokowa Geara wykonuje pięcio- krotnie większą liczbę kroków całkowania w porównaniu z metodami Radau IIA.
IRK Radau IIA 4 etapowy rzędu 7
Metoda Geara Automatyczny dobór
kroku i rzędu
Ilość iteracji
Czas [ms]
całkowania
Ilość iteracji
Czas [ms]
całkowania
10-1 9 322 65 18
10-2 9 199 55 18
10-3 13 460 68 19
10-4 16 209 65 19
10-5 8 368 59 19
10-6 19 387 59 18
10-7 19 85 59 17
10-8 20 79 59 18
10-9 21 209 59 17
10-10 22 314 59 17
Ze względu na rozwiązywanie znacznie mniejszego układu równań nielinio- wych w każdym kroku całkowania metoda Geara pomimo wykonania większej liczby kroków całkowania, wymaga kilkukrotnie mniejszego czasu obliczeń.
W przykładzie symulacyjnym z elektrotechniki wybrano układ prostownika pokazany na rys. 4. Do obliczeń w programie C# [7] zaimplementowano dwie metody: RADAU IIA oraz Geara.
Rys. 4. Układ prostownika
Równania opisujące układ prostownika:
) t ( u )]
t ( i [ u ) t ( u ) t ( u ) t (
e Rzr Lzr C (19)
)]
t ( u )]
t ( i [ u ) t ( i R ) t ( e L [
1 dt
) t ( di
C zr
zr
(20)
R ] ) t ( ) u t ( i C[
1 dt
) t ( du
0 C
C (21)
o parametrach: Rźródła = 0,5 [], Rprzewodzenia = 0,5 [],Rzaporowe = 10 500 [k], Robciążenia = 100 [], L = 5 [mH], C = 0,1 [F] oraz wymuszeniu napięciowym Emax = 100 [V] i częstotliwości f = 50 [Hz].
Obliczenia zrealizowano dla dwóch metod całkowania tj. niejawnej metody Radau IIA 4 i 5 etapowej oraz metody wielokrokowej Geara z automatycznym doborem kroku całkowania 4–6 rzędu przy założonym błędzie względnym i absolutnym 1,0e-7.
Wyniki tej symulacji pokazują, że ze wzrostem rezystancji obciążenia pro- blem rozwiązywania równań różniczkowych staje się coraz to bardziej sztywny.
Tabela 3. Porównanie ilości iteracji oraz czasów obliczeń dla metody niejawnej RADAU IIA i wielokrokowej Gear’a
Metoda RADAU IIA Wielokrokowa Geara
Rzap [] Niteracji tmax[ms] Niteracji tmax[ms]
10000 428 1750 5628 375
50000 420 854 176455 11232
100000 383 1020 175952 11167
Rys. 5. Interfejs projektu do badania stanów dynamicznych układu prostowniczego
Ze wzrostem rezystancji zaporowej układu wzrasta sztywność równań opisu- jących dynamikę układu. Przeprowadzone symulacje wykazały, że wraz ze wzrostem rezystancji zaporowej liczba iteracji w metodzie wielokrokowej Geara znacząco się zwiększa, natomiast czas obliczeń jest rzędu 11 s. W przypadku tego zadania bardziej skuteczna wydaje się metoda RADAU IIA 5 rzędu w któ- rej czas obliczeń jest dziesięciokrotnie mniejszy z zachowaniem takiej samej dokładności jak metoda Geara.
4. PODSUMOWANIE
W pracy przedstawiono wybrane dwie metody niejawne rozwiązywania rów- nań różniczkowych zwyczajnych spośród wielu metod zawartych w zorganizo- wanej obiektowo bibliotece numerycznej [11]. Bibliotekę wyposażono w szereg metod związanych bezpośrednio realizacją rozwiązania ze stałym lub zmiennym krokiem całkowania w oparciu o estymowany błąd ale również w metody doty- czące aproksymacji rozwiązania wielomianami umożliwiające pobór rozwiąza- nia w dowolnej chwili czasu z przedziału całkowania. Ponadto oprogramowanie wyposażono w metody umożliwiające śledzenie estymowanych błędów całko- wania zarówno lokalnych jak i globalnych co nie tylko umożliwiło badania te- stujące algorytmy lecz również jest pomocne przy wyborze najlepszej metody całkowania. Organizacja obiektowa biblioteki umożliwia ponadto tworzenie programów dedykowanych dla konkretnych zagadnień badania dynamiki ukła- dów.
LITERATURA
[1] Gear C.W. Numerical Initial Value Problems in Ordinary Differential Equations. Englewood Cliffs; Prentice–Hall, Inc., New York, 1971.
[2] Krupowicz A.: Metody numeryczne zagadnień początkowych równań różniczko- wych zwyczajnych. PWN, Warszawa 1984.
[3] Rattenbury N.: Almost Runge-Kutta methods for stiff and non–stiff problems, PhD Thesis, The University of Auckland, New Zealand, 2005
[4] Hairer E. ,Wanner G.: Solving Ordinary Differential Equations II: Stiff and Differential–algebraic Problems, 2nd revised ed., Springer, Berlin, 1996.
[5] Kennedy C. A. and Carpenter M. H.: Additive Runge-Kutta Schemes for Convection–Diffusion–Reaction Equations. Technical report, NASA, 2001.
NASA/TM–2001–211038.
[6] Dekker K., Verwer J.G.: Stability of Runge-Kutta methods for stiff nonlinear differential equations. Elsevier Science Publishers B. V., North-Holland Amsterdam–New York – Oxford 1984.
[7] Baron B. Metody numeryczne w C++. Wydawnictwo Helion, Gliwice 2004.
[8] Butcher J. C., and Chen D. J. L.: A new type of singlyimplicit Runge-Kutta method, Applied Numerical Mathematics, (2000), 34: 179–188.
[9] Kvaerno K.: Singly Diagonally Implicit Runge-Kutta Methods with an Explicit First Stage. BIT Numerical Mathematics, 2004.
[10] de Swart J. J.B., Sӧderlind G.: On the construction of error estimators for implicit Runge-Kutta methods, Journal of Computational and Applied Mathematics 86 (1997) 347–35.
[11] Baron B., Kolańska–Płuska J.: Metody numeryczne rozwiązywania równań róż- niczkowych zwyczajnych w języku C#, Oficyna Wydawnicza Politechniki Opol- skiej, 2015.
METHODS OF NUMERICAL INTEGRATION FOR STIFF PROBLEMS IN ELECTRICAL ENGINEERING
The paper presents selected two methods of solving implicit ordinary differential equations from a variety of methods defined in a object structured numerical library.
Library is equipped with a number of methods directly related to the implementation of solutions with constant or variable step integration on the basis of estimated error.
Additionary library includes methods for approximation of polynomial solutions which enable power solutions at any point in time with the integration interval. The software also includes methods for tracking the estimated integrating error both local and global which not only enabled the tests of the algorithms but also is helpful in choosing the best method of integration. Object organization library allows creation of programs dedicated to specific subjects research systems dynamics, as well.
(Received: 28. 02. 2017, revised: 4. 03. 2017)
