Rozwiązywanie układów równań

Równania o dwóch niewiadomych

Przykład 1: Rozwiąż układ równań liniowych

7 x -  4 y  =  3
2 x + 2 y = -3

Metody graficzne i numeryczne

Na początek zapiszemy dwa równania do zmiennych:

rownanie1=7x-4y==3

rownanie2=2x+2y==-3

Układ równań to lista pojedynczych równań:

uklad={rownanie1,rownanie2}

Graficzne rozwiązanie polega na narysowaniu dwóch funkcji uwikłanych na jednym wykresie. Współrzędne punktu przecięcia się linii to rozwiązanie układu 2 równań:

ContourPlot[uklad,{x,-12,12},{y,-12,12}]

I tu mała niespodzianka: pusty wykres. Początkującego użytkownika to może zmylić.  Aby ustalić przyczynę, wpisujemy ręcznie (lub za pomocą Ctrl+C + Ctrl+V) podany układ równań:

ContourPlot[{7 x-4 y==3,2 x+2 y==-3},{x,-12,12},{y,-12,12}]

Aby narysować wykres, musimy użyć polecenia Evaluate, które zamieni napis uklad, na równania zapisane pod taką nazwą:

ContourPlot[Evaluate[uklad],{x,-12,12},{y,-12,12}]

Aby dokładniej odczytać rozwiązanie, zmniejszamy skalę na osiach:

ContourPlot[Evaluate[uklad],{x,-2,2},{y,-2,2}]

Aby z dobrą dokładnością odczytać współrzędne wybieramy z  Menu -> Graphics -> Drawing Tools, i używając ,,celownika’’ kopiujemy współrzędne klawiszami Ctrl+C. Następnie wklejamy poniżej za pomocą Ctrl+V:

rozw1={-0.2716,-1.23}

Spodziewany wynik to pewne liczby wymierne. Próby ich identyfikacji na podstawie rozwinięcia dziesiętnego można dokonać poleceniem Rationalize:

Rationalize[rozw1]

Niestety, nie jest to prawidłowe rozwiązanie:

uklad/.{x->%[[1]],y->%[[2]]}

Spróbujmy zastosować FindRoot aby znaleźć lepsze przybilżenie dziesiętne:

FindRoot[uklad,{{x,rozw1[[1]]},{y,rozw1[[2]]}}]

{x,y}/.%

Rationalize[%]

Tym razem Rationalize poprawnie zidentyfikowało liczby wymierne, a uzyskane w taki sposób rozwiązanie symboliczne jest poprawne:

uklad/.{x->%[[1]],y->%[[2]]}

Powyższa procedura jest znana pod nazwą regresji symbolicznej.

Przykłady powyżej nadużywały znaku % oraz klawiszy Ctrl+C Ctrl+V. Zróbmy to samo za pomocą komend. Na początek zapisujemy odczytane z wykresu liczby:

{xSTART,ySTART}=rozw1

xSTART

ySTART

Teraz wywołujemy FindRoot:

rozw2={x,y}/.FindRoot[uklad,{{x,xSTART},{y,ySTART}}]

Także niewiadome można zapisać do pewnej zmiennej:

niewiadome={x,y}

I rozwiązać układ za pomocą np: NSolve, która to instrukcja nie potrzebuje wartości startowych:

NSolve[uklad,niewiadome]

Aby dostać same liczby, musimy wybrać pierwszy element z i listy symbolem [[1]] lub ewentualnie poleniem First:

rozw3={x,y}/.NSolve[uklad,niewiadome][[1]]

rozw3={x,y}/.First[NSolve[uklad,niewiadome]]

Rozwiązania numeryczne uzyskane za pomocą FindRoot oraz NSolve są w tym prostym przypadku identyczne:

rozw2-rozw3

Metody symboliczne

Przechodzimy do rozwiązania symbolicznego. Jest to bardzo proste, bo układ jest linowy:

Solve[uklad,{x,y}]

Lub, przy użyciu bardziej eleganckiej notacji:

Solve[uklad,niewiadome]

rozw4=niewiadome /. %[[1]]

Reduce daje taki sam wynik, ale w innej formie:

wynikReduce=Reduce[uklad,{x,y}]

Aby skonwertować taką postać to listy reguł transformacyjnych, możemy użyć Solve:

Solve[wynikReduce,{x,y}]

lub za pomocą dedykowanego polecenia ToRules:

ToRules[wynikReduce]

Dla porządku sprawdzamy jeszcze co zrobi FindInstance:

FindInstance[uklad,{x,y}]

