WYDZIAŁ ETI PG

Katedra Systemó w Elektróniki Mórskiej

Technika obliczeniowa i symulacyjna

- laboratorium

Numeryczne rozwiązywanie równań liniowych

MATLAB i SPICE jako narzędzia do obliczania

prądów, napięć i mocy w obwodzie elektrycznym

Odsłona 1:

Rozwiązywanie – w MATLAB-ie – równań opisujących obwód. Przypadek

oczkowego opisu obwodu liniowego i stałych pobudzeń (prądy, napięcia, moce na ele-

mentach, źródła zastępcze Thevenina i Nortona). Micro-Cap – analiza Dynamic DC

(Materiały dó ćwiczenia nr 2)

Opracówał: Witóld Szkudlin ski

Adaptacja: Czesław Stefański

Gdan sk 2016

2

1. Wstęp

W prostym obwodzie elektrycznym, na przykład składającym się z podobnej liczby elementów jak w ob-

wodzie przedstawionym na rys. 1a, poszukiwanie takich podstawowych wielkości elektrycznych jak natężenie
prądu (prąd) 𝒊, czy też różnica potencjałów (napięcie) 𝒖, może być względnie łatwym zadaniem.

Uznajemy, że dla obwodów z rys. 1 jest to szczególnie łatwe również z tego względu, że mamy w tym

przypadku do czynienia z tak zwanymi liniowymi obwodami prądu stałego.

i

0

i

1

i

2

i

3

E

0

R

1

R

2

R

3

u

0

u

2

u

3

u

1

J

4

u

4

i

4

i

0

i

1

i

2

i

3

E

0

R

1

R

2

R

3

u

0

u

2

u

3

u

1

Rys. 1a. Prosty obwód prądu stałego

Rys. 1b. Przykład strzałkowania skojarzonego

J

4

u

4

i

4

i

0

i

1

i

2

i

3

E

0

R

1

R

2

R

3

u

0

u

2

u

3

u

1

Rys. 1c. Przykład strzałkowania „dwuskojarzonego” (źródłowo-odbiorczego

1

)

Dla obwodu z rys. 1a, opierając się na regułach analitycznych opracowanych przez żyjących w XIX wieku

Georga Ohma, Wernera von Siemensa oraz Gustawa Kirchhoffa otrzymujemy prądy oraz napięcia:

𝑖

1

=

𝐸

𝑅

1

+

1

𝐺

2

+𝐺

3

,

𝑖

2

=

𝐺

2

𝐺

2

+ 𝐺

3

⋅ 𝑖

1

,

𝑖

3

=

𝐺

3

𝐺

2

+ 𝐺

3

⋅ 𝑖

1

,

𝑖

0

= −𝑖

1

,

(1)

𝑢

1

= 𝑅

1

⋅ 𝑖

1

,

𝑢

2

= 𝑅

2

⋅ 𝑖

2

,

𝑢

3

= 𝑅

3

⋅ 𝑖

3

,

𝑢

0

= 𝐸,

gdzie

𝐺

2

=

1

𝑅

2

i 𝐺

3

=

1

𝑅

3

.

Do wzorów (1) trzeba dodać, że zwrot spadku napięcia na każdym z elementów obwodu z rys. 1a jest

przeciwny do kierunku przepływającego przez ten element prądu. Sytuacja taka ma miejsce, gdy stosujemy
tzw. strzałkowanie skojarzone prądów i napięć na dwójnikach obwodu (na rys. 1b też mamy strzałkowanie
skojarzone).

Gdy na wszystkich dwójnikach obwodu, za wyjątkiem wszystkich niezależnych źródeł napięciowych i prą-

dowych, zastosowano strzałkowanie skojarzone, będziemy mówić o strzałkowaniu „dwuskojarzonym”. Nie-
którzy autorzy wolą mówić o strzałkowaniu źródłowo-odbiorczym

1

. Taką sytuację mamy na rysunku 1c.

Bywa, że rachując na papierze, strzałkujemy „oszczędnie”, to znaczy dla zmniejszenia liczby niewiado-

mych rezygnujemy z wielokrotnego oznaczania w istocie tego samego prądu, czy napięcia, np. dla rysunków
1a, b, c na źródłach niezależnych zrezygnowalibyśmy z oznaczania dodatkowo prądu 𝑖

0

źródła napięciowego

𝐸

0

i odpowiednio napięcia 𝑢

4

źródła prądowego 𝐽

4

, czy też z oznaczania jednego z napięć – 𝑢

2

, 𝑢

3

albo 𝑢

4

(np. gdyby na rysunku 1c nie oznaczono prądu 𝑖

0

, to jako prąd źródła 𝐸

0

można byłoby traktować prąd 𝑖

1

).

1

Stanisław Bolkowski, Elektrotechnika, WSiP, Warszawa 2005

3

Przy tworzeniu algorytmów analizy automatycznej (komputerowej) wygodnie jest trzymać się zasady, że

każdy dwójnik ma „swój” prąd i „swoje” napięcie, przeciwnie względem siebie zastrzałkowane (przypomi-
namy, że jest to tzw. strzałkowanie skojarzone). Przez to algorytmy stają się bardziej jednorodne i dzięki temu
zwykle są szybsze mimo większej niż przy strzałkowaniu „oszczędnym” liczby niewiadomych.

Ze wzrostem liczby i różnorodności elementów tworzących strukturę obwodu stopniowo zmniejsza się

możliwość – czy też opłacalność – ręcznego (analitycznego) poszukiwania prądów i napięć gałęziowych.

Omawiane ćwiczenie jest poświęcone wdrożeniu umiejętności obliczania prądów i napięć w obwodzie

przy wykorzystaniu narzędzia MATLAB z jego algorytmami numerycznego rozwiązywania układów równań
liniowych.

