Ministerstwo Edukacji Republiki Białoruś

instytucja edukacyjna

· kontur3(X, Tak, Z) – buduje linie konturowe dla powierzchni uzyskanej przez nałożenie na siebie trójwymiarowej figury kilkoma siecznymi płaszczyznami usytuowanymi równolegle do płaszczyzny odniesienia figury.

Poniższy przykład pokazuje, jak opisane polecenia mogą być użyte do wykreślenia powierzchni Z = sin(X) cos(X).

>> =siatka(-3:0.1:3;-3:0.1:3); Z=sin(X).*cos(X);

>> subplot(3,2,1), plot3(X,Y,Z) % Rysunek 4.3 (a)

>> subplot(3,2,2), mesh(X,Y,Z) % Rysunek 4.3 (b)

>> subplot(3,2,3), surf(X, Y, Z) % Rysunek 4.3 (c)

>> subplot(3,2,4), surfc(X, Y, Z) % Rysunek 4.3(d)

>> subplot(3,2,5),meshz(X,Y,Z) % Rysunek 4.3(e)

>> subplot(3,2,6),contour3(X, Y, Z) % Rysunek 4.3(e)

4.3. Formatowanie wykresu

System Matlab zapewnia możliwość dostosowywania i korygowania właściwości wykresów zarówno za pomocą interfejsu okna graficznego, jak i poprzez ustawienie odpowiednich poleceń i parametrów graficznych. Tabela 4.3.1. podano kilka prostych technik formatowania wykresów.

Tabela 4.2 - Formatowanie wykresu

Akcja

Narzędzia GUI

Przełącz na tryb edycji.

Kliknij przycisk EdytowaćIntrygować na pasku narzędzi okna graficznego.

Formatowanie linii i znaczników punkty odniesienia wykresy.

W trybie edycji kliknij dwukrotnie linię wykresu lewym przyciskiem myszy. W wyświetlonym oknie własnośćRedaktor-linia ustaw wszystkie niezbędne parametry linii (grubość, styl, kolor itp.).

działka(X,Y,S), działka3(X,Y,Z,S)

(Opis poleceń podano w paragrafie 4.1 i paragrafie 4.2)

Formatowanie osi wykresu.

W trybie edycji kliknij dwukrotnie oś wykresu. W wyświetlonym oknie własnośćRedaktor-osie ustaw wszystkie niezbędne parametry osi.

Tytuł wykresu i etykiety osi mogą być również używane z poleceniami WstawićTytuł, WstawićXetykieta,WstawićEtykieta Y menu główne okna graficznego.

osie– kontroluje właściwości osi.

krata- włącza i wyłącza siatkę współrzędnych.

xetykieta(S),etykieta y (S),zetykieta(S)– ustawia etykiety w pobliżu osi. Tutaj S jest ciągiem stałym lub zmienną.

tytuł(S)- wyświetla tytuł wykresu

Rysowanie napisów bezpośrednio na wykresie.

Kliknij przycisk WstawićTekst, ustal miejsce napisu jednym kliknięciem myszy i wprowadź żądany tekst.

tekst(x,Tak,S)– dodaje do wykresu dwuwymiarowego tekst określony ciągiem S, tak aby początek tekstu znajdował się w punkcie o współrzędnych (X, Y).

tekst(x,Tak,Z,S)- dodaje tekst do wykresu 3D.

Rysowanie linii i linii strzałkami bezpośrednio na wykresie

Kliknij jeden z przycisków WstawićStrzałka lub Wstawićlinia. Ustaw wskaźnik myszy w żądanym miejscu na wykresie i trzymając wciśnięty lewy przycisk myszy, narysuj linię.

Budowanie legendy

wstawić, a potem polecenie legenda.

legenda(S1,S2,S3,…)– dodaje do aktualnego wykresu legendę z objaśnieniami w postaci ciągów określonych na liście parametrów.

Skala kolorów wyjściowych

W menu głównym okna graficznego wybierz polecenie wstawić, a potem polecenie pasek kolorów.