Metody macierzowe

Większość studentów zapewne słyszała metodach macierzowych rozwiązywania układów równań liniowych (metodzie wyznaczników). Mathematica dysponuje poleceniami pozwalającymi na rozwiązywanie w postaci macierzowej. Oto nasz układ równań, w postaci tablicy:

uklad//TableForm

Z liczb stojących prze x i y tworzymy macierz główną m:

m={{7,-4},{2,2}}

MatrixForm[m]

Poniższy sposób pochodzi z wykładu prof. Longi i pozwala na równoczesne poprawne utworzenie macierzy i wyświetlenie jej w ładnej postaci:

m={{7,-4},{2,2}};
MatrixForm[m]

Analogicznie, tworzymy macierz niewiadomych X oraz macierz wyrazów wolnych b:

X={x,y}

b={3,-3}

Równanie ma teraz poniższą postać, gdzie symbol kropki dziesiętnej oznacza mnożenie macierzowe:

m.X==b

Taki układ mozna rozwiązać używając np: Solve:

Solve[m.X==b,X]

Jak dotąd poza zmianą notacji na trudniejszą nic nie uzyskaliśmy. W szczególności, system nie jest na tyle ,,mądry’’, aby odgadnąć czym jest niewiadoma X:

Clear[X]

Solve[m.X==b,X]

Aby rozwiązać system w postaci macierzowej należy użyć LinearSolve:

LinearSolve[{{7,-4},{2,2}},{3,-3}]

Ewentualnie, możemy obliczyć macierz odwrotną do m

macierzOdwrotna=Inverse[m]

i pomożyć ją przez macierz wyrazów wolnych b:

macierzOdwrotna.b

Możemy użyć także szkolnej metody, obliczając wyznaczniki:

w=Det[m]

wx=Det[{{3,-4},{-3,2}}]

wy=Det[{{7,3},{2,-3}}]

{wx/w,wy/w}

Uwaga końcowa: nie wyciągajmy armaty aby zabić komara. Metody macierzowe mają znaczenie w dwóch przypadkach: (1) system liniowy jest bardzo duży, dziesiątki i setki niewiadomych (2) analiza teoretyczna, w której postać macierzowa ułatwia zrozumienie problemu. W razie wątpliwości zastosuj Solve!

Zamiana na zagadnienie minimalizacji

Każdy układ równań można zamienić na jedno równania, stosując trik Feynmana: w każdym równaniu przenosimy wszystko na jedną stroną, podnosimy do kwadratu (ewentualnie wartość bezwzględną lub inną normę) a nastepnie dodajemy. Tak utworzone wyrażenie ma dwie własności:
(1) jest równe zeru, gdy każdy ze składników jest równy zeru
(2) jego minimum to zero, w punkcie, gdzie zachodzi każde z równań

Rozwiązanie układu równań jest więc równoważne rozwiązaniu równania lub znalezieniu minimum wyrażenia:

w=(7x-4y-3)^2+(2x+2y+3)^2

Minimum mozna znaleźć graficznie :

Plot3D[w,{x,-2,2},{y,-2,2}]

ContourPlot[w,{x,-1,1},{y,-2,0},Contours->64]

rozw5={-0.2637,-1.256}

rozw5-rozw1

Minimum szukamy numerycznie:

FindMinimum[w,{x,y}]

NMinimize[w,{x,y}]

Lub symbolicznie:

Minimize[w,{x,y}]

Problem jest też równoważny rozwiązaniu równania w=0

FindInstance[w==0,{x,y},Reals]

Tradycyjne polecenia Solve i Reduce, nie radzą sobie jednak z tego typu problemem (jedno równania, 2 niewiadome):

Solve[w==0,{x,y}]

Reduce[w==0,{x,y}]

Warto zapamiętać powyższy trik, gdyż pozwala na rozwiązanie także wyjątkowo dziwacznych równań. Zobacz np: ostatni przykład w istrukcji do NMinimize:

?NMinimize

Przykład 2: Rozwiąż system wielomianowy

2 (x^2)+2 (y^2)+2 (y^3)+x (y+(y^3))+(x^3) (-1-(y^2)+(y^3))=1
2 y+2 (y^3)-(x^3) y (1+y)+x y (-1+2 y+(y^2))+(x^2) (1+2 (y^2)-(y^3))=1

Metody graficzne i numeryczne: NSolve

Układy równań zawierające tylko wielomiany stanowią specjalną klasę w zasadzie zawsze rozwiązywalnych problemów. Szczegółowo opisano to w tutorialu RealPolynomialSystems. Wyniki mają postać na ogól niezrozumiałą dla studentów, ale w rzeczywistości są po prostu liczbami. Przykład ma na celu oswojenie z tego typu problemami.