Ponadto ćwiczący nabiera obycia w korzystaniu z symulacyjnego narzędzia Micro-Cap (jedna z „odmian”

SPICE’a). Narzędzie to, po „narysowaniu” obwodu w jego oknie graficznym przez użytkownika, potrafi auto-
matycznie, acz bez ujawniania ich użytkownikowi, utworzyć odpowiednie równania opisujące obwód, roz-
wiązać je (podobnymi jak MATLAB algorytmami) i udostępnić rozwiązania między innymi we wspomnianym
oknie z obwodem.

2. Poszukiwanie prądów, napięć i mocy w rozbudowanym obwodzie elektrycznym

Podczas tych (i następnych) zajęć obiektem naszej uwagi będzie obwód prądu stałego przedstawiony na

rysunku 2a. Naszym zadaniem będzie dokonać w miarę wszechstronnej analizy tego obwodu, to jest wyzna-
czyć prądy i napięcia, a także moce na każdym z jego elementów oraz porównać wyniki naszej analizy z wy-
nikami symulacji przeprowadzonej narzędziem Micro-Cap, a także sprawdzić bilansowanie się mocy.

Zacznijmy od tego, że chociaż w strukturze badanego obwodu mamy łącznie 11 elementów, to tylko

oznaczone symbolem 𝐸

1

napięciowe źródło niezależne albo niezależne źródło prądu 𝐽

3

jest niezbędne, aby

uzyskać przepływ jakiegokolwiek prądu. Usunięcie obu tych elementów lub jednoczesne wyzerowanie 𝐸

1

i 𝐽

3

implikuje wyzerowanie wszystkich prądów i napięć w obwodzie.

E

1

R

1

R

2

R

3

J

3

R

6

R

4

R

7

I

9

R

8

R

5

U

1

𝑈

1

I

3

I

6

I

5

I

8

I

7

I

2

U

2

U

8

U

5

Iin

ZPSN

Uin

ZPSN

U

7

U

4

I

4



k

l



I

o3

I

o1

I

o4

I

o2

𝐼

9

=

𝑱

𝟔

=

𝐼

𝑜𝑢

𝑡

𝑍𝑃

𝑆𝑁

=

𝑔

⋅𝑈

𝑖𝑛

𝑍𝑁

𝑆𝑃

𝑈

𝑖𝑛

𝑍𝑁

𝑆𝑃

=

𝑈

5

I

1

J

6

𝐼

3

𝐼

1

𝐼

5

𝐼

6

U

6

U

9

U

3

Rys. 2a. Badany obwód prądu stałego

4

Oprócz źródeł niezależnych 𝐸

1

i 𝐽

3

jest w obwodzie zależne (sterowane) źródło wytwarzające prąd 𝐼

9

liniowo zależny od napięcia 𝑈

5

(czyli ZPSN – źródło prądowe 𝐽

6

sterowane napięciem 𝑈

5

).

Na rysunku 2a zaznaczono prądy gałęziowe (𝐼

1

, 𝐼

2

, … , 𝐼

9

) i na początek niech to one stanowią poszuki-

wane wielkości obwodowe. Ogólnie rzecz biorąc, aby znaleźć 9 niewiadomych trzeba ułożyć 9 równań, jednak
w tym przypadku, w oparciu o koncepcję tzw. prądów oczkowych (𝐼

𝑜1

, 𝐼

𝑜2

, 𝐼

𝑜3

, 𝐼

𝑜4

), potrzebujemy tych rów-

nań tylko 4.

Podobnie w oparciu o koncepcję tzw. potencjałów węzłowych możemy też, w tym przypadku, ułożyć

tylko cztery równania z niewiadomymi potencjałami 𝑉

1

, 𝑉

2

, 𝑉

3

, 𝑉

4

i za ich pomocą wyznaczyć wszystkie napię-

cia i prądy na elementach obwodu. To będzie stanowić materiał na kolejne zajęcia. Ale po kolei!

Metoda prądów oczkowych
Klasyczna, niezmodyfikowana metoda prądów oczkowych „nienawidzi” sytuacji, gdy w obwodzie wystę-

pują inne elementy, niż w tym zdaniu wymienione, czyli: niezależne źródła napięciowe, źródła napięciowe
sterowane prądem oraz rezystory.

Mając powyższe na uwadze przekształcamy niezależne rzeczywiste źródło prądowe {𝑱

𝟑

‖ 𝑹

𝟑

} w rzeczy-

wiste źródło napięciowe {𝑬

𝟑

=𝐽

3

𝑅

3

¦ 𝑹

𝟑

} oraz rzeczywiste źródło prądowe sterowane napięciem

{𝑱

𝟔

=𝑔𝑼

𝟓

‖ 𝑹

𝟔

} w rzeczywiste źródło napięciowe sterowane napięciem {𝑬

𝟔

=𝑔𝑅

6

𝑼

𝟓

¦ 𝑹

𝟔

}, a to z kolei już

łatwiutko w rzeczywiste źródło napięciowe sterowane prądem {𝑬

𝟔

=𝑔𝑅

5

𝑅

6

𝑰

𝟓

¦ 𝑹

𝟔

}, bo 𝑈

5

= 𝑅

5

⋅ 𝐼

5

.

Obwód po przekształceniach pokazano na rys. 2b.

E

1

R

1

R

2

R

3

R

6

R

4

R

7

E

6

=gR

5

R

6

I

5

=gR

5

R

6

(I

o1

-I

o2

)

R

8

R

5

𝐼

3

𝐼

6

𝐼

5

𝐼

1

U

1

𝑈

1

I

5

I

8

I

7

I

2

U

2

U

8

U

5

U

7

U

4

I

4



k

l



I

o3

I

o1

I

o4

I

o2

I

1

E

3

=R

3

J

3

U

6

(=U

9

)

U

3

Rys. 2b. Badany obwód prądu stałego dostosowany do metody prądów oczkowych

Dla obwodu z rysunku 2b algorytm układania równań według metody prądów oczkowych prowadzi do

następującego układu równań

5





























































































































































2

6

5

1

6

