Rozwiązywanie równań z jedną niewiadomą: 
przypadek z parametrem

Równanie z parametrem.

Przykład 2: Rozwiąż w zależności od rzeczywistego parametru p równanie:

(Sqrt[x-p])=x+p

Rysowanie wyrażeń z parametrami

Zgodnie z podaną wcześniej procedurą powinniśmy rozpocząć od narysowania wykresów lewej i prawej strony równania :

Plot[{Sqrt[x-p],x+p},{x,-4,4}]

Powyższy błąd jest zadziwiająco częsty. Program komputerowy potrafi narysować jedynie wykres funkcji, za pomocą której może obliczyć liczby rzeczywiste, przeskalowane następnie do współrzędnych na ekranie komputera. Oznacza to, że rysować można jedynie numeryczne wartości, wyłącznie rzeczywiste! Który pixel na ekranie LCD miałby zostać podświetlony, jeżeli jego współrzędne to, powiedzmy, x=Sqrt[-1],  y=p ?
Aby narysować wykres, instrukcji Plot musimy podać wyrażenie, które posiada rzeczywistę wartość numeryczną. Powyższy, błędny przypadek, nie spełnia tego warunku. Po wstawnieniu dowolnej wartości x z przedziału {-4,4} otrzymamy:

{Sqrt[x-p],x+p}/.{x->2}

Dlatego, aby narysować wykres musimy nadać parametrowi p pewną wartość liczbową, powiedzmy p=2:

p=2

Teraz powtarzamy instrukcję Plot :

Plot[{Sqrt[x-p],x+p},{x,-4,4}]

Inna wartość p:

p=-1

Plot[{Sqrt[x-p],x+p},{x,-4,4}]

Formalnie, możemy analizę opisaną w przypadku bez parametrów powtorzyć dla dowolnej wartości parametru p. O ile na pierwszy rzut oka taka procedura może wydawać się niesłychanie uciążliwa, dla osoby która potrafi programować w Mathematice nie ma różnicy czy analizujemy dwa, 100, czy milion przypadków. Robi to za nas komputer.

Manipulate

Na początek poznamy wygodne polecenie Manipulate, stwotrzone z myślą o analizowaniu rozmaitych zagadnień w których występują zmienne parametry. Zaczniemy od obliczenia Pi z zadaną dokładnością:

Manipulate[N[Pi,n+1],{n,1,100,1}]

Wzory skróconego mnożenia:

Manipulate[Expand[(a+b)^n],{{n,2},0,12,1}]

,,Ładniejsza’’ wersja:

Manipulate[Superscript["(a+b)",n]==Expand[(a+b)^n],{{n,2},0,12,1}]

To samo ze wzoru na dwumian Newtona :

Manipulate[Sum[Binomial[n,k] a^k b^(n-k),{k,0,n}],{{n,2},0,12,1}]

Binomial to symbol Newtona. Aby wyświetlić go w postaci tradycyjnej zaznaczamy i naciskamy Ctrl+Shift+T

Binomial[n,k]

Taki sam efekt daje polecenie TraditionalForm:

TraditionalForm[Binomial[n,k]]

Polecenia zakończone frazą Form, np: FullForm, TreeForm, TraditionalForm, TableForm mają na celu wyświetlanie na ekranie wyrażeń w dogodnej formie. 

UWAGA: po zastosowaniu jednego z powyższych poleceń nie jest możliwe wykonywanie poprawnych dalszych obliczeń. Służą one tylko do prezentacji graficznej!

Wracamy do równania. Zastosujmy Manipulate do narysowania wykresu prawej i lewej strony interesującego nas równania :

Clear[p](* Usuwam wcześniej przypisane wartości p *)

Manipulate[
Plot[{Sqrt[x-p],x+p},{x,-4,4}]
,{p,-1,1}]

Posługując się suwakiem łatwo ustalamy, że nasze równanie nie ma rozwiązań dla p>0.12, dwa rozwiązania (obydwa w przedziale od 0 do 2) dla 0<=p<0.12 i jedno rozwiązanie dla p<0.  Wiedza ta jest cenna, i pozwoli nam uniknąć katastrofalnych błędów.

Rozwiązanie numeryczne równania z parametrem