pasek kolorów('zielony'),pasek kolorów('poziomo”)– Wyświetla pionową lub poziomą skalę kolorów.

Obrót wykresu

Kliknij przycisk Obróć 3D i obracaj wykres za pomocą myszy (może być również używany do wykresów dwuwymiarowych).

obróć3d– ustawia obrót figury trójwymiarowej.

z warunkami brzegowymi tak(t 0 , mieć skłonności, p) = tak, gdzie mieć skłonności, t 0 punktów początkowych i końcowych interwałów. Parametr t(zmienna niezależna) niekoniecznie oznacza czas, chociaż najczęściej rozwiązania DE poszukuje się w dziedzinie czasu. System DE w postaci Cauchy'ego jest napisany podobnie do (1.1), ale poniżej tak w tym przypadku zakładany jest wektor kolumnowy zmiennych zależnych. Wektor p ustawia warunki początkowe.

Aby rozwiązać DE drugiego i wyższy porządek muszą zostać zredukowane do systemu DE pierwszego rzędu.

Możliwe są funkcje sterujące, które nie są dozwolone w odniesieniu do pochodnej:

F(t, tak, dy/dt) = 0.

(1.2)

Równań (1.2) zazwyczaj nie da się sprowadzić analitycznie do postaci (1.1). Jednak numeryczne rozwiązanie szczególnych trudności nie powoduje wystarczającego określenia f(tak, t) rozwiąż (1.2) numerycznie w odniesieniu do pochodnej dla danego tak oraz t.

solwery ODE

Aby rozwiązać systemy ODE w MATLAB, zaimplementowano różne metody numeryczne. Ich implementacje noszą nazwę rozwiązywacze ODA.

W tej sekcji nazwa generyczna solver oznacza jedną z możliwych metod numerycznych rozwiązywania ODE: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c lub pdepe.

Solvery implementują następujące metody rozwiązywania systemów DE:

Ode45 jednoetapowe jawne metody Runge-Kutta czwartego i piątego rzędu zmodyfikowane przez Dormanda i Prinza. Jest to klasyczna metoda zalecana do wstępnego wypróbowania rozwiązania. W wielu przypadkach daje to dobre wyniki, jeśli układ równań do rozwiązania nie jest sztywny.

Ode23 jednoetapowe jawne metody Rungego-Kutty II i IV rzędu w modyfikacji Bogackiego i Champina. Przy umiarkowanej sztywności systemu ODE i niskich wymaganiach dotyczących dokładności, ta metoda może zapewnić wzrost szybkości rozwiązania.

Ode113 wieloetapowa metoda Adamsa-Bashwortha-Multona klasy predyktora-korektora zmiennego rzędu. Jest to metoda adaptacyjna, która może zapewnić wysoką dokładność rozwiązania.

Ode15s wielostopniowa metoda kolejności zmiennych (od 1 do 5, domyślnie 5) z wykorzystaniem formuł numerycznych „wstecznego różnicowania”. Jest to metoda adaptacyjna i powinna być stosowana, jeśli solver ode45 nie zapewnia rozwiązania, a system sterowania jest sztywny.

Ode23s to jednoetapowa metoda wykorzystująca zmodyfikowaną formułę Rosenbrocka drugiego rzędu. Może zapewnić dużą szybkość obliczeniową przy niskiej dokładności rozwiązywania sztywnego układu DE.

Ode23t ukryta metoda trapezowa z interpolacją. Metoda ta daje dobre rezultaty przy rozwiązywaniu problemów opisujących układy oscylacyjne z prawie harmonicznym sygnałem wyjściowym. W przypadku umiarkowanie twardych systemów zdalne sterowanie może zapewnić wysoką dokładność rozwiązania.

Niejawna metoda Ode23tb Runge Kutta na początku rozwiązania, a następnie metoda wykorzystująca formuły „wstecznego różnicowania” drugiego rzędu. Mimo stosunkowo małej dokładności metoda ta może być bardziej wydajna niż ode15s.

Bvp4c służy do problemu wartości brzegowych układów układów sterowania postaci tak′ = f(t, tak), F(tak(a), tak(b), p) = 0 (długa forma układ równań Cauchy'ego). Rozwiązywane przez niego problemy nazywane są dwupunktowymi problemami z wartościami brzegowymi, ponieważ rozwiązanie poszukuje się poprzez ustalenie warunków brzegowych zarówno na początku, jak i na końcu przedziału rozwiązania.

Wszystkie solwery mogą rozwiązywać układy jawnych równań tak′ = F(t, tak), a do rozwiązywania sztywnych układów równań zaleca się używanie tylko specjalnych solverów ode15s, ode23s, ode23t, ode23tb.

Korzystanie z solwerów ODE

tspan wektor definiujący przedział całkowania [ t 0 t ostateczna]. Aby uzyskać rozwiązania w określonych punktach w czasie t 0 , t 1 , …, t ostateczna(uporządkowane w kolejności malejącej lub rosnącej) należy stosować tspan = [t 0 t 1 … t ostateczna];

tak 0 wektor warunków początkowych;

Argument opcji generowany przez funkcję odeset (inna funkcja odeget lub bvpget (tylko dla bvp4c) pozwala wyświetlić opcje ustawione domyślnie lub za pomocą funkcji odeset/bvpset);

p 1, p 2,… dowolne parametry przekazywane do funkcji F;

T, Tak macierz decyzyjna Tak, gdzie każdy wiersz odpowiada czasowi zwróconemu w wektorze kolumny T.

Przejdźmy do opisu składni funkcji rozwiązywania układów zdalnego sterowania (nazwa solver oznacza dowolną z przedstawionych powyżej funkcji).

[T,Tak]=rozwiązujący(@ F,tspan,tak 0) integruje system zdalnego sterowania formularza tak′ = F(t, tak) na interwale tspan z warunkami początkowymi tak 0 . @F deskryptor funkcji ODE (można również określić funkcję w postaci " F"). Każdy wiersz w tablicy rozwiązań Tak odpowiada wartości czasu zwróconej w wektorze kolumny T.

[T,Tak]=rozwiązujący(@ F,tspan,tak 0 ,opcje) daje rozwiązanie podobne do powyższego, ale z opcjami określonymi przez wartości argumentu options utworzonego przez funkcję odeset. Powszechnie używane parametry obejmują dopuszczalny błąd względny RelTol (domyślnie 1e3) oraz wektor dozwolone wartości błąd absolutny AbsTol (domyślnie wszystkie komponenty to 1e6).

[T,Tak]=rozwiązujący(@ F,tspan,tak 0 ,opcje, p 1 ,p 2...) daje rozwiązanie podobne do opisanego powyżej, przechodząc Dodatkowe opcje p 1 , p 2 , … w m-plik F ilekroć jest nazywany. Użyj options=, jeśli nie podano żadnych opcji.

Rozwiązanie ODE pierwszego rzędu

PROCEDURA PRACY

· Strona tytułowa;

dane początkowe wariantu;

· rozwiązanie problemu;

wyniki rozwiązania problemu.

Przykład

Znajdź rozwiązanie równania różniczkowego na przedziale, dla którego w(1,7) = 5,3.

Utwórz funkcję użytkownika w oknie poleceń

[e-mail chroniony](x,y);

W składni funkcji @(x,y) x zmienna niezależna, tak zmienna zależna, x-sałata( tak/Liczba Pi) po prawej stronie DE.

Proces rozwiązywania odbywa się poprzez wywołanie solvera w oknie poleceń z następującą instrukcją:

Ode23(g,,);

Konstrukcję grafu z siatką realizują następujące operatory:

Wynik pokazano na ryc. 1,1

Ryż. 1.2.1. Wizualizacja rozwiązania numerycznego

ĆWICZENIE

1. Znajdź rozwiązania DE pierwszego rzędu , spełniające warunki początkowe y(x 0 ) = y 0 w przedziale [ a,b].

2. Konstruować wykresy funkcji.

Opcje zadań.

numer opcji		y(x 0 )=y 0	[a,b]
		tak 0 (1,8)=2,6
		tak 0 (0,6)=0,8
		tak 0 (2,1)=2,5
		tak 0 (0,5)=0,6
		tak 0 (1,4)=2,2
		tak 0 (1,7)=5,3
		tak 0 (1,4)=2,5
		tak 0 (1,6)=4,6
		tak 0 (1,8)=2,6
		tak 0 (1,7)=5,3
		tak 0 (0,4)=0,8
		tak 0 (1,2)=1,4

Laboratorium #2

Rozwiązanie systemów ODE

CEL PRACY

Aby sformułować pomysły uczniów na temat korzystania z systemów zdalnego sterowania w różne pola; zaszczepić umiejętność rozwiązania problemu Cauchy'ego dla systemów zdalnego sterowania.

PROCEDURA PRACY

1. Przestudiuj część teoretyczną. Wykonaj zadania odpowiadające numerowi twojej opcji i pokaż je nauczycielowi.

2. Przygotuj raport z laboratorium, który powinien zawierać:

· Strona tytułowa;

dane początkowe wariantu;

· rozwiązanie problemu;

wyniki rozwiązania problemu.

Przykład

Rozwiąż system

za pomocą solvera ode23().

Rozwiązanie:

1. Utwórz w edytorze plik m funkcji do obliczania właściwych części DE.

Niech nazwa w edytorze plików to sisdu.m, wtedy funkcja może wyglądać tak:

funkcja z=sisdu(t,y)

z1=-3*y(2)+cos(t)-exp(t);

z2=4*y(2)-cos(t)+2*exp(t);

>> t0=0;tf=5;y0=[-3/17,4/17];

>> =ode23("sisdu",,y0);

>>wykres(t,y)

Ryż. 1.3.1. Wizualizacja rozwiązania numerycznego otrzymanego za pomocą funkcji ode23.

1. Co to znaczy rozwiązać problem Cauchy'ego dla układu sterowania?

2. Jakie istnieją metody rozwiązywania systemów zdalnego sterowania?

ĆWICZENIE

1. Znajdź rozwiązanie do zdalnego sterowania

spełnienie warunków początkowych na przedziale ;

2. Buduj wykresy funkcji.

Na przykład podana jest funkcja decyzyjna 8. opcji:

funkcja z=ssisdu(t,y)

% opcja 8

z1=-a*y(1)+a*y(2);

z2=a*r(1)-(a-m)*r(2)+2*m*r(3);

z3=a*r(2)-(a-m)*r(3)+3*m*r(4);

z4=a*r(3)-3*m*r(4);

>> =ode23("ssisdu",,);

>>wykres(t,100*r)

Ryż. 1.3.2. Wizualizacja rozwiązania numerycznego otrzymanego za pomocą funkcji ode23.

Opcje zadań.

numer opcji	Zadania
a	m
	0,1	1,2
	0,2	1,5
	0,3	1,7
	0,4	1,9
	0,5
	0,6	1,9
	0,7	2,3
	0,8	2,7
	0,9
	0,1	1,5
	0,2	1,1
	0,3

Laboratorium #3

1.4 Rozwiązanie ODE n-tego rzędu

CEL PRACY

Formułowanie pomysłów uczniów na temat wykorzystania systemów sterowania wyższego rzędu w różnych dziedzinach; zaszczepić umiejętność rozwiązania problemu Cauchy'ego dla DE wyższego rzędu przy użyciu programy użytkowe; rozwijać umiejętności sprawdzania uzyskanych wyników.

PROCEDURA PRACY

1. Przestudiuj część teoretyczną. Wykonaj zadania odpowiadające numerowi twojej opcji i pokaż je nauczycielowi.

2. Przygotuj raport z laboratorium, który powinien zawierać:

· Strona tytułowa;

dane początkowe wariantu;

· rozwiązanie problemu;

wyniki rozwiązania problemu.

Przykład 1

Rozwiąż drugie zamówienie DE w danych warunkach początkowych .

Rozwiązanie:

Najpierw wprowadzamy DE do systemu:

1. Utwórz plik m funkcji do obliczania właściwych części DE.

Niech nazwa pliku to sisdu_3.m, wtedy funkcja może wyglądać tak:

funkcja z=sisdu_3(x,y)

z2=6*x*exp(x)+2*y(2)+y(1);

2. Wykonaj następujące czynności:

>> x0=0;xf=10;y0=;

>> =od23("sisdu_3",,y0);

>>wykres(x,y(:,1))

Ryż. 1.4.1. Wizualizacja rozwiązania numerycznego otrzymanego za pomocą funkcji ode23.

PRZYKŁADOWE PYTANIA DOTYCZĄCE OCHRONY PRACY

1. Co to znaczy rozwiązać problem Cauchy'ego dla DE wyższego rzędu?

2. Jak zabrać DU m zamówienie do systemu zdalnego sterowania?

ĆWICZENIE

1. Znajdź rozwiązanie DE spełniające warunki początkowe na przedziale .

2. Buduj wykresy funkcji.

Opcje zadań.

numer opcji	Zadania
Równania	Warunki początkowe

Laboratorium #4 – 5

Systemy dynamiczne (DS)

CEL PRACY

Zapoznanie studentów z podstawowymi pojęciami DS, ich klasyfikacją, przestrzenią fazową DS, kinematyczną interpretacją DS, ewolucją DS. Równanie ruchu wahadła. Dynamika oscylatora Van der Pol.

2. System dynamiczny (DS) obiekt matematyczny odpowiadający rzeczywistym systemom (fizycznym, chemicznym, biologicznym itp.), którego ewolucja jest jednoznacznie zdeterminowana przez stan początkowy. DS jest określany przez układ równań (różnicowy, różnicowy, całkowy itp.), który pozwala na istnienie unikalnego rozwiązania dla każdego warunku początkowego w nieskończonym przedziale czasu.

Stan DS opisuje zbiór zmiennych dobranych ze względu na naturalność ich interpretacji, prostotę opisu, symetrię itp. Zbiór stanów DS tworzy przestrzeń fazową, każdy stan odpowiada punktowi w niej, a ewolucja jest reprezentowana przez trajektorie (fazy). Aby określić bliskość stanów, w przestrzeni fazowej DS wprowadza się pojęcie odległości. Zbiór stanów w ustalonym punkcie w czasie charakteryzuje się objętością fazową.

Opis DS w sensie ustalenia prawa ewolucji pozwala również na dużą różnorodność: jest realizowany za pomocą równań różniczkowych, odwzorowań dyskretnych, za pomocą teorii grafów, teorii łańcuchów Markowa itp. Wybór jednej z metod opisu wyznacza określony typ modelu matematycznego odpowiedniego DS.

Model matematyczny DS uważa się za dane, jeśli wprowadzane są zmienne dynamiczne (współrzędne) układu, które jednoznacznie określają jego stan, oraz określone jest prawo ewolucji stanu w czasie.

W zależności od stopnia aproksymacji do tego samego systemu można przypisać różne modele matematyczne. Badanie rzeczywistych systemów przebiega ścieżką badania odpowiednich modeli matematycznych, których ulepszenie i rozwój determinuje analiza wyników eksperymentalnych i teoretycznych podczas ich porównywania. W związku z tym przez system dynamiczny będziemy rozumieć właśnie jego model matematyczny. Badając ten sam DS (np. ruch wahadła), w zależności od stopnia uwzględnienia różnych czynników, otrzymamy różne modele matematyczne.

Interpretacyjny język programowania systemu MATLAB jest zaprojektowany w taki sposób, aby dowolne (czasem bardzo złożone) obliczenia można było wykonywać w trybie obliczeń bezpośrednich, czyli bez przygotowania programu przez użytkownika. Jednocześnie MATLAB pełni funkcje superkalkulatora i pracuje w trybie wiersz poleceń.

Praca z systemem ma charakter interaktywny i odbywa się zgodnie z zasadą „zadałem pytanie – otrzymałem odpowiedź”. Użytkownik wpisuje wyliczone wyrażenie na klawiaturze, edytuje je (jeśli to konieczne) w wierszu poleceń i uzupełnia dane wejściowe, naciskając klawisz ENTER. Jako przykład na rysunku pokazano najprostsze i dość oczywiste obliczenia.

Nawet z tak prostych przykładów można wyciągnąć pouczające wnioski:

* symbol >> służy do wskazania wprowadzania danych początkowych;

* dane wprowadza się za pomocą prostego edytora linii;

* aby zablokować wyjście wyniku obliczeń jakiegoś wyrażenia po nim, musisz ustawić znak; (średnik);

* jeśli nie podano zmiennej dla wartości wyniku obliczenia, to MATLAB przypisuje takiej zmiennej nazwę ans;

* znak przypisania jest znanym matematykom znakiem równości =, a nie znakiem złożonym: =, jak w wielu innych językach programowania i systemach matematycznych;

* funkcje wbudowane (na przykład sin) są pisane małymi literami, a ich argumenty są podane w nawiasach;

* wynik obliczeń wyświetlany jest w wierszach wyjściowych (bez znaku >>);

* dialog odbywa się w stylu „zadałem pytanie – otrzymałem odpowiedź”.

Poniższe przykłady ilustrują, jak MATLAB może być używany do wykonywania prostszych operacji wektorowych. Rysunek pokazuje również okno przeglądarki system plików, który jest dostępny na karcie Bieżący katalog. W trybie poleceń wygodniej jest wywołać okno przeglądarki systemu plików z paska narzędzi, aktywując przycisk za listą katalogów systemowych MATLAB. Mogą wystąpić przypadki odmowy obliczeń, gdy bieżący katalog jest niepoprawnie ustawiony, jeśli m-pliki potrzebne do obliczeń nie zostaną znalezione.

W większości systemów matematycznych obliczenie sin(V) lub exp(V), gdzie V jest wektorem, skutkowałoby błędem, ponieważ sin i exp muszą mieć argument skalarny. Jednak MATLAB jest systemem macierzowym, a wektor jest rodzajem macierzy o rozmiarze 1 × n lub n × 1. Dlatego w naszym przypadku wynikiem obliczenia będzie wektor o tej samej wielkości co argument V, ale elementy zwracanego wektora będą sinusami lub wykładnikami elementów wektora V.

Macierz jest określona jako seria wektorów reprezentujących jej wiersze i jest ujęta w nawiasy kwadratowe. Spacja lub przecinek służy do oddzielenia elementów wektorów, a średnik służy do oddzielenia jednego wektora od drugiego. Do wybrania pojedynczego elementu macierzy M stosuje się wyrażenie w postaci M(j,i), gdzie M to nazwa macierzy, j to numer wiersza, a i to numer kolumny.

Do przeglądania zawartości tablic wygodnie jest używać przeglądarki Workspace. Każdy wektor i macierz w nim zawarte są reprezentowane jako kwadrat z komórkami, po prawej stronie którego wskazany jest rozmiar tablicy. Dwukrotne kliknięcie kwadratu myszą prowadzi do pojawienia się okna Array Editor. Praca z edytorem tablic jest dość oczywista - można nie tylko przeglądać elementy tablicy, ale także je edytować i zastępować.

Jak widać z powyższych przykładów, wprowadzanie wyrażeń początkowych do obliczeń w systemie MATLAB odbywa się w najbardziej powszechnym formacie tekstowym. Wyniki obliczeń wyświetlane są w tym samym formacie, z wyjątkiem graficznych. Oto przykłady rejestrowania obliczeń wykonywanych przez system MATLAB w wierszu poleceń:

Praca z edytorem tablic

Aby rozpocząć, wybierz „Pomoc MATLAB” z menu Pomoc.

>> wpisz grzech

grzech jest wbudowaną funkcją.

>> pomóc grzechowi

SIN(X) jest sinusem elementów X.

Przeciążone metody

>>V=

0.8415 0.9093 0.1411 -0.7568

Błąd przy użyciu ==> ^

Matryca musi być kwadratowa.

Możesz zwrócić uwagę na formę odpowiedzi podczas wykonywania prostych operacji bez określania zmiennej, do której przypisany jest wynik. W takich przypadkach MATLAB sam przypisuje zmienną ans, do której przypisywany jest wynik i której wartość jest następnie wyświetlana.

Formularz wyjściowy i łamanie wiersza w sesji

Należy zwrócić uwagę na cechy wyjścia w systemie MATLAB. Dane wyjściowe rozpoczynają się w nowym wierszu z wciętymi danymi liczbowymi i bez wcięcia danych tekstowych. Aby zaoszczędzić miejsce w tej książce, wyprowadzenie zostanie podane w dalszej części bez tłumaczenia na Nowa linia. Na przykład wyjście wektora wiersza

zostaną podane w formie:

Wyjątkiem jest wyjście wektorów kolumnowych i macierzy - tutaj zostanie zachowana bardziej opisowa i domyślna forma wyjścia MATLAB.

W niektórych przypadkach wprowadzone wyrażenie matematyczne może być tak długie, że jeden wiersz nie wystarczy. Następnie część wyrażenia można przenieść do nowej linii za pomocą znaku wielokropka „...” (3 lub więcej kropek), na przykład:

s = 1 - 1/2 + 1/3 - 1/4 + 1/5 - 1/6 + 1/7 ...

1/8 + 1/9 - 1/10 + 1/11 - 1/12;

Maksymalna liczba znaków w jednej linii trybu poleceń to 4096, aw pliku m nie jest ograniczona, ale praca z tak długimi liniami jest niewygodna. W wczesne wersje w jednym wierszu nie było więcej niż 256 znaków.

Uruchamianie przykładów aplikacji MATLAB z wiersza poleceń

MATLAB ma wiele przykładów aplikacji, z których część można uruchomić bezpośrednio z wiersza poleceń. Na przykład polecenie

uruchamia m-file bench.m demo testu systemu.

Obliczenia i aproksymacja danych w MATLAB

Rozwiązywanie układów liniowych równań algebraicznych i problemów spektralnych

Rozwiązywanie układów liniowych równań algebraicznych jest jednym z głównych problemów obliczeniowych, ponieważ sprowadza się do niego ogromna liczba problemów pojawiających się w różnych stosowanych dziedzinach. Metody numeryczne stosowane do przybliżonego rozwiązywania problemów w mechanice, do obliczania przepływów cieczy i gazów, innych procesów fizycznych w ośrodkach ciągłych, metody obliczania obwodów elektrycznych, metody aproksymacji danych prowadzą do konieczności rozwiązywania układów liniowych równań algebraicznych ( i nie jest to pełna lista źródeł występowania liniowych systemów algebraicznych). Problemy spektralne pojawiają się m.in. w analizie częstotliwościowej konstrukcji, w badaniu stabilności procesów w schematy elektryczne i systemów sterowania oraz w wielu innych zastosowaniach.

Przegląd możliwości MATLAB do rozwiązywania układów liniowych równań algebraicznych

Środowisko MATLAB zostało pierwotnie opracowane do pracy z macierzami (MATLAB jest skrótem od Matrix Laboratory), więc arsenał narzędzi MATLAB do rozwiązywania układów liniowych równań algebraicznych jest dość bogaty i obejmuje:

• rozwiązanie systemów z matrycami kwadratowymi i prostokątnymi;
• rozwiązywanie systemów metodami bezpośrednimi i iteracyjnymi (w tym z możliwością warunkowania wstępnego);
• rozszerzenia macierzy;
• przechowywanie dużych rzadkich macierzy w zwartej formie i specjalnych algorytmów rozwiązywania systemów z takimi macierzami.

Rozwiązywanie systemów za pomocą znaku odwrotnego ukośnika

Najprostszym sposobem rozwiązania systemów jest użycie znaku odwrotnego ukośnika. Załóżmy, że musimy rozwiązać system

Aby to zrobić, wypełnij macierz i wektor kolumn po prawej stronie (prawa strona musi być dokładnie tą kolumną, w przeciwnym razie zostanie wyświetlony błąd o niezgodności wymiarów)

i użyj odwrotnego ukośnika

X = A \ f x = 1 1 1

Zamiast odwrotnego ukośnika możesz wywołać funkcję mldivide

X = miliardvide(A, f)

Wynik będzie taki sam

Obliczając resztę upewniamy się, że rozwiązanie zostało znalezione poprawnie (ogólnie rzecz biorąc, nie zawsze małe wartości składowych rezydualnych wskazują na prawidłowo znalezione rozwiązanie)

F - A*x ans = 0 0 0

Bardzo ważne jest, aby nie mylić A i f, ponieważ podczas wykonywania operacji

nie będzie błędu, ale wyjście

X = 0,2707 0,3439 0,3854

co nie ma nic wspólnego z rozwiązaniem rozważanego systemu. Przeanalizujmy, dlaczego tak się dzieje, spisując układ z „macierzą” f i „prawą stroną” A. W tym przypadku będzie tylko jedna niewiadoma, ponieważ f rozmiar 3 na 1:

Jeśli po prawej stronie znajduje się macierz, to tyle układów zostanie rozwiązanych, ile jest kolumn w macierzy, a każda kolumna jest wektorem prawej strony odpowiedniego układu (tzn. operator odwrotnego ukośnika pozwala rozwiązać kilka systemy na raz). Więc systemy są rozwiązane

Ściśle mówiąc, żaden z tych systemów nie ma rozwiązania. Jeżeli jednak macierz układu jest prostokątna, a liczba jego wierszy jest większa niż liczba kolumn (czyli liczba równań jest większa niż liczba niewiadomych), to taki układ nazywamy naddeterminowanym (układami naddeterminowanymi). a jego przybliżone rozwiązanie poszukuje się minimalizując normę euklidesową reszty. Tak więc dla pierwszego systemu rozwiązaniem jest x, który dostarcza minimum do następującego wyrażenia

(4-7x 1) 2 +(3-11x 1) 2 +(2-12x 1) 2

Łatwo jest upewnić się, że minimum dostarcza tylko . Użyjmy na przykład Symbolic Math Toolbox

Syms x R = (7*x-4)^2 + (11*x-3)^2 + (12*x-2)^2 dRdx = diff(R) x = rozwiąż(dRdx) x = 85/314 podwójne(x) ans = 0,2707

Tak więc wektor uzyskany przez napisanie x = f \ A zawiera rozwiązanie trzech układów naddeterminowanych, z których prawa strona każdego z nich jest odpowiednią kolumną macierzy A.

Oczywiście macierz przedefiniowanego systemu niekoniecznie zawiera tylko jedną kolumnę. Jako przykład rozważ problem z doborem parametrów a oraz b model liniowy zgodnie z podanymi danymi tabelarycznymi

k	1	2	3	4	5
x k	1.0	1.5	2.0	2.5	3.0
y k	2.99	2.81	2.89	3.03	3.21

Problem ten sprowadza się do nadmiernie określonego systemu w odniesieniu do a oraz b, jeśli wymagamy równości y(x k)=y k dla k=1,2,...,5:

Składamy macierz układu i wektor prawej strony

X = (1:0.5:3)" A = ;f = ; c(1)/x + c(2)*log(x); fplot(fun, ) trzymaj na wykresie(x, y, "lub")

Jeśli układ ma więcej równań niż niewiadomych, to taki układ nazywamy układem niedookreślonym i jest również rozwiązywany w MATLAB za pomocą znaku odwrotnego ukośnika. Będzie miał nieskończenie wiele rozwiązań i znajdzie rozwiązanie zawierające jak najwięcej wartości zerowych. Rozważ ten przykład

A = ; f = ; x = A\f x = 1,3333 1,0000 0 1,6667

Oprócz znaku odwrotnego ukośnika używany jest również normalny ukośnik. Łączy je następująca zasada

B/A jest równoważne (A"\B")"
gdzie apostrof oznacza transpozycję. Możesz także użyć funkcji mrdivide zamiast normalnego znaku ukośnika.

Na razie nie będziemy zagłębiać się w algorytmy rozwiązywania systemów wbudowanych w operację odwrotnego ukośnika. Ta sekcja jest poświęcona temu. W szczególności rozkład LU służy do rozwiązywania układów z ogólną macierzą kwadratową. Dodatkowo sprawdzana jest warunkowość macierzy. Klasycznym przykładem źle uwarunkowanej macierzy jest macierz Hilberta, której elementy określa wzór. Aby utworzyć macierz Hilberta o danym rozmiarze w MATLAB, użyj funkcji hilb. Rozwiążmy na przykład układ 13. rzędu, którego prawa strona jest taka, że ​​wszystkie jednostki muszą być jego rozwiązaniem ( k-ty element po prawej stronie to suma elementów k wiersz macierzy). A = hilb(13); f = suma(A, 2); x = A\f

W rezultacie otrzymujemy źle uwarunkowaną wiadomość

Ostrzeżenie: Matryca jest zbliżona do liczby pojedynczej lub źle przeskalowana. Wyniki mogą być niedokładne. RCOND = 2,339949e-018.

(RCOND to oszacowanie odwrotności numeru stanu) i błędna decyzja

X = 1,0000 1,0000 1,0004 0,9929 1,0636 0,6552 2,2023 -1,7860 5,3352 -3,4773 3,9431 -0,1145 1,1851

O niebezpieczeństwach, jakie niesie ze sobą zła warunkowość matrycy, napisano w rozdziale. Przejdziemy teraz do przeglądu rozszerzeń macierzy dostępnych w MATLAB.

Rozszerzenia macierzowe Cholesky, LU i QR

MATLAB posiada funkcje dla następujących rozszerzeń:

• Rozszerzenia Cholesky'ego - funkcja chol;
• Rozkłady LU - funkcja lu;
• Rozszerzenia QR - funkcja qr.

Rozszerzenia pomagają skutecznie rozwiązać cały zestaw układów liniowych równań algebraicznych z tą samą macierzą i kilkoma wektorami po prawej stronie. Rozważ następujący zestaw systemów:

Oś=f(k), k=1, 2,...,N.

Załóżmy, że macierz A jest reprezentowany jako iloczyn dwóch macierzy i systemów z macierzami A=BC i systemy z macierzami B oraz C rozwiązany znacznie szybciej niż w przypadku matrycy A. Następnie rozwiązania oryginalnych systemów uzyskuje się w następujący sposób. Ponieważ Oś=f(k) oraz A=BC, następnie BCx=f (k) i rozwiązywanie sekwencyjnie By=f(k) oraz Cx=y znaleźć rozwiązanie k systemu. Dla wypełnionych macierzy dekompozycja jest wykonywana w O(n 3) operacje, gdzie n- wielkość macierzy (tj. w tym samym czasie co rozwiązanie systemu) oraz rozwiązanie systemu z każdym współczynnikiem rozszerzenia w O(n 2) operacje. Dlatego decyzja N układy c z predekompozycji trwa O(n 3)+O(n 2)N. Rozwiązanie każdego układu bez rozkładu macierzy wymagałoby O(n 3)N działania arytmetyczne.

Rozkład Choleskiego dla danej macierzy A jest znalezienie takiej macierzy górnego trójkąta R z dodatnimi elementami przekątnymi, czyli A=RT R. Wiadomo, że jeśli macierz A jest symetryczny i dodatnio określony (tj. dla dowolnego wektora x jest prawdziwe: x T Ax0, czyli wszystkie wartości własne macierzy są dodatnie), to rozkład Cholesky'ego istnieje, a ponadto jest unikalny.

Rozważmy jako przykład rozkład macierzy Choleskiego

A = ; R = chol(A) R = 2,0000 0,5000 0,5000 0 1,9365 0,3873 0 0 1,8974

Proste sprawdzenie upewnia się, że rozwinięcie jest poprawne (w granicach błędów, które pojawiają się podczas operacji na liczbach rzeczywistych)

A - R"*R ans = 1.0e-015 * 0 0 0 0 0 0 0 0 0,4441

Jeśli macierz nie jest dodatnio określona, ​​stanie się to jasne w procesie dekompozycji, ponieważ będzie potrzeba wydobycia korzenia z Liczba ujemna i wyświetlane jest odpowiednie ostrzeżenie:

A = ; R = chol(A) ??? Błąd przy użyciu ==> chol Matrix musi być dodatnio określony.

Funkcja chol nie sprawdza macierzy pod kątem symetrii. Przy obliczaniu wykorzystuje się elementy macierzy oryginalnej znajdujące się na przekątnej i powyżej:

A = ; R = chol(A) R = 2,0000 0,5000 1,0000 0 1,9365 0,2582 0 0 1,7127 A - R"*R ans = 0 0 0 0 0 0 -1 0 0

Dla dowolnej macierzy kwadratowej można zrobić LU-rozkład, czyli znajdź dolną trójkątną macierz L i górna trójkątna matryca U takie, że A=LU. Bardziej rygorystyczne twierdzenie jest sformułowane w następujący sposób.

Jeśli Ak jest głównym minorem macierzy kwadratowej A rozmiar n, składający się z pierwszego k wiersze i kolumny (tj. A(1:k, 1:k)) oraz det( K)0 dla k=1,2,..., n-1, to istnieje unikalna dolna trójkątna macierz L, którego przekątna składa się z jedynek i jedynej górnej trójkątnej macierzy U takie, że A=LU.

Podczas obliczania LU Dekompozycja LU może wymagać permutacji wierszy, aby zapewnić stabilność numeryczną procesu dekompozycji, dzięki czemu funkcja lu może zwrócić macierz L, który jest trójkątny dolny aż do permutacji linii, na przykład:

A = ; = lu(A) L = 0 1.0000 0 1.0000 0 0 0.5000 0.5000 1.0000 U = 2 3 1 0 1 1 0 0 4

Praca LU gdzie L- trójkątne aż do permutacji wierszy i będzie równe macierzy A(ogólnie rzecz biorąc, biorąc pod uwagę błąd w operacjach na liczbach rzeczywistych).

A-L*U ans = 0 0 0 0 0 0 0 0 0

Jeśli dla matrycy A wykonywana jest dekompozycja Choleskiego lub dekompozycja LU, następnie rozwiązanie układu z macierzą A, jak wspomniano powyżej, sprowadza się do rozwiązania dwóch układów z macierzami trójkątnymi. Rozwiązanie to można wykonać za pomocą operacji odwrotnego ukośnika, ponieważ osadzony w nim algorytm określa macierze trójkątne i stosuje efektywny sposób rozwiązywania z nimi układu, wymagający O(n 2) operacje arytmetyczne. Rozważ przykład rozwiązania systemu

ze wstępnym LU rozkład macierzy.

A = ; f = ;

Wykonujemy LU-rozkład

Lu(A);

i rozwiązać kolejno dwa układy z macierzami trójkątnymi, najpierw z L, potem z U

Y = L\f; x = U\y x = 1 1 1

Rozwiązanie dwóch systemów można zapisać jednym wyrażeniem

X = U\(L\f)

a wynik będzie taki sam. Należy zwrócić uwagę na znaczenie stosowania nawiasów do określenia kolejności rozwiązywania układów z macierzami trójkątnymi. Wyrażenie bez nawiasów

X = U\L\f

doprowadzi do zupełnie innego wyniku, ponieważ najpierw U/L, co jest równoważne (jak omówiono powyżej) do rozwiązywania systemów z macierzą U i kolumnami L jako wektorami po prawej stronie. W efekcie otrzymuje się macierz, której każda kolumna jest rozwiązaniem odpowiedniego układu, następnie układ jest rozwiązywany z tą macierzą i wektorem prawej strony f. Ten proces oczywiście nie ma nic wspólnego z rozwiązaniem pierwotnego systemu.

Rozważ ostatni rozkład macierzy - QR-dekompozycja wykonywana przez funkcję qr. QR-dekompozycji można dokonać dla macierzy prostokątnych, dokładnie jeśli A macierz rozmiarów m na n, to jest macierz ortogonalna Q rozmiar m na m(tj. takie, że Q -1 = Q T) i macierz R rozmiar m na n ze wszystkimi elementami zerowymi pod główną przekątną, czyli A=QR.

QR-ekspansja charakteryzuje się dobrą stabilnością obliczeniową i jest stosowana w szczególności przy aproksymacji danych metodą najmniejszych kwadratów. Jako przykład rozważ rozwiązanie poprzedniego systemu z wykorzystaniem rozkładu QR:

Qr(A) Q = -0,7428 0,6000 -0,2971 -0,5571 -0,8000 -0,2228 -0,3714 0,0000 0,9285 R = -5,3852 -5,3852 -5,0138 0 -5,000 0,4000 0 0 6,6108

Po zakończeniu QR-rozkład, rozwiązanie układu Topór=f można znaleźć jako rozwiązanie trójkątnego systemu matrycy R, dlatego ponieważ QR=f, następnie Rx=Q T f:

>> x = R\(Q"*f) x = 1.0000 1.0000 1.0000

Aby uzyskać więcej informacji na temat rozwinięć macierzy odpowiadających funkcjom MATLAB chol, lu i qr oraz ich argumentom, zobacz sekcję Rozkłady macierzy.

Algorytm rozwiązywania układu równań liniowych za pomocą znaku odwrotnego ukośnika

Rozwiązując układ liniowych równań algebraicznych w MATLAB za pomocą znaku odwrotnego ukośnika
x = A\b
działa algorytm, który w zależności od rodzaju matrycy rozwiązuje system przy użyciu najbardziej odpowiedniej metody zaimplementowanej w jednej z procedur pakietu LAPACK lub UMFPACK. Tutaj b może być również macierzą, której liczba wierszy jest równa liczbie wierszy macierzy A. Następnie zwracana jest macierz x, każda i-ta jej kolumna zawiera rozwiązanie układu Ax (i) = b (i) , i=1,2,... ,k, gdzie b (i) jest i-tą kolumną macierzy b = [ b (1) | b(2) ...| b(k)].

Algorytm zaimplementowany w operacji \ składa się z następujących kroków:

• 1. W trywialnym przypadku, jeśli A jest sparse i diagonalne (sekcja Sparse Matrices jest poświęcona macierzom sparse w MATLAB), rozwiązanie znajduje się za pomocą prostego wzoru x k = b k /a kk , gdzie k=1,2,. ..n.

2. Jeśli A jest macierzą kwadratową, rzadką i pasmową, używany jest solwer macierzy pasmowej. Mogą być na to dwie opcje:

a. Jeżeli A jest macierzą trójkątną, a b jest jedną kolumną liczb rzeczywistych, to układ jest rozwiązywany przez eliminację Gaussa (jego operacje są wykonywane tylko na elementach na przekątnych). Jeśli podczas procesu eliminacji wymagane są permutacje wierszy w celu utrzymania stabilności lub jeśli macierz A nie jest trójprzekątna, działa następujący element.
b. Jeśli A jest pasmowy z gęstością elementów niezerowych większą niż parametr bandden, domyślnie równe 0,5, lub warunki z poprzedniego paragrafu nie są spełnione, to w zależności od typu A i b ( podwójnie lub pojedynczy) nazywane są następujące procedury biblioteki LAPACK:

A i b prawdziwy typ podwójne - procedury DGBTRF, DGBTRS
Podwójny typ złożony A i b - procedury ZGBTRF, ZGBTRS
A lub b pojedynczy typ rzeczywisty - procedury SGBTRF, SGBTRS
A lub b złożony typ pojedynczy - procedury CGBTRF, CGBTRS

Gęstość elementów niezerowych rozumiana jest jako stosunek liczby elementów niezerowych w taśmie do liczby wszystkich elementów w taśmie matrycowej np. dla matrycy 7 na 7, której szablon jest podany poniżej liczba elementów niezerowych

nz = 1 (wykres#-2) + 6 (wykres#-1) + 7 (wykres#0) + 6 (wykres#1) + 1 (wykres#2) + 1 (na wykresie nr 3) = 22,

i liczba wszystkich elementów w taśmie

pasmo = 5 (wykres#-2) + 6 (wykres#-1) + 7 (wykres#0) + 6 (wykres#1) + 5 (wykres#2) + 4 (na wykresie nr 3) = 33

a gęstość elementów niezerowych wyniesie 2/3. Ponieważ 2/3 > 0,5, zostanie zastosowany solwer macierzy pasm
Parametr bandden, podobnie jak inne parametry sterujące algorytmami macierzy rzadkich w MATLAB, ustawia się za pomocą funkcji spparms.



Znaczenie tego kroku algorytmu jest takie, że rozwiązanie systemu z macierzą pasmową nie wymaga działań poza pasmem. Jeśli taśma jest wystarczająco pełna, nie będzie zbyt wielu akcji z zerowymi elementami. Wręcz przeciwnie, jeśli taśma jest wystarczająco rzadka, to podejścia podane w dalszych krokach algorytmu mogą przynieść większy efekt.

• 3. Jeżeli A jest górną lub dolną macierzą trójkątną, to stosuje się metodę odwrotnego podstawienia lub odpowiednio metodę bezpośredniego podstawienia, w której składnik rozwiązania znajduje się z ostatniego (lub pierwszego równania), a następnie znalezione składniki są podstawione do poniższych równań, aby znaleźć kolejne rozwiązania składowe układów równań.
• 4. Jeśli A - można zredukować permutacjami do postaci trójkątnej, to wykonuje się odpowiednią redukcję, a następnie układ równań rozwiązuje się jak w pkt 3.

5. Jeżeli A jest macierzą symetryczną (lub w przypadku złożonym hermitowskim), to jej przekątna zawiera tylko dodatnie elementy rzeczywiste, to w zależności od tego, czy A jest rzadki, czy nie, spełniony jest punkt a) lub b).

a. Jeśli A nie jest rzadki, to podejmuje się próbę wykonania rozkładu Choleskiego A = R T R, a następnie rozwiązania układów z macierzami trójkątnymi R T i R: R T y = b i Rx = y. Wywołuje to następujące procedury biblioteki LAPACK:

A i b są prawdziwe i podwójne - DLANGE, DPOTRF, DPOTRS, DPOCON
A i b są złożone i podwójne - ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
A i b są prawdziwe i typu single - SLANGE, SPOTRF, SPOTRS, SDPOCON
A i b złożone i jednorodne - CLANGE, CPOTRF, CPOTRS, CPOCON

Rozkład Cholesky'ego nie powiedzie się, jeśli macierz nie będzie dodatnio określona. Ale nie przeprowadza się wstępnej kontroli na dodatnią jednoznaczność i określa się ją dokładnie w procesie dekompozycji, ponieważ dość często wystarczy wykonać kilka etapów rozkładu Choleskiego, aby dowiedzieć się, że macierz nie jest dodatnio określona ( podczas obliczeń konieczne staje się wyodrębnienie pierwiastka kwadratowego z liczby ujemnej). Jeśli rozkład Choleskiego się nie powiedzie, przechodzimy do następnego kroku.

b. Jeśli A jest rzadkie, to symetryczne permutacje wierszy i kolumn są wstępnie wykonywane zgodnie z algorytmem symetrycznego minimalnego stopnia (za pomocą funkcji symmmd) w celu zmniejszenia wypełnienia współczynnika rozszerzenia Choleskiego, tj. aby zmniejszyć liczbę nowych niezerowych elementów, które pojawiają się podczas procesu napełniania:

p = symmmd(A) - wektor p zawiera permutacje

R = chol(A(p,p));

i rozwiązano dwa systemy z czynnikami Cholesky'ego, pierwszy z transponowanym czynnikiem i odpowiednimi permutacjami w wektorze po prawej stronie

a drugi ze współczynnikiem rozszerzalności i wprowadzeniem składników rozwiązania do odpowiednich pozycji w wektorze x
x(p, :) = R \ y

• 6. Jeśli A jest macierzą Hessenberga, to zostaje sprowadzona do macierzy górnego trójkąta (z odpowiednią modyfikacją prawej strony), a następnie układ jest rozwiązywany przez podstawienia

7. Jeżeli A jest macierzą kwadratową, która nie spełnia warunków ust. 1-6, to w zależności od tego, czy jest rzadka, czy nie, wykonywane są następujące czynności

a. Jeżeli A nie jest macierzą rzadką, to za pomocą eliminacji Gaussa z wyborem elementu wiodącego (w celu zapewnienia stabilności rozwinięcia) dokonuje się dekompozycji LU macierzy A = LU, gdzie

U - górna trójkątna matryca
L - macierz, którą można sprowadzić do trójkątnej przez permutacje wierszy

a rozwiązanie układu Ax = b znajdujemy przez kolejne rozwiązania układów z macierzami trójkątnymi Ly = b, Ux = y.
Aby przeprowadzić dekompozycję LU, wywoływane są następujące procedury biblioteki LAPACK:

A i b są prawdziwe i podwójne - DLANGE, DGESV, DGECON
A i b są złożone i typu podwójne - ZLANGE, ZGESV, ZGECON
A i b są prawdziwe i typu single - SLANGE, SGESV, SGECON
A i b złożone i jednorodne - CLANGE, CGESV, CGECON

b. Jeśli A jest macierzą rzadką, kolumny są permutowane w celu zmniejszenia wypełnienia czynników L i U w procesie znajdowania rozkładu. Ponadto, permutując wiersze w procesie rozkładu LU, zapewniona jest stabilność obliczeniowa, po czym układy z macierzami trójkątnymi są ponownie rozwiązywane. Aby przeprowadzić dekompozycję LU macierzy rzadkiej, stosuje się procedury biblioteki UMFPACK

8. Jeśli poprzednie punkty algorytmu nie zadziałały, a zatem A nie jest macierzą kwadratową, to wszystko jest określone przez jej rzadkość i działa jeden z następujących punktów:

a. Jeśli A nie jest macierzą rzadką, to przeprowadzany jest rozkład QR AP = QR, gdzie P jest macierzą permutacji kolumnowej, Q jest macierzą ortogonalną (Q T Q = I), a R jest górnym trójkątem. Jeśli A to m na n, to Q to m na m, a R to m na n. Ponadto rozwiązanie systemu jest następujące:

x = P*(R \ (Q" * b))

ponieważ z Ax = b i AP = QR wynika to: QRx = bP i Rx = Q T bP.

b. W przypadku rzadkiej i prostokątnej macierzy A nie ma sensu obliczanie macierzy Q, ponieważ najprawdopodobniej zostanie ona wypełniona. Dlatego algorytm rozkładu QR oblicza c = Q T b (tj. zmodyfikowaną prawą stronę) i rozwiązuje układ Rx = c z macierzą trójkątną, aby uzyskać rozwiązanie oryginalnego układu Ax = b.

Oprócz wszystkich powyższych, powyższy algorytm ocenia warunkowość macierzy i wyświetla ostrzeżenie o możliwym błędzie w rozwiązaniu, jeśli macierz jest źle uwarunkowana (patrz rozdział Wpływ warunkowości macierzy na dokładność rozwiązywania system z nim). Powyższy algorytm zawiera wszelkiego rodzaju kontrole, które mogą zająć znaczną ilość dodatkowego czasu. Poniższa sekcja, Funkcja linsolve do rozwiązywania układów równań liniowych, zawiera funkcję linsolve, która pozwala określić typ macierzy, skracając w ten sposób czas obliczeń.

funkcja linsolve do rozwiązywania układów równań liniowych

Zamiast rozwiązywać układ (lub układy z wieloma prawymi stronami określonymi w macierzy) liniowych równań algebraicznych za pomocą znaku odwrotnego ukośnika, można użyć funkcji linsolve, która nie wykonuje wszystkich kontroli macierzy właściwych algorytmowi \ operacja (patrz poprzedni rozdział).

Funkcja linsolve, wywoływana w najprostszej formie z dwoma argumentami wejściowymi i jednym argumentem wyjściowym
x = linsolve(A, b)
rozwiązuje układ Ax = b na jeden ze sposobów, w zależności od tego, czy macierz jest kwadratowa, czy nie.

• Jeśli A jest macierzą kwadratową, to jej rozkład LU jest wstępnie obliczany, a następnie rozwiązywane są dwa układy z macierzami trójkątnymi L i U.
• Jeśli A jest macierzą prostokątną, to jej rozkład QR jest wstępnie obliczany, a następnie rozwiązywany jest układ z macierzami trójkątnymi.

(Więcej szczegółów na temat działań z wynikającymi z nich współczynnikami ekspansji znajduje się w sekcji). W takim przypadku, jeśli macierz A jest źle uwarunkowana lub A jest macierzą o niepełnej randze, to wyświetlane jest odpowiednie ostrzeżenie, na przykład:

A = ; b = ; x = linsolve(A,b) Ostrzeżenie: brak rangi, ranga = 1, tol = 4.4686e-015. x = 0 0 2.0000

Aby wyłączyć wyjście tej wiadomości w oknie poleceń, wywołaj funkcję linsolve z dwoma argumentami wyjściowymi
= linsolve(A, b)
wtedy otrzymane rozwiązanie zostanie zapisane w x, a albo rząd macierzy (jeśli A jest macierzą prostokątną) albo odwrotność oszacowania dla numeru jej warunku (patrz rozdział ), na przykład

A = ; b = ; = linsolve(A,b) x = -1.0000e+000 9.9999e-001 3.3307e-006 r = 6.9444e-013

Główne zalety funkcji linsolve nad operacją \ pojawiają się przy określaniu typu macierzy układu równań. W tym celu przed wywołaniem funkcji linsolve należy wypełnić strukturę informacjami o macierzy z następującymi polami

• SYM - symetryczny;
• LT - dolny trójkątny;
• UT - górny trójkątny;
• UHESS - Hessenberg;
• POSDEF - symetryczny, dodatnio określony;
• RECT - prostokątny;
• TRANSA - czy konieczne jest rozwiązanie układu z macierzą transponowaną na daną.

Każde pole może mieć wartość prawda lub fałsz. Oczywiście nie wszystkie kombinacje wartości pól są dozwolone, np. macierz nie może być trójkątna i jednocześnie symetryczna i dodatnio określona. Prawidłowe kombinacje wartości pól zebrano w poniższej tabeli

LT	UT	UHESS	SYM	POSDEF	RECT	TRANSA
PRAWDA	fałszywy	fałszywy	fałszywy	fałszywy	prawda fałsz	prawda fałsz
fałszywy	PRAWDA	fałszywy	fałszywy	fałszywy	prawda fałsz	prawda fałsz
fałszywy	fałszywy	PRAWDA	fałszywy	fałszywy	fałszywy	prawda fałsz
fałszywy	fałszywy	fałszywy	PRAWDA	PRAWDA	fałszywy	prawda fałsz
fałszywy	fałszywy	fałszywy	fałszywy	fałszywy	prawda fałsz	prawda fałsz

Jeśli macierz układu jest dodatnio określona, ​​należy to wziąć pod uwagę przy rozwiązywaniu, ponieważ w przypadku macierzy dodatnio określonych rozwiązanie opiera się na rozkładzie Choleskiego, który wymaga mniej operacji niż rozkład LU stosowany przy rozwiązywaniu układów z ogólnymi macierzami kwadratowymi . Łatwo to zweryfikować na poniższym przykładzie, w którym tworzona jest symetryczna macierz dodatnio określona (do transponowanej do niej dodawana jest macierz liczb losowych, a do przekątnej dodawane są duże liczby) oraz układ równań z tą macierzą najpierw rozwiązywany jest jako układ z macierzą ogólną (opts.SYM=false i opts.POSDEF=false), a następnie jak z macierzą symetryczną i dodatnio określoną (opts.SYM=true i opts.POSDEF= PRAWDA).

% ustawia wszystkie pola struktury ops, z wyjątkiem SYM i POSDEF opts.TRANSA = false; opts.UHESS = fałsz; opts.RECT ​​= fałsz; opts.UT = fałsz; opts.LT = fałsz; % utwórz wektor o rozmiarach macierzy N = 2.^(8:12); % utwórz puste tablice do rejestrowania czasów rozwiązania TimeGeneral = ; CzasPosSym = ; % utwórz macierze w pętli i porównaj czasy rozwiązania dla n = N % utwórz macierz symetryczną i dodatnio określoną % oraz wektor prawej strony A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % rozwiąż układ jako układ z ogólną macierzą opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % rozwiąż system jako system z symetrią i macierzą pozycji opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve( A,b, opts); Tend = cputime; TimePosSym = ; end

Wynikające z tego koszty obliczeniowe tych sposobów rozwiązania systemu przedstawiono na poniższym wykresie.

Oczywiście, jeśli macierz jest trójkątna, to bardzo ważne jest, aby to wskazać, ponieważ rozwiązanie układu z macierzą trójkątną wykonuje się w operacjach O(n 2), a jeśli algorytm rozwiązania przeznaczony dla macierz ogólna jest stosowana do układu z macierzą trójkątną, a następnie kolejność operacji O(n 3).

Funkcja linsolve nie sprawdza, czy macierz spełnia określone właściwości, co może skutkować błędem. Na przykład, jeśli rozwiązując układ z następującą macierzą i prawą stroną

A = ; b = ; ustaw pole UT na true (i wszystkie inne ustaw na false) opts.UT = true; opts.TRANSA = fałsz; opts.LT = fałsz; opts.UHESS = fałsz; opts.SYM=fałsz; opts.POSDEF = fałsz; opts.RECT ​​= fałsz; wtedy przy rozwiązywaniu układu funkcja linsolve uwzględni macierz jako górną trójkątną i wybierze potrzebne elementy z górnego trójkąta A x = linsolve(A,b, opts) Da to rozwiązanie x = 1 1 1 odpowiadające macierzy A = ;

Wpływ warunkowości macierzy na dokładność rozwiązywania z nią układu

W tej części rozważamy, jak błędy w elementach wektora prawostronnego oraz w macierzy układu równań liniowych Ax = b mogą wpływać na rozwiązanie tego układu. Rozważmy najpierw przypadek wprowadzenia perturbacji do prawostronnego wektora b. Więc rozwiązujemy dwa systemy

ponadto zakłada się, że rozwiązujemy dokładnie każdy z układów, a głównym pytaniem, które się nasuwa, jest to, jak bardzo rozwiązanie x układu (1) będzie się różnić od rozwiązania układu (2) z zaburzeniem po prawej stronie δb. Jest to dość ważne pytanie, ponieważ elementy prawej strony można otrzymać w wyniku niektórych pomiarów, zawierających odpowiednio błąd δb k w każdej składowej b k . Pożądane byłoby, aby przy pomiarze tej samej wielkości (za każdym razem z własnymi, małymi błędami) odpowiednie rozwiązania układów z nieco różnymi prawymi stronami również nie różniły się zbytnio od siebie. Niestety, nie zawsze tak jest. Powodem tego jest charakterystyka macierzy, zwana jej kondycjonowanie. Zostanie to omówione dalej.

Najpierw musimy wprowadzić miarę bliskości wektorów i x, czyli wektor błędu . Miarą wielkości wektora jest norma (można ją definiować na różne sposoby). Na razie weźmy za normę zwykłą normę euklidesową wektora (pierwiastek kwadratowy z sumy kwadratów jego składowych), tj.

Odpowiednio

gdzie n to liczba niewiadomych i równań w systemie. Oprócz normy wektorowej do szacowania rozmiaru wektora, potrzebujemy również normy macierzy do szacowania rozmiaru macierzy. Weźmy normę macierzową, która jest zdefiniowana następująco (nazywa się to normą spektralną):

tych. norma macierzy widmowej jest równa pierwiastkowi kwadratowemu z maksymalnej wartości własnej macierzy A T A. MATLAB ma funkcję normy do obliczania norm macierzy i wektorów, która w szczególności może obliczać normy podane powyżej. Dlaczego wybraliśmy takie normy wektorów i macierzy, wyjaśniono szczegółowo w sekcji, w której podano niektóre obliczenia i definicje. Wiąże się to z oszacowaniem, którego użyjemy dla błędu w rozwiązaniu systemu (wyprowadzenie tej nierówności jest również podane we wspomnianym punkcie):

Tutaj oznacza numer warunku macierzy, który definiuje się następująco:

Z powyższej nierówności wynika, że ​​im większa liczba warunków macierzy systemu, tym większy może być błąd względny w rozwiązaniu z małym błędem po prawej stronie .

Rozważ klasyczny przykład źle uwarunkowanej macierzy – macierz Hilberta – i dowiedz się, jak błąd po prawej stronie układu wpływa na błąd w rozwiązaniu. Macierz Hilberta jest zdefiniowana w następujący sposób

Aby utworzyć macierz Hilberta, MATLAB udostępnia funkcję hilb, której argument wejściowy określa rozmiar macierzy. Weźmy małą macierz 6 na 6:

N = 6; H = hilb(n) H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.090

X = jedynki (n, 1); b = H*x b = 2,4500 1,5929 1,2179 0,9956 0,8456 0,7365

Widzimy, że ani w macierzy, ani w wektorze po prawej stronie nie ma nic „podejrzanego”, wszystkie liczby niewiele się od siebie różnią. Teraz tworzymy zaburzoną prawą stronę b + δb, dodając małe liczby rzędu 10 -5 do wektora b i rozwiązujemy układ z zaburzoną prawą stroną, aby otrzymać wektor .

delta_b = 1e-5*(1:n)"; x_tilda = H\(b + delta_b) x_tilda = 0,9978 1,0735 0,4288 2,6632 -1,0160 1,8593

Widać, że wynikowe rozwiązanie jest dalekie od dokładnego, w którym powinny znajdować się wszystkie jednostki. Obliczmy błąd względny w rozwiązaniu i po prawej stronie (funkcja norm domyślnie oblicza normę euklidesową wektora):

Delta_x = x - x_tylda; LEWY = norm(delta_x)/norm(x) LEWY = 1.1475 PRAWY = norm(delta_b)/norm(b) PRAWY = 2.7231e-005

Tak więc błąd jest w rozwiązaniu rzędu jednego, chociaż zmiany po prawej stronie były rzędu 10 -5 . To w pełni pasuje do powyższej nierówności błędu. Rzeczywiście, obliczamy numer warunku cond(H) za pomocą funkcji MATLAB zwanej cond i domyślnie oblicza dla normy widmowej macierzy

C = przew(H) c = 1,4951e+007

LEWY 1.1475 mniej

C*PRAWO? 1.4951e+07 * 2.7231e-05 ≈ 407

a nierówność jest zaspokojona (nawet z pewnym marginesem).

Wraz ze wzrostem wielkości macierzy Hilberta błąd w rozwiązaniu będzie tylko rósł (łatwo to sprawdzić, ustawiając n = 7, 8, ...). Ponadto dla n = 12 wyświetli się komunikat, że matryca jest źle uwarunkowana i rozwiązanie może okazać się niepoprawne

Ostrzeżenie: Matryca jest zbliżona do liczby pojedynczej lub źle przeskalowana.
Wyniki mogą być niedokładne. RCOND = 2.409320e-017.

RCOND jest tutaj wybrany jako miara warunkowości. równy jeden podzielone przez oszacowanie numeru warunku (liczba warunku jest szacowana przy użyciu dość szybkiego algorytmu, który jest znacznie szybszy niż cond, jak opisano bardziej szczegółowo w pomocy funkcji rcond). Jeśli wartość RCOND jest niewielka, to matrycę uważa się za źle uwarunkowaną.

Wzrost błędu w rozwiązaniu wraz ze wzrostem wielkości macierzy Hilberta wynika z faktu, że liczba warunków macierzy Hilberta rośnie bardzo szybko wraz z jej rozmiarem. Łatwo to zweryfikować za pomocą prostej pętli i funkcji semilogii (skala wzdłuż osi y jest logarytmiczna):

N = 1:20; C = zera (1, 20); dla n = NH = hilb(n); C(n) = kond(H); koniec semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Nawiasem mówiąc, ponieważ funkcja cond znajduje numer warunku za pomocą metody numerycznej (a mianowicie dekompozycji na wartości osobliwe w celu znalezienia wartości osobliwych), liczba warunku po n = 12 jest już obliczana niepoprawnie, w rzeczywistości powinna dalej rosnąć, jak może być widocznym za pomocą obliczeń symbolicznych w MATLAB i operacji z podaną liczbą cyfr znaczących

N = 1:20; C = zera (1, N); cyfry(60) dla n = N H = vpa(sym(hilb(n))); % obliczenie macierzy Hilberta do 60 cyfry sigma = svd(H); % znajdź wartości osobliwe macierzy Hilberta % oblicz numer stanu macierzy Hilberta C(n) = max(double(sigma))/min(double(sigma)); koniec semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Rozważmy teraz trzy ważne punkty dotyczące tego przykładu.

Pierwszy z nich - rozwiązanie układu Hx = b (z wektorem po prawej stronie odpowiadającym rozwiązaniu jednostek) za pomocą operatora odwrotnego ukośnika nie da dokładnych jednostek, błąd będzie już w dziesiątej cyfrze (choć w MATLAB wszystkie obliczenia są wykonywane domyślnie z podwójną precyzją)

Format długi H\b ans = 0.99999999999916 1.00000000002391 0.99999999983793 1.00000000042209 0.9999999995366 1.00000000018389

Wynika to z faktu, że obliczając wektor b przy mnożeniu macierzy H przez wektor wszystkich jednostek, wprowadziliśmy już do niej pewien błąd. Dodatkowo pewną rolę odegrały również błędy zaokrąglania w procesie rozwiązywania systemu, a zła warunkowość macierzy (nawet niewielki rozmiar) prowadziła do takich błędów w rozwiązaniu. Rzeczywiście, dla małych macierzy o małej liczbie warunków efekt ten nie będzie obserwowany. Weź macierz liczb losowych 6 na 6, oblicz jej numer warunku

A = rand(n); przew(A) odp = 57,352452799075771

Następnie tworzymy prawą stronę odpowiadającą dokładnemu rozwiązaniu ze wszystkich jednostek

X = jedynki (n, 1); b = A*x;

i rozwiązać system, co skutkuje dobrą dokładnością

>> A\b ans = 1.00000000000000 1,000000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000

Drugi ważny punkt związane z wyznacznikiem macierzy. Oblicz wyznacznik macierzy Hilberta 6 na 6, za pomocą której rozwiązaliśmy układ

>> det(hilb(6)) ans = 5,3673e-018

Wyznacznik okazuje się na tyle mały, że można z tego wyciągnąć błędny wniosek, że źródłem problemu dużego błędu w rozwiązaniu z małym zaburzeniem prawej strony układu jest małość wyznacznika. W rzeczywistości tak nie jest, jeśli wyznacznik układu jest mały, to nie oznacza to, że macierz jest źle uwarunkowana i w konsekwencji występuje duży błąd w rozwiązaniu. Rzeczywiście, weź macierz

A = ;

Jej wyznacznik jest bardzo mały, rzędu 10 -121

Det(A) ans = 9.9970e-121

(słowo „bardzo” jest oczywiście dowolne, ale jest mniejsze niż wyznacznik macierzy Hilberta 6 na 6, przy rozwiązywaniu układu, z którym mieliśmy problemy). Jednak małość wyznacznika nie wpłynie na błąd w rozwiązaniu, gdy prawa strona układu jest zaburzona, co łatwo wykazać tworząc układ ze znanym rozwiązaniem, wprowadzając perturbacje po prawej stronie i rozwiązując układ:

X = ; b = A*x; delta_b = 1e-5*; x_tilda = A\(b+delta_b); delta_x = x - x_tylda; LEWY = norm(delta_x)/norm(x) PRAWY = norm(delta_b)/norm(b) LEWY = 2.1272e-005 PRAWY = 2.1179e-005

Tak więc błąd względny w rozwiązaniu jest prawie równy błędowi względnemu po prawej stronie, ponieważ numer warunku powyższej macierzy A jest nieco większy niż 1, a mianowicie:

C = przew(A) c = 1,0303

I na koniec rozważmy trzecie pytanie, dotyczące osiągnięcia znaku równości w nierówności dla błędu w rozwiązaniu

Innymi słowy, dowiedzmy się: czy może zaistnieć taka sytuacja, gdy zrobimy małe perturbacje względne po prawej stronie układu, powiedzmy 10 -5 , numer warunku macierzy układu wynosi​​10 10 , aw roztworze otrzymuje się błąd względny 10 5 . Możesz przejść przez wiele przykładów, iterując po różnych macierzach, nie uzyskując równości i błędnie dochodząc do wniosku, że jest to tylko przeszacowanie błędu w rozwiązaniu. Tak jednak nie jest, o czym przekonuje poniższy przykład.

w którym względna perturbacja prawej strony wynosi 10 -5

PRAWO = norma(delta_b)/norm(b) PRAWO = 1.0000e-005

Błąd względny w rozwiązaniu układu okazuje się wynosić 10 5

X = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tylda; LEWO = norma(delta_x)/norm(x) LEWO = 1.0000e+005

A dzieje się tak, ponieważ
1) w tym przykładzie numer warunku macierzy A wynosi 10 10 ;

C = przew(A) c = 1.0000e+010

2) nierówność błędu w rozwiązaniu zamienia się w równość.

Jeśli ustawimy format na long, zobaczymy, że LEFT, RIGHT i numer warunku nie są odpowiednio równe 10 5 , 10 -5 i 10 10 , ale jest to spowodowane błędami zaokrągleń. Po rozwiązaniu w arytmetyce dokładnej równość okazałaby się dokładna, co można pokazać analitycznie (patrz rozdział ).

W poprzednich przykładach założyliśmy, że macierz nie zawiera błędów, a błędy mogą znajdować się tylko po prawej stronie układu równań liniowych. W praktyce może się okazać, że elementy macierzy systemu są określone z błędami, tj. zamiast systemu

|| Topór || V = || b || V ⇒ || || M || x || V ≥ || b || V

gdzie || || M - jakaś norma macierzowa. Aby móc to zrobić, norma macierzowa || || M i norma wektorowa || || V musi spełniać następującą nierówność

|| || M || x || V ≥ || Topór || V

dla dowolnych macierzy A i wektorów x. Jeśli ta nierówność się utrzymuje, to norma macierzowa || || M nazywa się Zgoda z normą wektorową || || V. Wiadomo na przykład, że norma spektralna

(pierwiastek kwadratowy z maksymalnej wartości własnej macierzy AT) jest zgodny z normą wektora euklidesowego

Dlatego w dziale przeprowadziliśmy wszystkie eksperymenty z tymi normami.

Po uzyskaniu nierówności || || M || x || V ≥ || b || V dalej zauważyć, że od Ax = b i wynika z tego . Ponieważ macierz A jest nieosobliwa, wynika z tego, że δx = A -1 δb oraz || δx || V = || A-1 δb || V. Ponownie wykorzystujemy zgodność norm i uzyskujemy nierówność || δx || V ≤ || A-1 || M || b || V. Co więcej, w dwóch uzyskanych nierównościach

|| || M || x || V ≥ || b || V i || δx || V ≤ || A-1 || M || b || V

dzielimy mniejszą wartość jednej z nierówności przez większą wartość drugiej nierówności, a odpowiednio większą przez mniejszą

i poprzez proste przekształcenie w końcu uzyskujemy wymaganą nierówność

gdzie przew(A) = || || M * || A-1 || M.

Zauważ, że wyprowadzając nierówność dla błędu, wykorzystaliśmy tylko fakt, że norma macierzowa jest zgodna z wektorową. Dotyczy to nie tylko normy spektralnej macierzy i normy euklidesowej wektora, ale także innych norm. Na przykład maksymalna norma macierzy kolumn, obliczona ze wzoru

zgodny z pierwszą normą wektorową

oraz maksymalna norma macierzy wierszy, obliczona ze wzoru

zgodny z normą wektorów nieskończonych

Funkcja norm oblicza nie tylko normę euklidesową wektora i normę spektralną macierzy, ale także normy wektorowe i macierzowe wymienione powyżej. Aby to zrobić, musisz zadzwonić za pomocą dodatkowa sekunda argument:

• q = norm(A, 1) - maksymalna norma kolumny macierzy A
• q = norm(A, inf) - maksymalna norma wiersza macierzy A
• q = norm(a, 1) - pierwsza norma wektora a
• q = norm(a, inf) - nieskończona wektorowa norma a

Numer warunku macierzy cond(A) = || || M * || A-1 || M w odniesieniu do różnych norm macierzowych można obliczyć za pomocą funkcji cond. Jeśli cond zostanie wywołane z jednym argumentem wejściowym (macierzą), to liczba warunku jest obliczana w odniesieniu do normy macierzy widmowej. Przy wywołaniu funkcji cond z dodatkowym argumentem zwracane są numery warunków w odniesieniu do określonej normy macierzy:

• c = cond(A, 1) - numer warunku względem maksymalnej normy kolumny macierzy
• c = cond(A, inf) - numer warunku względem maksymalnej normy wiersza macierzy

Jako przykład pokazujący znaczenie zgodności dziury macierzy z normą wektora w nierówności błędu podamy przykład z macierzą A, wektorem prawostronnym b i wektorem błędu po prawej stronie delta_b .

A = ; b = [ -5,7373057243726833e-001 -1,5400413072907607e-001 -5,3347697688693385e-001 -6.0209490373259589e-001]; delta_b = [-0,71685839091451e-5 0,54786619630709e-5 0,37746931527138e-5 0,208503222383081e-5];

Obliczamy prawo i lewa strona estymacje z wykorzystaniem pierwszej normy wektorowej oraz liczby warunku macierzy względem normy macierzy widmowej

PRAWY = norma(delta_b,1)/norm(b,1); c = przew(A); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tylda; LEWO = norma(delta_x,1)/norm(x,1); sformatuj krótki e disp()

Otrzymujemy następujące wartości dla lewej i prawej części nierówności
9.9484e+004 9.9323e+004

Lewa strona nie powinna przekraczać prawej strony (jak wykazano powyżej), ale okazała się większa ze względu na wybór niespójnych norm.

Przeanalizujmy teraz, dlaczego liczba warunku macierzy nie może być mniejsza niż jeden. Numer warunku macierzy A jest zdefiniowany jako cond(A) = || || M * || A-1 || M , gdzie || || M - jakaś norma macierzowa A. Norma macierzowa (czyli reguła, według której każdej macierzy przypisywana jest liczba) nie może być dowolna, musi spełniać następujące cztery aksjomaty:

A1. || || ≥ 0 dla dowolnej macierzy A i || || = 0 wtedy i tylko wtedy, gdy A jest macierzą zerową.

A2. || αA || = | α | * || || dla dowolnej macierzy A i liczby α.

A3. || A+B || ≤ || || + || B || dla dowolnych macierzy A i B

A4. || AB || ≤ || || * || B || dla dowolnych macierzy A i B.

Z ostatniego aksjomatu wynika, że ​​norma jest definiowana tylko dla macierzy kwadratowych (choć w zasadzie w powyższych wzorach na obliczanie różnych norm nie ma takiego ograniczenia). Ponadto z ostatniego aksjomatu wynika, że ​​każda norma macierzy tożsamości I jest nie mniejsza niż jeden, w istocie

|| || = || I*I || ≤ || || 2 || || ≥ 1.

Następnie, ponownie korzystając z czwartego aksjomatu, otrzymujemy, że liczba warunku macierzy jest zawsze większa niż jeden (prawda dla liczby warunku macierzy w odniesieniu do dowolnej normy macierzy)

1 ≤ || || = || AA-1 || ≤ || || * || A-1 || = przew(A).

Dokończmy ta sekcja badając przyczynę pojawienia się dokładnej równości w szacowaniu błędu w rozwiązaniu pod perturbacją prawej strony:

i budowanie odpowiednich przykładów w MATLAB. (Rozumowanie podane poniżej zawarte jest np. w książce J. Forsythe, K Mohler. Numeryczne rozwiązywanie układów liniowych równań algebraicznych. M: "Mir", 1969.)

Istotną rolę w rozumowaniu odgrywa twierdzenie o rozkładzie macierzy na wartości osobliwe, zgodnie z którym dla każdej rzeczywistej macierzy o rozmiarze on istnieją dwie macierze ortogonalne U i V o rozmiarze n przez n (U T U=UU T i V T V = VV T) tak, że iloczyn D = U T AV jest macierzą diagonalną i możemy wybrać U i V tak, że

gdzie μ 1 ≥ μ 2 ≥… ≥ μ r ≥ μ r+1 =…= μ n =0,

a r jest rządem macierzy A. Liczby μ k nazywane są liczbami widmowymi macierzy A. W przypadku macierzy nieosobliwej A prawdziwe jest:

μ 1 ≥ μ 2 ≥ … ≥ μ n ≥ 0.

Kolejnym ważnym faktem jest to, że mnożenie przez macierz ortogonalną nie zmienia normy euklidesowej wektora, tj. dla dowolnego wektora x o n elementach i dowolnej macierzy ortogonalnej U o rozmiarze n na n mamy

|| Ux || = || x ||.

Skoro mnożenie przez macierz ortogonalną nie zmienia normy widmowej, to

dlatego dla warunku numer macierzy równość jest prawdziwa

Mamy dwa systemy: Ax = b (z dokładne rozwiązanie x) i (z dokładnym rozwiązaniem). Jest oczywiste, że błąd δx jest rozwiązaniem układu, którego prawą stroną jest zaburzenie δb, czyli systemy Aδx = δb. Niech D = U T AV będzie osobliwym rozkładem macierzy A, to UDV T = A ze względu na to, że U i V są macierzami ortogonalnymi. Dalej

Ax = b ⇒ UDV T x = b.

Wprowadźmy notację

x" = V T x, b" = U T b,

wtedy następujące systemy są równoważne

Oś = b ⇔ Dx" = b"

Dokładnie w ten sam sposób rozpatrujemy system pod kątem błędu

Aδx = δb ⇒ UDV T δx = δb

Wprowadzamy notację

δx" = V T δx, δb" = U T δb,

w których następujące systemy są równoważne

Aδx = δb ⇔ Dδx" = δb"

Tych. otrzymaliśmy proste układy równoważne z macierzą diagonalną, której przekątna zawiera liczby widmowe macierzy A.

Wybierzmy teraz w szczególny sposób prawe strony tych systemów, a mianowicie niech

gdzie β > 0, to odpowiednim rozwiązaniem układu Dx" = b" będzie

gdzie μ 1 jest maksymalną wartością osobliwą macierzy A. Podobnie postąpimy z układem Dδx" = δ b" , czyli niech

gdzie γ > 0, to odpowiednim rozwiązaniem układu Dδx" = δb" będzie

Ponieważ norma wektora jest zachowana, gdy jest pomnożona przez macierz ortogonalną, wynika z tego, że

β/μ 1 = || x" || = || V T x || = || x || i γ/μ n = || δx" || = || V T δx || = ||δx ||.

W dokładnie ten sam sposób uzyskujemy równości:

β = || b" || = || U T b || = || b || i γ = || δb" || = || U T δb || = || b ||.

i ponieważ

wtedy w końcu dostajemy

Tak więc dla każdej macierzy A można skonstruować wektor prawej strony układu i jego zaburzenia tak, aby błąd w rozwiązaniu był iloczynem liczby warunku macierzy i błędu po prawej stronie. W MATLAB odpowiednia konstrukcja może być wykonana bez użycia dekompozycji według wartości osobliwych (chociaż można ją znaleźć za pomocą funkcji svd).

Najpierw określamy n i otrzymujemy dwie macierze ortogonalne U i V, wykonując rozkład QR macierzy n przez n z liczb losowych:

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

Następnie konstruujemy macierz A używając macierzy diagonalnej D oraz macierzy ortogonalnych U i V

A=U*D*V”;

Numer warunku macierzy A będzie pokrywał się z numerem warunku macierzy D, który w naszym przykładzie wynosi 10 10

beta=1; gamma = 1e-5;

i skonstruuj wektory b" i δb"

B1 = "; db1 = ";

dzięki której znajdujemy wektory b i δb

X = A\b; x_tilda = A\(b+db);

obliczyć lewą i prawą część nierówności

dx = x - x_tylda; PRAWO = norma(db)/norm(b); W LEWO = norma(dx)/norm(x);

i wyjmij je

sformatuj krótki e disp()

Dostajemy równość

1. Program można uruchomić za pomocą ikony na pulpicie lub poprzez przycisk Start (w lewym dolnym rogu ekranu). Otworzy się przestrzeń robocza pakietu MatLab.

2. Kliknij w polu Okno poleceń, aby je uaktywnić.

Wpisz 1+2 w wierszu z ikoną » i migającym pionowym kursorem i naciśnij Enter. W rezultacie w oknie poleceń MatLaba wyświetlane jest:

Wynik obliczenia sumy 1+2 jest zapisywany do zmiennej specjalnej ans i jej wartość równa 3 jest wyświetlana w oknie poleceń. Pod odpowiedzią znajduje się linia poleceń z migającym kursorem, co oznacza, że ​​MatLab jest gotowy do dalszych obliczeń. Możesz wpisać nowe wyrażenia w wierszu poleceń i znaleźć ich wartości.

3. Aby kontynuować pracę z poprzednim wyrażeniem, na przykład do obliczenia (1 + 2) / 4,5, wskazane jest użycie już istniejącego wyniku, który jest przechowywany w zmiennej ans.

Wpisz ans/4.5 (kropka jest używana przy wprowadzaniu liczb dziesiętnych) i naciśnij Enter, okazuje się:

4. Rekord „ans = 0,6667” nazywany jest echem.

Wykonywaniu każdej komendy w MatLabie towarzyszy podobne echo, co często utrudnia dostrzeżenie działania programu.

Aby wyłączyć echo, wstaw znak po poleceniu; (średnik). Na przykład:

Tutaj wynik pomnożenia ans przez 3, który jest pośredni, nie jest wyświetlany. Wyświetlana jest tylko ostateczna odpowiedź.

5.Zapisz wartości zmiennych. Dla tego:

- wybierz Zapisz obszar roboczy jako z menu Plik;

— w wyświetlonym oknie Save Workspace Variables podaj katalog i nazwę pliku (domyślnie sugerowane jest zapisanie pliku w podkatalogu work głównego katalogu MatLaba). Wyniki pracy zostaną zapisane w pliku z rozszerzeniem maty.

6. Zamknij MatLab.

7. Uruchom ponownie MatLab. Aby przywrócić wartości zmiennych z poprzedniej sesji, otwórz zapisany plik za pomocą podpunktu Otwórz w menu Plik. Zapisane zmienne mogą być używane w nowo wprowadzanych poleceniach.

8. Zapisz polecenia, które wykonujesz, i wyniki w plik tekstowy, który następnie można odczytać lub wydrukować z Edytor tekstu. Dla tego:

- wprowadź polecenie dziennika;

- określ nazwę pliku, w którym będzie przechowywany dziennik pracy jako argument polecenia dziennika.

Przykład podano w paragrafie 1.3.

9. Aby zakończyć pracę w systemie MatLab, wpisz komendę quit.

1. Przestudiuj część teoretyczną.

2. Uzyskaj opcję zadania.

3. Wykonaj przykład podany w paragrafie 2.

4. Wykonaj obliczenia zgodnie z wybraną opcją.

5. Prześlij zgłoszenie w w formie elektronicznej.

6. Chroń Praca laboratoryjna odpowiadając na pytania nauczyciela.

Opcje

Artykuły do ​​przeczytania:

Podstawy pracy z Mathcadem. Najprostsze obliczenia. Lekcja 4