5

3

3

1

6

3

1

4

3

2

1

7

6

4

4

7

4

4

3

2

2

7

8

7

5

5

2

5

5

2

1

0

0

0

0

0

0

o

o

o

o

o

o

I

R

gR

I

R

gR

R

J

E

E

E

E

I

I

I

I

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

(2)

Po przesunięciu z prawej na lewą stronę składowych z prądami oczkowymi otrzymujemy



































































































































0

0

0

0

0

3

3

1

4

3

2

1

7

6

4

4

5

6

7

5

6

4

4

3

2

2

7

8

7

5

5

2

5

5

2

1

R

J

E

I

I

I

I

R

R

R

R

R

gR

R

R

gR

R

R

R

R

R

R

R

R

R

R

R

R

R

R

R

o

o

o

o

(3)

Wracając do prądów gałęziowych (𝐼

1

, 𝐼

2

, . . . , 𝐼

9

), można je wyrazić przez otrzymane z równania (3) prądy

oczkowe. Mamy bowiem (w zgodzie z rysunkami 2a i 2b):

𝐼

1

= −𝐼

𝑜1

,

𝐼

2

= 𝐼

𝑜3

− 𝐼

𝑜1

,

𝐼

3

− 𝐽

3

= 𝐼

3

= 𝐼

𝑜3

,

𝐼

4

= 𝐼

𝑜3

− 𝐼

𝑜4

,

𝐼

5

= 𝐼

𝑜1

− 𝐼

𝑜2

,

𝐼

6

+ 𝐼

9

= 𝐼

6

= 𝐼

𝑜4

,

𝐼

7

= 𝐼

𝑜4

− 𝐼

𝑜2

,

𝐼

8

= 𝐼

𝑜2

,

𝐼

9

= 𝑱

𝟔

= 𝑔𝑅

5

𝑰

𝟓

.

Organizując obliczenia na podstawie powyższych wzorów i równania (3), na przykład w MATLAB-ie, po-

winniśmy trzymać się zasady, że do programu kierowane są dane nieprzetworzone, a obliczenia są realizo-
wane w programie. Zatem w tym przypadku jako dane do programu powinny posłużyć wartości parametrów
elementów obwodu, na przykład zadane jako R1=…, R2=…, …, g=…, E1=…, J3=… . Instrukcje obli-
czeniowe natomiast powinny być realizacją zależności (3) oraz po niej następujących zależności na prądy.

Bilans mocy w obwodzie
Bilans mocy w obwodzie układamy przyjmując, że iloczyn prądu i napięcia przy przeciwnej orientacji

napięcia na dwójniku i płynącego przez ten dwójnik prądu oznacza moc traconą na tym dwójniku; zawsze na
rezystorach moc tak liczona jest dodatnia, zaś bardzo często liczona w ten sam sposób moc na źródłach jest
ujemna; ujemną moc traconą należy interpretować (po zmianie znaku) jako dodatnią moc dostarczoną. Dla
przykładowego obwodu z rysunku 2a mamy

"

𝑃

rezystorów

tracona

= ∑ 𝐼

𝑘

⋅ 𝑈

𝑘

8

𝑘=1

= ∑ 𝑅

𝑘

⋅ 𝐼

𝑘

2

8

𝑘=1

.

(10a)

Moc dostarczana przez źródła (jest równa minus mocy traconej przez źródła):

"

𝑃

źródeł

dostarczana

= −(𝐼

1

⋅ 𝐸

1

− 𝐽

3

⋅ 𝑈

3

+ 𝐼

9

⋅ 𝑈

9

) = −(𝐼

1

⋅ 𝐸

1

− 𝐽

3

⋅ 𝐼

3

⋅ 𝑅

3

+ 𝐼

9

⋅ 𝐼

6

⋅ 𝑅

6

) .

(10b)

Bilans mocy winien być oczywiście zerowy (moc dostarczona – moc tracona = zero):

"

𝑃

źródeł

dostarczana

−

"

𝑃

rezystorów

tracona

= 0.

(11)

Koncepcja źródeł zastępczych Thevenina i Nortona
Jeżeli nasze potrzeby ograniczają się do znalezienia prądu (napięcia, mocy) tylko w kilku zgrupowanych

elementach (w skrajnym przypadku w jednym elemencie) obwodu, to całą pozostałą część struktury możemy
zastąpić tzw. źródłem Thevenina lub źródłem Nortona. Przyjmując, że wybranym ważnym elementem jest
rezystor R

5

, możemy zgodnie z powyższą koncepcją wieloelementowy obwód zredukować do jednej z przed-

stawionych na rysunku 3a,b

postaci.

N

R

5

R

5

I

)

b

T

E

T

R

5

R

5

I

)

a

5

U

5

U

(y)

(x)

Obwód liniowy

zastępowany

źródłem

Thevenina

lub Nortona

q

-t

y

w

ar

ia

n

t

o

b

ci

ąż

en

ia

o

b

w

o

d

u

, q

=1

lu

b

2

i

k

k

i

J

N

-I

5

E

T

-U

5

J

N

)

c

q

U

q

I

6

Rys. 3. Obwód zastępczy Thevenina a), obwód zastępczy Nortona b), schemat ogólny do wyznaczania parametrów obwodu za-

stępczego c).

Poza 𝑅

5

, parametry pozostałych elementów z rysunków 3a,b muszą zostać dopiero obliczone na tej pod-

stawie, że przez rezystor 𝑅

5

winien przepływać taki sam prąd 𝐼

5

, jaki przepływał przez ten element w ramach

obwodu z rysunku 2a.

Wprowadzone przez Thevenina oraz Nortona definicje dają nam następujące zapisy, które odnoszą się

do rysunku 2a i rysunków 3a,b:

𝐸

𝑇

= 𝑈

5

|

𝑅

5

→∞

, 𝐽

𝑁

= 𝐼

5

|

𝑅

5

=0

, 𝑅

𝑇