rownanie1=2 (x^2)+2 (y^2)+2 (y^3)+x (y+(y^3))+(x^3) (-1-(y^2)+(y^3))==1;

rownanie2=2 y+2 (y^3)-(x^3) y (1+y)+x y (-1+2 y+(y^2))+(x^2) (1+2 (y^2)-(y^3))==1;

ukladWielomianowy = {rownanie1,rownanie2}

Tradycjnie, zaczynamy od wykresów:

ContourPlot[Evaluate[ukladWielomianowy],{x,-12,12},{y,-12,12}]

W rozważanym przykładzie, ustalenie który z wykresów odpowiada pierwszemu równaniu nie jest możliwe. Dodajmy więc jawnie kolory linii:

ContourPlot[Evaluate[ukladWielomianowy],{x,-4,4},{y,-12,12},ContourStyle->{{Red,Thick},{Black,Thick}}]

Przykład jest dosyć trudny, na pojedynczym wykresie nie mogę doliczyć się nawet ilości rowiązań. Musimy wykonać kilka wykresów:

SetOptions[ContourPlot, ContourStyle->{{Red,Thick},{Black,Thick}}];

ContourPlot[Evaluate[ukladWielomianowy],{x,-3/2,-1/2},{y,-24,24}]

ContourPlot[Evaluate[ukladWielomianowy],{x,-1,3},{y,-2,2}]

{{-1.113,6.56},{-0.575,0.4736}}

{{-0.6193,0.1963},{0.6314,0.3093},{1.003,-0.02966},{2.011,0.5029}}

Wychodzi, że rozwiązań jest 6, ale mam wątpliwości to do tych z x≈-0.5:

ContourPlot[Evaluate[ukladWielomianowy],{x,-1,0},{y,-1,1}]

{{-0.622,0.2178}}

Ostatecznie mamy więc 5 rozwiązań:

rozw0={{-1.113,6.56},{-0.622,0.2178},{0.6314,0.3093},{1.003,-0.02966},{2.011,0.5029}}

Spróbujemy nieco zautomatyzować proces szukania rozwiązań numerycznych. Najpierw utworzymy polecenie jednolinijkowe (ang: one-liner), działające dla pewnej wartości, np: 1

{{x,y},rozw0[[1]]}

Transpose[%]

FindRoot[ukladWielomianowy,Transpose[{{x,y},rozw0[[1]]}]]

Teraz zamieniamy liczbę 1 na symbol k, a następnie wstawiamy całość to Table:

Table[FindRoot[ukladWielomianowy,Transpose[{{x,y},rozw0[[k]]}]],{k,1,5}]

rozw1={x,y}/.%

Trzy pierwsze rozwiązania są pewnymi liczbami rzeczywistymi (zmiennoprzecinkowymi), natomiast dwa ostatnie wyglądają
na {1,0}  - bo liczba rzędu (10^(-16))(  $MachineEpsilon ) jest prawdopodobnie zerem, oraz {2,1/2} . Sprawdźmy to:

ukladWielomianowy/.{x->1,y->0}

ukladWielomianowy/.{x->2,y->1/2}

Poleceniem numerycznym przeznaczonym do rozwiązywania układów wielomianowych jest NSolve. Nie wymaga podawania punktów startowych, i znajduje wszystkie rozwiązania!

NSolve[ukladWielomianowy,{x,y}]

rozw2={x,y}/.%

Otrzymaliśmy wszyskie rozwiązania, w tym zespolone. Na przyszlość przyda się sposób wybrania przypadków (ang: cases) w których obie liczby są rzeczywiste. Polecenie Cases wybiera z listy elementy pasujące do danego wzorca (ang. pattern). Poniżej działający wzorzec wybierający pary liczb rzeczywistych:

rozw3=Cases[rozw2,{_Real,_Real}]

Chcielibyśmy porównań wyniki FindRoot oraz NSolve, ale są w różnej kolejności. Najprościej jest je posotrtować:

Sort[rozw1]

Sort[rozw3]

%-%%

Wyniki zgadzają się przynajmniej do 11 cyfr znaczących (po przecinku). To 100 tysięcy razy mniej niż dokładność procesora, ale ciągle dobra zgodność. Gdyby elementy które chcemy porównań nie były możliwe do posortowania, lub chcieibyśmy zachować oryginalną kolejność, możemy wygenerować macierz:

Table[Norm[rozw1[[i]]-rozw3[[j]]],{i,1,5},{j,1,5}]//MatrixForm//Chop

Wyżej widzimy jak należy przyporządkować rozwiązania: 1->1, 2->5, 3->4, 4->2, 5->3. Układ możemy od razu rozwiązywać w domenie liczb rzeczywistych:

rozw4 ={x,y}/.NSolve[ukladWielomianowy,{x,y},Reals]

Sort[rozw4]-Sort[rozw3]

Wyniki powyżej są nieco zaskakujące. Różnice pomiędzy dwoma wywołaniami NSolve są spore, znacznie większe od $MachineEpsilon. Na koniec podaję, jak znaleźć rozwiązania z dowolną precyzją (dla CPU jest to 16):

rozw5 ={x,y}/.NSolve[ukladWielomianowy,{x,y},Reals,WorkingPrecision->64]

Dygresja: pierwiastki wielomianów o współczynnikach całkowitych

Aby studenci nie doznali szoku na widok rozwiązań symbolicznych systemu wielomianowego, pokażę kilka prostych przykładów na pierwiastki wielomianów. Zaczniemy od dobrze wszystkim znanego równania kwadratowego, które rozwiązujemy na kartce ,,deltą’’ (Δ):

Solve[x^2+2x-3==0,x]

Znalezienie rozwiązań jest równoważne rozłożeniu wielomianu na czynniki (faktoryzacja), do czego służy polecenie Factor:

Factor[x^2+2x-3]

Zastosujmy to samo rozumowanie do równania kwadratowego z niewymiernymi pierwiastkami:

Solve[x^2-3x-3==0,x]

Factor[x^2-3x-3]

Wielomian nie został rozłożony na czynniki. Wynika to z tego, że wielomian po rozłożeniu musiałby zawierać liczby niecałkowite, a konkretnie liczbę (Sqrt[Δ])=(Sqrt[21]).

Aby rozłożyć taki wielomian na czynniki nie wystarczy znajomość liczb całkowitych. Musimy znać pierwiastki, a w tym konkretnym przypadku musimy wprowadzić dodatkowy symbol, pierwiastek z delty:

Factor[x^2-3x-3,Extension->Sqrt[21]]

Tego typu wyrażenia nazywamy liczbami algebraicznymi. Są one utworzone z liczb wymiernych, dodawania, mnożenia i rozwiązania pewnego wielomianu (w rozważanym przykładzie jest to x^2-21==0 ).

x/.Solve[x^2-3x-3==0,x]

ToNumberField[%]

Kilka podobnych przykladów:

Factor[x^3-2]

Factor[x^3-2,Extension->2^(1/3)]

TraditionalForm[%]

Factor[x^4+1]

Factor[x^4+1,Extension->Sqrt[2]]

Factor[x^4+1,Extension->I]

Lekcja z powyżsych przykładów jest taka, że aby zapisać rozwiązania, musimy wprowadzić pewne symbole, które nie występują w ,,treści’’ zadania. Np: rozwiązanie równania kwadratowego wymaga wprowadzenia symbolu pierwiastka. Co dzieje się w przypadku równań wyższych stopni? (Zobacz też: tutorial/PolynomialsOverAlgebraicNumberFields ). Jeżeli pierwiastki są całkowite, sprawa jest prosta:

w3a=x^3+x^2-10x+8

Solve[w3a==0,x]

Factor[w3a]

Zagadnienie w ogólnym przypadku komplikuje się:

w3b=x^3-6x+2

Solve[w3b==0,x]

Mamy trzy rozwiązania, wyglądające na zespolone. Faktoryzacja jest możliwa po wprowadzeniu pewnego pierwiastka trzeciego stopnia:

Factor[w3b, Extension->{(-1+I Sqrt[7])^(1/3)}]

Spróbujmy jawnie zażądać rozwiązania rzeczywistego:

Solve[w3b==0,x,Reals]

Reduce[w3b==0,x]

Reduce[w3b==0,x,Cubics->True]

Cokolwiek byśmy nie robili, otrzymujemy wyrażenia zawierające jednostkę urojoną, lub obiekt Root. Naszym celem jest oswojenie się z nim, i zrozumienie czym jest. Jest to ważne, bo będzie się on pojawiał w waszych wynikach bardzo często. Jako wstęp rozwiążemy równanie 3 stopnia numerycznie:

NSolve[w3b,x]

Każdy widzi, że wynikiem są trzy liczby rzeczywiste. Logika podpowiada, że Solve także daje 3 liczby rzeczywiste:

Solve[w3b==0,x,Reals]

N[%]

Root jest więc zapisem liczby będacej n-tym pierwiastkiem wielomianu. Wewnątrz wielomian zapisany jest za pomocą funkcji szkieletowej (ang: pure function), w której rolę zmiennej x pełni znak # (hash). Drugi argument stanowi numer kolejny pierwiastka.

w3b

r=Root[#^3-6#+1&,1]

Konwersji do postaci ,,pierwiatkowej’’ dokonujemy poleceniem ToRadicals (ang. radical - pierwiastek). W ogólności, dla stopnia wielomianu 5 lub więcej, taka konwersja jest niemożliwa.

ToRadicals[r]

ComplexExpand[Re[%]]

FullSimplify[%]

Pomimo, że same pierwiastki muszą być wyrażone za pomocą symbolu Root , wynik operacji na nich może okazać się prosty. Np. znane wzory Viete’a:

Root[#^3-6#+1&,1]+Root[#^3-6#+1&,2]+Root[#^3-6#+1&,3]

Simplify[%]

Root[#^3-6#+1&,1]*Root[#^3-6#+1&,2]*Root[#^3-6#+1&,3]

Simplify[%]

Dla treningu, zapiszmy za pomocą funkcji Root rozwiązanie równania kwadratowego:

Root[#^2-3#-3&,1] (* pierwszy pierwiastek *)

Root[#^2-3#-3&,2] (* drugi pieriwastek *)

Root[#^2-3#-3&,3] (* trzeciego nie ma *)

Ilustracją skuteczności symbolu Root, niech będzie analiza równania 5 stopnia:

w5=x^5+x+1

Solve[w5==0,x,Reals]

Factor[w5,Extension->Root[1-((#1)^2)+((#1)^3)&,1]]

%/.Root[1-((#1)^2)+((#1)^3)&,1]->x0

w5/(x^2+x+1)

Simplify[%]

ToRadicals[Root[1-((#1)^2)+((#1)^3)&,1]]

Zwracam także uwagę, że istnieje rozszerzenie polecenia Root na równania przestępne.

FindInstance[x^x==1-x,x,Reals]

Taka postać na ogół nie jest możliwa do wprowadzenia z klawiatury:

Root[{-1+#1+((#1)^(#1))&,0.30365912703}]

Końcowa uwaga: symbol Root jest po prostu oznaczeniem liczby. Nie powinien budzić większego zdziwienia, niż symbol (Sqrt[2]).

Rozwiązanie symboliczne systemu wielomianowego

Wracamy do systemu wielomianowego.

ukladWielomianowy

Jego rozwiązanie nie mieści się na ekranie:

Solve[ukladWielomianowy,{x,y}]

Aby w ogóle pojąć co jest wynikiem, obliczny jego wartość numeryczną:

N[%]

Dimensions[%]

Powyżej widzimy, że mamy 17 rozwiązań, wiekszość zespolonych. Dodajemy opcję Reals. Nauczka dana poprzednio skłania do zakończenia polecenia średnikiem. Oglądanie tego typu rezultatu i tak nic nie da.

Solve[ukladWielomianowy,{x,y},Reals];

Poleceniem Dimensions sprawdzamy ile jest rozwiązań (podaje ono wymiary listy):

Dimensions[%]

Mamy 5 rozwiązań, tak jak wynika z wykresu, oraz z obl. numerycznych.

rozw6={x,y}/.%%;

rozw6//Dimensions

rozw6[[1]]

rozw6[[2]]

rozw6[[3,1]]

Sprawdzamy dokładność numeryczną. Zgadza się perfekcyjnie aż do zadanych 32 cyfr po przecinku.

Sort[N[rozw6,32]]-Sort[N[rozw5,32]]

Gigantyczny output w pierszym momencie odpycha, ale to tylko kwestia notacji.

rozw6[[3,2]]

Sprowadźmy powyższe wyrażenie do ,,ludzkiej’’ postaci:

%/.Root[-24+32 #1+76 ((#1)^2)-48 ((#1)^3)-100 ((#1)^4)-131 ((#1)^5)+170 ((#1)^6)+116 ((#1)^7)-4 ((#1)^8)-15 ((#1)^9)-25 ((#1)^10)-10 ((#1)^11)-25 ((#1)^12)+25 ((#1)^13)-6 ((#1)^14)+3 ((#1)^15)&,1]->r

Simplify[%]

Ostatnie wyrażenie jest już na tyle krótkie, że można by pokusić się nawet o przepisanie go na kartkę, oczywiście w razie potrzeby.