Rozważmy np: przypadek p=0.2 o którym wiemy, że nie ma żadnego rozwiązania i zastosujmy FindRoot:

p=0.2

FindRoot[Sqrt[x-p]==x+p,{x,1/2}] (* Dla tej wartości startowej otrzymujemy błędny wynik ! *)

Sqrt[x-p]==x+p/.%  (* UWAGA! Błędny wynik!!! *)

FindRoot[Sqrt[x-p]==x+p,{x,0}] (* Dla innej wartości startowej otrzymujemy dobry wynik, ale zespolony! *)

Sqrt[x-p]==x+p/.% (* Sprawdzenie rozwiązania. *)

Teraz weźmy p = 0.1, dla którego z wykresu i Manipulate spodziewamy się dwóch rozwiązań :

p=0.1

FindRoot[Sqrt[x-p]==x+p,{x,0}]

FindRoot[Sqrt[x-p]==x+p,{x,1}]

FindRoot[Sqrt[x-p]==x+p,{x,1/3}]

Należy się parę słów komentarza do tych wyników. W pierwszym przypadku wynik jest zespolony, ale część urojona jest bardzo mała. Część rzeczywista odpowiada spodziewanemu rozwiązaniu. Jest to wynikiem rozpoczęcia iteracji od zespolonej wartości startowej:

Sqrt[x-p]/.x->0

Niezależnie od tego, uzyskany wynik powinien być rzeczywisty (zespolony z częścią urojoną zero). Z powodu skończonej precyzji obliczeń wykonywanych przez procesor, część urojona jest zwykle bardzo mała, ale różna od zera. Orientacyjnie, wynik rzędu $MachineEpsilon lub mniejszy (co do wartości bezwzględnej):

$MachineEpsilon

jest prawdopodobnie równy zeru. Do obcinania liczb będących w większości przypadków skutkiem błędów zaokrąglania używamy polecenie Chop.

FindRoot[Sqrt[x-p]==x+p,{x,0}]

x/.%

Chop[%]

UWAGA: punktów które posiadają tego typu artefakty w części urojonej, nie da się narysować instrukcją Plot ! Należy pozbyć się ich poleceniem Chop (obetnij), Re (weź część rzeczywistą) lub Abs (oblicz wartość bezwzględną).

W przypadku, gdy istnieją dwa bliskie rozwiązania, FindRoot zwykle losowo daje jedno lub drugie, w zależności od warunków początkowych. Co gorsza, ze względu na skończoną dokładność obliczeń, nie są one zwykle identyczne. Tego typu problem można elegancko rozwiązać metodami bioinformatycznymi, np. poleceniem FindClusters (więcej na dalszych wykładach).

Dla p<0 spodziewamy się jednego rozwiązania i powyższy problem nie wystepuje. Można więc  zautomatyzować proces obliczeń numerycznych.
Na początek ,,rozgrzewamy się’’ przypadkiem p=-0.5.

Clear[p]

p=-0.5

FindRoot[Sqrt[x-p]==x+p,{x,2}]

x/.FindRoot[Sqrt[x-p]==x+p,{x,2}] (* W ten sposób jako wynik dostajemy samą liczbę *)

Teraz jesteśmy przygotowani do wygenerowania rozwiązań dla dużej liczby wartości p. Najwygodniejszym poleceniem jest tu Table. Składnia jest podobna do Manipulate. Table w odróżnieniu generuje wszystkie możliwe wartości p.

Manipulate[p,{p,-10,0,2}]

Table[p,{p,-10,0,2}]

Teraz możemy umieścić FindRoot wewnątrz Manipulate:

Manipulate[FindRoot[Sqrt[x-p]==x+p,{x,2}],{p,-10,10,2}]

W ten sposób Manipulate generuje same liczby będące rozwiązaniem:

Manipulate[x/.FindRoot[Sqrt[x-p]==x+p,{x,2}],{p,-10,0,2}]

Teraz zamieniamy po prostu Manipulate na Table:

Table[x/.FindRoot[Sqrt[x-p]==x+p,{x,2}],{p,-10,0,2}]

Otrzymane rozwiązania możemy nanieść na wykres (bardzo niedoskonały):

ListPlot[%]

Punkty wydają się układać na linii prostej. Spróbujmy zagęścić wykres:

Table[x/.FindRoot[Sqrt[x-p]==x+p,{x,2}],{p,-10,0,0.1}]

ListPlot[%]

A jednak nie... Końcówka wykresu ,,zagina się’’. Aby uniknąć wyświetlania na ekranie dużych wyników, polecenie kończymy średnikiem:

Table[x/.FindRoot[Sqrt[x-p]==x+p,{x,2}],{p,-10,0,0.01}];

ListPlot[%]

Na niedoskonałych wykresach wcześniej, na osi poziomej był numer kolejny obliczanych punktów. Chcielibyśmy, aby była to wartość p. Uzyskujemy to generując w Table pary {p, x}. Pierwsza z liczb będzie na osi X, druga na osi Y:

Table[{p,x/.FindRoot[Sqrt[x-p]==x+p,{x,2}]},{p,-10,0,0.2}]

ListPlot[%]

Do póżniejszego użytku warto zapisać wyniki obliczeń pod zmienną o wygodnej dla nas (oraz współpracowników!) nazwą, np: wyniki:

wyniki=Table[{p,x/.FindRoot[Sqrt[x-p]==x+p,{x,2}]},{p,-10,0,0.2}];

Także wykres jest pewnym wyrażeniem w Mathematice, i możemy go zapisać pod pewną zmienną:

wykres=ListPlot[wyniki]

Temat powyższy będzie szczegółowo omawiany na wykładzie dotyczącym analizy danych numerycznych i doświadczalnych. Na razie potrzebny jest nam tylko wykres do porównania wyników numerycznych i symbolicznych.

Podsumowanie  zdobytych graficznie i numerycznie informacji

a) dla p<0 równanie posiada jedno rozwiązanie, które bez problemu potrafimy obliczyć numerycznie; jego wykres w zależności od p jest prawie liniowy, za wyjątkiem p blisko zera

b) w przedziale 0<=p<=p* mamy dwa rozwiązania; można je obliczać numerycznie, ale są pewne problemy z ich rozróżnieniem i urojonymi artefaktami; p* ≈ 0.12; rozwiązania są w przedziale 0<x<2

c) dla p>p* nie ma rozwiązań rzeczywistych, FindRoot podaje błędne lub zespolone wyniki

d) dla p=p* powinno być jedno rozwiązanie

Dopiero teraz możemy zabrać się za metody symboliczne!

Metody symboliczne

Clear[p]

Zaczynamy od Solve:

rozwiazanie0=Solve[Sqrt[x-p]==x+p,x]

Zgodnie z oczekiwaniami, Solve podało naiwny wynik, uzyskany poprzez podniesienie obustronne do kwadratu i rozwiązanie  równania kwadratowego ,,deltą’’. Na szczęście, przeprowadziliśmy wcześniej analizę graficzną i numeryczną, dzięki której wiemy, że nie może to być całkowicie poprawne rozwiązanie! Sprawdźmy np: dla p=-2:

x/.rozwiazanie0

%/.p->-2

Jako początkujący użytkownicy, nie wnikamy w szczegóły składniowe i kopiujemy najpierw pierwszy (z plusem) , a potem drugi wynik (z minusem) znaną kombinacją klawiszy Ctrl+C Ctrl+V (metoda Kopiego-Pasta :-)

Sqrt[x-p]==x+p/.{p->-2,x->(1/2) (5-(Sqrt[17]))}

Simplify[%]

FullSimplify[%] (* To wyrażenie to tzw. zagnieżdżone pierwiastki (nested radicals) - tego typu operacje bywają bardzo trudne - jak nie wierzysz to spróbuj uprościć (nie znając lewej strony) używając kartki i ołówka *)

Jak widać, dla p<0,  jedno z rozwiązań jest błędne. Spróbujmy explicite zażądać rozwiązania rzeczywistego:

rozwiazanie1 = Solve[Sqrt[x-p]==x+p,x,Reals]

Powyższa odpowiedź jest znacznie bardziej precyzyjna. Mamy dwa rozwiązania, ale jedno z nich jest zdefiniowane warunkowo ( ConditionalExpression - wyrażenie warunkowe ) tylko w przedziale 0<p<1/8. Drugie podobnie, określono dla 0<p<1/8 oraz dla p<0. Wszystko się więc zgadza z wcześniejszą analizą. Pozostaje mimo to mały niedosyt: nie ma przypadku p=0 a także p=1/8 :