= 𝑅

𝑁

=

𝐸

𝑇

𝐽

𝑁

.

(12a)

Inaczej sprawę ujmując, prąd Nortona otrzymamy na podstawie rozwiązania równań (3) po wyzerowa-

niu oporności R

5

w macierzy współczynników z (3). W tej sytuacji wymaga to tylko zmiany wartości elementu,

bez zmiany kształtu równania (3).

Napięcie Thevenina otrzymujemy wtedy, kiedy w obwodzie z rysunków 2 usuniemy element R

5

(odpo-

wiada to przyjęciu, że 𝑅

5

→ ∞). To niestety wymaga przekonstruowania równań typu (3) opisujących metodą

oczkową obwód (bowiem, na skutek usunięcia opornika 𝑅

5

, liczba oczek zmniejszyła się o jedno). Natomiast

napięcie 𝑈

5

, gdy już nie ma opornika 𝑅

5

, musimy składać jako sumę napięć na innych elementach.

Na szczęście do wyznaczenia parametrów źródeł Thevenina i Nortona można podejść bardziej jednolicie,

mianowicie parametry te można obliczyć z ogólniejszych zależności, przy oznaczeniach z rysunku 3c:

𝐸

𝑇

=

𝐼

2

⋅ 𝑈

1

− 𝐼

1

⋅ 𝑈

2

𝐼

2

− 𝐼

1

, 𝐽

𝑁

=

𝑈

2

⋅ 𝐼

1

− 𝑈

1

⋅ 𝐼

2

𝑈

2

− 𝑈

1

, 𝑅

𝑇

= 𝑅

𝑁

=

𝑈

1

− 𝑈

2

𝐼

2

− 𝐼

1

(12b)

Przyjmując, że pierwszy wariant obciążenia obwodu to 𝑅

5

= 𝑅

5

1

, a drugi to 𝑅

5

= 𝑅

5

2

oraz uwzględnia-

jąc, że 𝑈

=

𝑞

𝑅

5

𝑞

⋅ 𝐼

𝑞

możemy dla obwodu z rysunków 2 przepisać wzory (12b) w postaci:

𝐸

𝑇

=

𝐼

2

⋅ 𝐼

1

𝐼

2

− 𝐼

1

⋅ ( 𝑅

5

1

− 𝑅

5

2

) , 𝐽

𝑁

=

𝐼

2

⋅ 𝐼

1

𝑅

5

2

⋅ 𝐼

2

− 𝑅

5

1

⋅ 𝐼

1

⋅ ( 𝑅

5

2

− 𝑅

5

1

), 𝑅

𝑇

= 𝑅

𝑁

=

𝑅

5

1

⋅ 𝐼

1

− 𝑅

5

2

⋅ 𝐼

2

𝐼

2

− 𝐼

1

. (12c)

Gdyby dodatkowo założyć, że 𝑅

5

2

= 0, co odpowiada obciążeniu obwodu zwarciem, to uzyskalibyśmy,

że:

𝐸

𝑇

=

𝐼

2

⋅ 𝐼

1

𝐼

2

− 𝐼

1

⋅ 𝑅

5

1

, 𝐽

𝑁

= 𝐼

2

, 𝑅

𝑇

= 𝑅

𝑁

=

𝐼

1

𝐼

2

− 𝐼

1

⋅ 𝑅

5

1

.

(12d)

Oznacza to, że aby wyznaczyć parametry źródeł zastępczych wystarczy policzyć lub pomierzyć prąd

𝑞

𝐼

(𝑞 = 1, 2) płynący przez obciążenie w warunkach dwóch znanych i różnych obciążeń i skorzystać ze wzoru
(12c) (albo (12d), gdy drugie z obciążeń było zwarciem).

Do wyznaczenia źródeł Thevenina lub Nortona, gdy korzystamy z wzorów (12b-c) wystarczy manipulo-

wać obciążeniem (byleby było niezerowe i skończone), a nie strukturą całego obwodu. Wzór (12d) też ma te
zalety, ale tylko przy analizie metodą oczkową (jak się później okaże – w metodzie węzłowej przyjęcie, że
oporność obciążenia wynosi zero – powoduje zmianę struktury obwodu, bo – z punktu widzenia metody –
„skleja” ze sobą dwa węzły w jeden).

Zanim skupimy się na rozwiązywaniu równań w MATLAB-ie policzymy jeszcze moce w obwodach 3a-b.

Moce (liczone jako tracone (!!!)) w obwodzie z rysunku 3a

𝑃

rezystorów

tracona

= (𝑅

𝑇

+ 𝑅

5

) ⋅ 𝐼

5

2

,

(13)

𝑃

źródła napięciowego

tracona

= −𝐼

5

⋅ 𝐸

𝑇

.

Moce (liczone jako tracone (!!!)) w obwodzie z rysunku 3b

𝑃

rezystorów

tracona

= 𝑅

𝑁

⋅ (𝐽

𝑁

− 𝐼

5

)

2

+ 𝑅

5

⋅ 𝐼

5

2

,

(14)

𝑃

źródła prądowego

tracona

= −

𝐽

𝑁

2

𝐺

𝑁

+𝐺

5

= −

𝐽

𝑁

2

𝑅

𝑁

𝑅

5

𝑅

𝑁

+𝑅

5

.

Warto zastanowić się, czy (i dlaczego) suma mocy traconych na wszystkich, oprócz opornika 𝑅

5

, elemen-

tach obwodu 2a jest równa sumie mocy traconych na wszystkich, oprócz opornika 𝑅

5

, elementach obwodu

z rysunku 3a, jak i 3b.

3. Wspomagane programem MATLAB rozwiązywanie układu równań liniowych

7

Zajmiemy się teraz szczegółami numerycznego poszukiwania rozwiązania układu równań liniowych na

przykładzie poniższego układu równań zamieszczonego w materiałach pomocniczych

1

do wykładu

MATLAB

-a

