Celem zadania jest numeryczne rozwiązanie układu dwóch równań z dwiema niewiadomymi. Poniżej, jako przykład, będziemy rozwiązywać następujący układ:

x+ 6 y - 2 == 0
exp(-(x^2)+y)-1==0

Dla wygody zapisujemy rozwiązywany układ w postaci jednej funkcji przyjmujacej jeden argument. Funkcja przyjmuje jako wartość dwuelementowy wektor (listę), i zwraca dwuelementowy wektor. Dzięki temu napisany program nie będzie różnił się zbytnio od programu rozwiązującego równania z jedną niewiadomą.

fun[{x_,y_}]:={x+6y-2,Exp[-x^2+y]-1}

fun[{1,2}]

fun[{1.0,2.0}]

Dla wyrobienia poglądu na temat badanego układu rysujemy wykresy zadanych funkcji ( y(x) w postaci uwikłanej f(x,y)=0 ):

wykres1=ContourPlot[fun[{x,y}][[1]]==0,{x,-5,5},{y,-5,5},ContourStyle->Red]

wykres2=ContourPlot[fun[{x,y}][[2]]==0,{x,-5,5},{y,-5,5},ContourStyle->Blue]

Po nałożeniu wykresy mogą (lub nie) przecinać się. Punkty przecięcia są rozwiązaniami graficznymi układu równań:

Show[wykres1,wykres2]

Show[wykres1,wykres2,PlotRange->{{-1,1},{0,1}}](* powiększenie interesującego obszaru *)

Aby rozwiązać układ numerycznie, musimy wybrać wartośc startową w pobliżu rozwiązania.

poprzednia={0.5,0.5}

Obliczenie kolejnego przybliżenia wymaga znajomości jakobianu, czyli macierzy 2x2 zawierającej pierwsze pochodne funkcji:

jakobian = {{D[fun[{x,y}][[1]],x],D[fun[{x,y}][[1]],y]},{D[fun[{x,y}][[2]],x],D[fun[{x,y}][[2]],x]}}

Obliczamy numeryczną wartość jakobianu w wybranym punkcie:

jakobianNum=jakobian/.{x->poprzednia[[1]],y->poprzednia[[2]]}

Nowe, dokładniejsze przybliżenie otrzymujemy odejmując od poprzedniej wartości macierz odwrotną do jakobianu (Inverse) pomnożoną przez wektor zawierający wartości funkcji:

nowa=poprzednia-Inverse[jakobianNum].fun[poprzednia]

poprzednia=nowa (* czynnośći powtarzamy do skutku, wracając 2 linijki wyżej *)

Automatycznie można to zrobić np. funkcją Nest. W tym celu tworzymy jednoargumentową funkcję obliczającą macierz odwrotną do jakobianu:

odwrotnosc[{x_,y_}]:=Evaluate[Inverse[jakobian]]

odwrotnosc[{0.5,0.5}]

Wykorzystując funkcję Nest możemy zautomatyzować iteracje na kilka sposobów podanych poniżej:

Nest[(N[#-odwrotnosc@#.fun@#])&,{0.5,0.5},16](* wykonaj podaną funkcję 16 razy, zaczynając od {0.5,0.5} *)

NestList[(N[#-odwrotnosc@#.fun@#])&,{0.5,0.5},16] (* jak wyżej, ale zapisuj wyniki pośrednie do listy *)

NestWhile[(N[#-odwrotnosc@#.fun@#])&,{0.5,0.5},Norm[fun@#]>0.001&] (* wykonuj dopóki jest spełniony warunek *)

punkty=NestList[(N[#-odwrotnosc@#.fun@#])&,{-1,-1},12]

Manipulate[Show[wykres1,wykres2,ListPlot[punkty[[1;;num]]],PlotRange->{{-2,2},{-2,2}}],{num,1,12,1}]