rozwiazanie1/.{p->1/8}

Użyjmy więc ,,mocniejszego’’ polecenia Reduce:

Reduce[Sqrt[x-p]==x+p,x]

Wynik jest podejrzany, za sprawą komunikatów o błędach. Pojawia się określenie ,,branch cuts’’, czyli tzw. ,,cięcia’’, które dotyczą zachowania się wielowartościowych funkcji zespolonych. Nas interesują tyko rozwiązania rzeczywiste:

rozwiazanie2=Reduce[Sqrt[x-p]==x+p,x,Reals]

Teraz odpowiedź jest kompletna, aczkolwiek trudna do odczytania. Reduce podaje wynik w postaci zdań logicznych. Symbol && to logiczne ,,i’’ (And), natomiast symbol II to logiczne ,,lub’’ (Or).

rozwiazanie2//TraditionalForm

Dla porządku wyjaśniam, że FindInstance nie rozwiązuje równań z parametrami.

FindInstance[Sqrt[x-p]==x+p,x]

Może to zrobić dopiero po przypisaniu konkretnej wartości p:

p=1/8

FindInstance[Sqrt[x-p]==x+p,x]

p=0

FindInstance[Sqrt[x-p]==x+p,x,3]

Na zakończenie porównajmy wynik numeryczny i symboliczny dla p<0. Wykres numeryczny, uzyskany poprzednio:

wykres

I rozwiązanie symboliczne dla p<0:

wykres2=Plot[(1/2) (Sqrt[1-8 p])+(1/2) (1-2 p),{p,-10,0}]

Dwa wykresy wyświetlamy na ekranie równocześnie poleceniem Show:

Show[wykres,wykres2]

Jak widać wszytko zgadza się w tym przedziale. Spróbujmy jeszcze narysować kompletny wykres, Rysujemy rozwiązanie w przedziale 0<p<1/8 :

wykres3=Plot[{(1/2) (1-2 p)+(1/2) (Sqrt[1-8 p]),(1/2) (1-2 p)-(1/2) (Sqrt[1-8 p])},{p,0,1/8}]

Polecenie Show pozwala wyświetlić dwa wykresy jednocześnie:

Show[wykres2,wykres3]

Show[wykres2,wykres3,PlotRange->{{-1,1/4},{0,3}}]

Metoda funkcji uwikłanej / wykresu konturowego

Poniższe rozważania stanowią wstęp do rozwiązywanie graficznego układów równań. Dociekliwy student potrafi zauważyć, że równanie o jednej niewiadomej z parametrem, jest równoważne zadaniu pewnej funkcji x(p) z postaci uwikłanej. Takie funkcje rysujemy poleceniem ContourPlot:

ContourPlot[Sqrt[x-p]==x+p,{p,-1,1/4},{x,0,3}]

Aby wykres bardziej przypominał wcześniejszy, usuwamy ramkę (Frame)  i dodajemy osie (Axes) :

ContourPlot[Sqrt[x-p]==x+p,{p,-1,1/4},{x,0,3},Frame->False,Axes->True]

Jest to dokładnie taki sam wykres jak poprzednio!

Użycie rozkładu cylindrycznego

Algorytmem o fundamentalnym znaczeniu dla algebry symbolicznej jest algebraiczny rozkład cylindryczny (ang: cylindrical algebraic decomposition, CAD) zbiorów liczb rzeczywistych. Zbiory te mogą być zdefiniowane w rozmaity sposób, a algorytm ma na celu doprowadzenie ich do standardowej postaci sumy rozłącznych przedziałów. W każdym z nich zbiór musi być określony za pomocą jawnych formuł matematycznych. Typowe zagadnienia, w których studenci ręcznie dokonuja takiej procedury to rozwiązywanie nierówności oraz obliczanie całek oznaczonych w celu obilczania pól figur i objętości brył. Tego algorytmu można użyć bezpośrednio:

Clear[p]

CylindricalDecomposition[Sqrt[x-p]==x+p,{p,x}]

Wynik powyższy jest lepszy od tego uzyskanego za pomocą Reduce, gdyż wartości parametru p są uporządkowane, od najmniejszych do największych. Łatwiej jest więc prześledzić odpowiedź.