[

2 −3

4

5

4

4 −8

2

3 −3 −4

2

2 −2

3 −5

] ⋅ [

𝑥

1

𝑥

2

𝑥

3

𝑥

4

] = [

28

−4

−7

−13

] .

(15)

W uogólnionym zapisie (przy notacji bez nawiasów kwadratowych) równanie to wygląda następująco

𝑨 ∙ 𝒙 = 𝒃

(16)

W tym przypadku rozwiązanie jest w naszym zasięgu nawet bez wspomagania komputerowego, przy

czym często korzystamy z przejrzystych zapisów metody Cramera

𝑥

1

=

Δ

1

Δ

, 𝑥

2

=

Δ

2

Δ

, 𝑥

3

=

Δ

3

Δ

, 𝑥

4

=

Δ

4

Δ

, Δ ≠ 0 .

(17)

gdzie Δ jest wyznacznikiem macierzy 𝑨, zaś każdy z Δ

𝑛

to wyznacznik macierzy powstałej z 𝑨 przez zastąpienie

w niej wektorem b kolumny o numerze „𝑛-tym”.

Rozwiązanie analityczne (a więc dokładne) równania (15) metodą Cramera daje wyniki

𝑥

1

=

1450
1450

= 1, 𝑥

2

=

2900
1450

= 2, 𝑥

3

=

4350
1450

= 3, 𝑥

4

=

5800
1450

= 4

(18)

Współcześnie, spośród wielu opracowanych metod numerycznego rozwiązywania układów równań li-

niowych, duże znaczenie mają różne warianty metody eliminacji Gaussa. Algorytm tej metody prześledzimy
po uogólnieniu zapisu równania (15). Otrzymujemy

[

𝑎

11

𝑎

12

𝑎

13

… . 𝑎

1𝑁

𝑎

21

𝑎

22

𝑎

23

… . 𝑎

2𝑁

𝑎

31

𝑎

32

𝑎

33

… 𝑎

3𝑁

…

…

…

…

…

𝑎

𝑁1

𝑎

𝑁2

𝑎

𝑁3

… 𝑎

𝑁𝑁

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

…

𝑥

𝑁

]

=

[

𝑏

1

𝑏

2

𝑏

3

…

𝑏

𝑁

]

(19)

Eliminacja niewiadomych w wyrażeniu (19) przebiega następująco.

 W pierwszym kroku od wierszy 2, 3, … . , 𝑁 odejmujemy wiersz pierwszy pomnożony przez taki

współczynnik, który zapewnia jako wynik odejmowania nowy wiersz nie zawierający już niewia-
domej 𝑥

1

.

Po pierwszym kroku (krok 𝑠 = 1; ang.: step) otrzymujemy

[

𝑎

11

𝑎

12

𝑎

13

… 𝑎

1𝑁

0 𝑎

22

(1)

𝑎

23

(1)

… 𝑎

2𝑁

(1)

0 𝑎

32

(1)

𝑎

33

(1)

… 𝑎

3𝑁

(1)

… …

…

… …

0 𝑎

𝑁2

(1)

𝑎

𝑁3

(1)

… 𝑎

𝑁𝑁

(1)

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

…

𝑥

𝑁

]

=

[

𝑏

1

𝑏

2

(1)

𝑏

3

(1)

…

𝑏

𝑁

(1)

]

(20a)

gdzie

𝑎

𝑖𝑘

(1)

= 𝑎

𝑖𝑘

− 𝑙

𝑖1

⋅ 𝑎

1𝑘

,

𝑏

𝑖

(1)

= 𝑏

𝑖

− 𝑙

𝑖1

⋅ 𝑏

1

, przy 𝑙

𝑖1

=

𝑎

𝑖1

𝑎

11

oraz 𝑖, 𝑘 = 2, 3, … , 𝑁

 W drugim kroku od wierszy 3, 4, … . , 𝑁 odejmujemy wiersz drugi pomnożony przez taki współ-

czynnik, który zapewnia jako wynik odejmowania nowy wiersz nie zawierający już niewiadomej
𝑥

2

.

Po drugim kroku (krok 𝑠 = 2) otrzymujemy

[

𝑎

11

𝑎

12

𝑎

13

… 𝑎

1𝑁

0 𝑎

22

(1)

𝑎

23

(1)

… 𝑎

2𝑁

(1)

0 0

𝑎

33

(2)

… 𝑎

3𝑁

(2)

… …

…

… …

0 0

𝑎

𝑁3

(2)

… 𝑎

𝑁𝑁

(2)

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

…

𝑥

𝑁

]

=

[

𝑏

1

𝑏

2

(1)

𝑏

3

(2)

…

𝑏

𝑁

(2)

]

(20b)

gdzie

𝑎

𝑖𝑘

(2)

= 𝑎

𝑖𝑘

(1)

− 𝑙

𝑖2

⋅ 𝑎

2𝑘

(1)

,

𝑏

𝑖

(2)

= 𝑏

𝑖

(1)

− 𝑙

𝑖2

⋅ 𝑏

2

(1)

, przy 𝑙

𝑖2

=

𝑎

𝑖2

(1)

𝑎

22

(1)

oraz 𝑖, 𝑘 = 3, 4, … , 𝑁

1

R. Salamon: MATLAB. Podstawy i zastosowania. PGWETI KSEM Gdańsk 2008

8

 Powyższe postępowanie stosujemy do otrzymanego w drugim kroku równania (20b) i kolejno w

następnych krokach do następnych równań ze zwiększającą się ilością zer pod przekątną główną
przetwarzanej macierzy; po (N-1) krokach otrzymujemy

[

𝑎

11

𝑎

12

𝑎

13

… 𝑎

1𝑁

0

𝑎

22

(1)

𝑎

23

(1)

… 𝑎

2𝑁

(1)

0

0

𝑎

33

(2)

… 𝑎

3𝑁

(2)

…

…

…

… …

0

0

0

… 𝑎

𝑁𝑁

(𝑁−1)

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

…

𝑥

𝑁

]

=

[

𝑏

1

𝑏

2

(1)

𝑏

3

(2)

…

𝑏

𝑛

(𝑁−1)

]

(21)

gdzie

𝑎

𝑖𝑘

(𝑠)

= 𝑎

𝑖𝑘

(𝑠-1)

− 𝑙

𝑖𝑠

⋅ 𝑎

𝑠𝑘

(𝑠-1)

, 𝑏

𝑖

(𝑠)

= 𝑏

𝑖

(𝑠-1)

− 𝑙

𝑖𝑠

⋅ 𝑏

𝑠

(𝑠-1)

,

przy 𝑙

𝑖𝑠

=

𝑎

𝑖𝑠

(𝑠-1)

𝑎

𝑠𝑠

(𝑠-1)

, 𝑠 = 1, 2, … , 𝑁-1, 𝑖 = 𝑠+1, 𝑠+2, … , 𝑁, 𝑘 = 𝑖-1, 𝑖, … 𝑁 (przy czym 𝑎

𝑝𝑞

(0)

= 𝑎

𝑝𝑞

).

W przypadku macierzy z przykładu liczbowego (15) po kroku pierwszym, drugim i trzecim otrzymujemy

kolejno:

[

2 −3

4

5

0

10 −16

−8

0 −

3
2

−10 −

11

2

0

1

−1 −10]

⋅

[

𝑥

1

𝑥

2

𝑥

3

𝑥

4

]

=

[

28

−60

−49

−41]

,

[

2 −3

4

5

0

10 −16

−8

0

0 −

38

5

−

43
10

0

0

3
5

−

46

5

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

𝑥

4

]

=

[

28

−60

−40

−35]

,

[

2 −3

4

5

0 10 −16

−8

0

0 −

38

5

−

43
10

0

0

0 −

725

76

]

⋅

[

𝑥

1

𝑥

2

𝑥

3

𝑥

4

]

=

[

28

−60

−40

−

725

19

]

.

(22)

Równanie skalarne reprezentowane przez ostatni wiersz ostatniego z równań macierzowych (22) jest

równaniem z jedną niewiadomą (konkretnie równaniem

−

725

76

𝑥

4

=

−

725

19

), co pozwala na łatwe obliczenie wiel-

kości 𝑥

4

i następnie wdrożenie procedury wstecznego obliczania kolejno 𝑥

3

, 𝑥

2

, 𝑥

1

. Jednak zajmiemy się dalej

pewnym wariantem metody Gaussa noszącym nazwę metody LU. Polega ona na przyjęciu lewostronnej ma-
cierzy ze wzoru (21) jako macierzy 𝑼 (trójkątna górna) i utworzeniu dodatkowej macierzy 𝑳 (trójkątnej dolnej)
zbudowanej ze współczynników 𝑙

𝑖𝑠

. Dla przypadku macierzy o wymiarach 4x4 (𝑁 = 4) otrzymujemy

𝑼 =

[

𝑎

11

𝑎

12

𝑎

13

𝑎

14

0

𝑎

22

(1)

𝑎

23

(1)

𝑎

24

(1)

0

0

𝑎

33

(2)

𝑎

34

(2)

0

0

0

𝑎

44

(3)

]

, 𝑳 =

[

1

0

0

0

𝑙

21

1

0

0

𝑙

31

𝑙

32

1

0

𝑙

41

𝑙

42

𝑙

43

1]

(23)

Macierze te oprócz intrygującej formy graficznej mają następujące pożyteczne właściwości

𝑳 ∙ 𝑼 = 𝑨, |𝑳| = 1, |𝑨| = |𝑼|

(24)

Wzory (23) i (24) wskazują na korzystniejszy od tradycyjnego sposób obliczania wyznacznika macierzy 𝑨;

dodajmy, sposób wykorzystywany w operacji det(A) MATLAB-a.

Poszukiwanie wektora rozwiązań 𝒙 możemy teraz przedstawić w postaci

𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃 ,

(25)

𝒚 = 𝑳

−𝟏

∙ 𝒃 ,

𝒙 = 𝑼

−𝟏

∙ 𝒚 ,

gdzie 𝒚 to pomocniczy wektor kolumnowy w procedurze (25).

Dla rozważanego przykładu liczbowego (15) mamy

𝑨

=

[

2 −3

4

5

4

4 −8

2

3 −3 −4

2

2 −2

3 −5

]

,

"

𝒃

= [

28

−4

−7

−13

] ,

"

𝒙

= [

𝑥

1

𝑥

2

𝑥

3

𝑥

4

],

𝑨 ⋅ 𝒙 = 𝒃

,

𝑨 = 𝑳 ⋅ 𝑼,

(26)

9

"

𝑼 =

[

2 −3

4

5

0

10 −16

−8

0

0

−38

5

−43

10

0

0

0

−725

76

]

i 𝑼

−𝟏

=

[

1

2

3

20

−1

19

4

25

0

1

10

−4

19

8

725

0

0

−5

38

43

725

0

0

0

−76

725

]

, 𝑳 =

[

1

0

0

0

2

1

0

0

3

2

3

20

1

0

1

1

10

−

3

38

1

]

i 𝑳

−𝟏

=

[

1

0

0

0

−2

1

0

0

−6

5

−3

20

1

0

−17

19

−17

152

3

38

1

]

𝒚 = 𝑳

−𝟏

∙ 𝒃 =

[

28

−60
−40

−725

19

]

,

𝒙 = 𝑼

−𝟏

∙ 𝒚 = [

1
2
3
4

]

Wykonane w MATLAB-ie przeliczenia zgodne z (25), ale o ograniczonej dokładności, dają w rezultacie

𝒚 = [

28

−60
−40

−38.1560

]

,

𝒙 = [

1.0003
2.0000
3.0001
3.9998

]

.

(27)

Widoczna jest niewielka różnica w porównaniu do rozwiązania (18) i (26) otrzymanego na drodze anali-

tycznej bez stosowania ograniczonej długości reprezentacji dziesiętnej ułamków.

Eliminację Gaussa prowadzącą do otrzymania macierzy LU można także zastosować nie do oryginalnego

równania (15), tylko do układu równań z przestawionymi wierszami. Najbardziej oczywista potrzeba przesta-
wienia wierszy może w ramach algorytmu eliminacji wynikać z potrzeby ominięcia sytuacji dzielenia przez
zero, ale możemy także w ten sposób zwiększać dokładność obliczeń.

Wspomnianą optymalizację dokładności obliczeń w odniesieniu do układu równań (19) otrzymujemy

w MATLAB-ie po zastosowaniu do dekompozycji funkcji lu MATLAB-a z argumentami wyjściowymi 𝑳, 𝑼, 𝑷:

[L,U,P] = lu(A).

(28)

Funkcja ta w miejsce macierzy 𝑨 zwraca nam trzy macierze, które możemy podstawić do równania (19). Po
podstawieniu, zapis tego równania przyjmuje postać

𝑷 ∙ 𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃 .

(29)

Zawierająca informację o przestawieniach wierszy macierz 𝑷 ma następujące właściwości

𝑷 = 𝑷

−𝟏

,

|𝑷| = 1 .

(30)

Stąd rozwiązanie równania (19) możemy w tym przypadku uzyskać na podstawie wyrażenia

𝑳 ∙ 𝑼 ∙ 𝒙 = 𝒃

𝒑

(31)

gdzie

𝒃

𝒑

= 𝑷 ∙ 𝒃.

W celu ilustracji rozważań powróćmy jeszcze raz do przykładu liczbowego (15).
Po dekompozycji [L, U, P], wykonanej matlabowym poleceniem lu(A), otrzymujemy

𝑷 = [

0 1 0 0
0 0 1 0
1 0 0 0
0 0 0 1

] , 𝑳 = [

1

0

0

0

0.75

1

0

0

0.5

0.83333

1

0

0.5

0.66667

0.89474

1

], 𝑼 = [

4

4

−8

2

0

−6

2

0.5

0

0

6.33333

3.58333

0.5

0.6667

0.89474

−9.53948

].

(32a)

Wtedy według równań (25) po uwzględnieniu w nich zmian (31) otrzymujemy:

10

𝒃

𝒑

=

[

−4
−7

28

−13

]

,

𝒚

𝒑

= 𝑳

−𝟏

∙ 𝒃

𝒑

= [

−4
−4

33.33333

−38.15790

]

,

𝒙 = 𝑼

−𝟏

∙ 𝒚

𝒑

= [

1
2
3
4

]

(32b)

Zauważmy, że po porównaniu macierzy 𝒃

𝒑

z macierzą 𝒃 dysponujemy informacją, że w równaniu (15)

w trakcie działania procedury wiersz numer 1 został przesunięty w miejsce wiersza numer 3, co oczywiście
wynika także bezpośrednio z postaci macierzy permutacji 𝑷 (ściśle rzecz biorąc, macierz ta mówi nam, w tym
konkretnym przypadku, że nowym pierwszym wierszem ma być stary wiersz drugi, nowym drugim wierszem
ma być stary wiersz trzeci, nowym trzecim wierszem ma być stary wiersz pierwszy, wreszcie nowym czwartym
wierszem ma pozostać stary wiersz czwarty).

Realizując poleceniem

lu(A)

opisaną powyżej procedurę stosujemy optymalizowaną dekompozycję

[L, U, P] w sposób jawny, ale oprócz tego dekompozycja taka jest automatycznie wykonywana niejawnie
w ramach następujących poleceń MATLAB-a:

 w przypadku równania postaci (16), gdy je rozwiązujemy poleceniem

x=A\b

(33)

 w przypadku obliczania wyznacznika macierzy A

DA = det(A)

(34)

4. Pliki MATLAB-a do wspomagania analizy obwodów

function [x]=cramer(A,b)
%Metoda Cramera rozwiązywania układów równań liniowych
%Zapis układu równań A*x=y
N=rank(A);
for k=1:N
K=A(:,k); % Kolumna o numerze k w macierzy A na starcie
A(:,k)=b; % Macierze A

1

, A

2

,…,A

m

DK(k)=det(A); % Wyznaczniki D

1

, D

2

,…,D

m

A(:,k)=K; %Powrót do startowej postaci macierzy A
end
D=det(A); % Wyznacznik główny
x=DK./D; % Tablicowe dzielenie wyznaczników
%Koniec cramer.m

function [LA,UA]=dekompozycja(A)

% Dekompozycja (rozkład) macierzy A na dwie macierze
%„trójkątną-dolną” LA i „trójkątną- górną” UA

N=rank(A); % Określa liczbę niezależnych równań
el=eye(N); % Generuje macierz z "jedynkami" na głównej przekątnej
a=zeros(N); % Generuje macierz wypełnioną zerami

a=A; % Przepisuje elementy macierzy

for s=1:(N-1) % Litera s oznacza numer kroku w metodzie Gaussa
% Litery i oraz k oznaczają numery wiersza i kolumny
n=s+1;
for i=n:N

el(i,s)=a(i,s)/a(s,s)

for k=s:N

a(i,k)=a(i,k)-el(i,s)*a(s,k);

end

end
end

LA=el; % Macierz trójkątna dolna
UA=a; % Macierz trójkątna górna
% Koniec dekompozycja.m

Cw2_1.m

%***Program wspomagający poszukiwanie prądów oczkowych***
%
R=[1,2,3,4,5,6,7,8];E1=10;J3=5; g=5; % dane o obwodzie (opory,źródła)
%
Roo= [R(1)+R(2)+R(5) , -R(5) , -R(2) , 0 ;

11

-R(5) , R(5)+R(7)+R(8) , 0 , -R(7) ;
-R(2) , 0 , R(2)+R(3)+R(4), -R(4) ;
-g*R(6)*R(5) , -R(7)+g*R(6)*R(5), -R(4) ,R(4)+R(6)+R(7)];
%
% Roo to macierz współczynników (macierz oporności oczkowych)
%
E=[E1;0;-J3*R(3);0];

%wektor pobudzeń (wektor prądów źródłowych)

%
Io_cramer=cramer(Roo,E)

%metoda wyznaczników Cramera

%
[LRoo,URoo]=dekompozycja(Roo) %nasz wariant LU metody Gaussa
y=(LRoo^(-1))*E;

%wektor pomocniczy y

Io_gauss_lu=(URoo^(-1))*y

%wektor prądów oczkowych

%
Io_lup=Roo\E

%matlabowy [L,U,P] wariant metody Gaussa

5. Symulacja obwodu elektrycznego za pomocą programu SPICE
W narzędziach symulacyjnych rodziny SPICE poszukiwanie rozwiązania układu równań liniowych jest

oparte na tej samej metodzie numerycznej LU, która jest opisana w rozdziale 3. Jednak w odróżnieniu od
właściwości narzędzia MATLAB, nie ma w programie SPICE dla WINDOWS potrzeby zewnętrznego formuło-
wania równań opisujących obwód. Na przykład w edytorze narzędzia Micro-Cap (jest to „odmiana” SPICE’a)
tworzymy wirtualny obwód złożony: ze źródeł niezależnych, źródeł sterowanych, rezystorów, cewek induk-
cyjnych i kondensatorów i w miarę potrzeby z innych elementów pochodzących z banku elementów tego
narzędzia. Wczesne wersje narzędzia SPICE i jego odmian do wprowadzania danych o obwodzie, w tym o
jego strukturze, zamiast edytora graficznego wykorzystywały tzw. tablicę połączeń. W ten sposób SPICE i po-
krewne jemu narzędzia – w zewnętrznym oglądzie – raczej symulują właściwości obwodu niż wspomagają
poszukiwanie rozwiązań równań opisujących taki obwód. Pomijamy sposób, w jaki to czynią, zauważymy
tylko, że nie wykorzystują równań oczkowych, a stosują równania węzłowe uogólnione (o czym w trakcie
następnego ćwiczenia).

Z wielu wariantów SPICE’a w laboratorium jest do dyspozycji narzędzie Micro-Cap w wersjach demon-

stracyjnych 10 i 11. Obwód prądu stałego z rysunku 2a, wprowadzony do programu Micro-Cap wygląda jak na
rysunku 4.

Rysunek 2a i rysunek 4 różnią się przede wszystkim dodatkowymi fragmentami, które pojawiły się z tego

względu, że w banku elementów narzędzia Micro-Cap mamy rysunki źródeł sterowanych w postaci czwórni-
kowej. Edytor nie pozwala na wygodny sposób określania kierunków odniesienia dla przepływu prądów w ga-
łęziach oraz nadawania numerów porządkowych węzłów, jednak z tym poradziliśmy sobie, jak widać. We
wprowadzonym do edytora programu obwodzie można poprzez kliknięcie myszką w obrębie wybranego ele-
mentu zmienić w otwierającym się wtedy okienku wartość liczbową parametru (E1, J3, g, R1, R2,…, R8) (na
rysunku 4 mamy przykładowo R1=1 [Ω], g6= 5 [S]

1

).

Rys. 4. Schemat obwodu z rys. 2a utworzony w programie Micro-Cap

1

g6=5 [S] to parametr zwany transkonduktancją i oznaczany wcześniej przez g dla ZPSN, oznaczonego na rysunku 4 przez I9ofU5,

rozpiętego między zaciskami (0, 3) – wejście i (4, 5) – wyjście (uwaga: węzeł 0, to węzeł masy).

12

Jeżeli z kilku wariantów badania obwodu DC (prądu stałego) wybierzemy Dynamic DC otrzymamy na-

pięcia węzłowe (V1, V2, V3, V4), czyli między innymi potencjały każdego z ponumerowanych węzłów wzglę-
dem napięcia odniesienia (masa, ziemia) oraz prądy gałęziowe (I1, I2, …, I9). Jeżeli dodatkowo zrealizujemy
w polach tekstowych pewne działania arytmetyczne, możemy uzyskać zestaw wyników liczbowych, widoczny
na rysunku 4 w prawej jego części, potrzebny do sprawdzenia poprawności wyników naszych obliczeń wspo-
maganych narzędziem MATLAB.

Pytania kontrolne

1) Na czym polega tzw. metoda LU rozwiązywania układu równań liniowych?
2) Podać szczegóły tworzenia tzw. źródła Thevenina dla wieloelementowego obwodu liniowego.
3) Podać szczegóły tworzenia tzw. źródła Nortona dla wieloelementowego obwodu liniowego.
4) Podać szczegóły metody rozwiązywania układu równań liniowych związanej z nazwiskiem Cramera.
5) W jaki sposób wprowadzamy dane badanego obwodu do programu Micro-Cap?
6) Na czym polega tzw. dekompozycja LU macierzy współczynników obwodu?
7) Przez element zwierający zaciski dwójnika liniowego płynął prąd 2 A, a gdy ten element zwierający

zastąpiono opornikiem 1 Ω, to przez ten opornik płynął prąd 1 A. Podaj parametry zastępczego źródła
Thevenina dla tego dwójnika. Odpowiedź uzasadnij.

8) Co to jest strzałkowanie skojarzone?
9) Objaśnij jak jest rozumiane MATLAB-owe „dzielenie” A\b oraz objaśnij czym są (macierze, wektory,

wymiary) poszczególne elementy (A, b, wynik) tego działania.

10) Rozwiąż metodą eliminacji Gaussa równanie [

1 1
4 2

] ⋅ [

𝑥

1

𝑥

2

] = [3

8

]. Zapisz kolejne kroki obliczeń.