(* Treść Zad. 2 zapisujemy pod nazwą "rownanie". Proszę zwrócić uwagę na różnicę pomiędzy pojedynczym = (przypisanie równania E^(p/x)==Gamma[x] do nazwy "rownanie") oraz podwójnym == (równość logiczna True/False) *)

rownanie = E^(p/x)==Gamma[x]

(* Do obejrzenia pełnej formy równania, tak jak widzi je Kernel programu Mathematica, można użyć polecenia FullForm, wywołanego w notacji ,,postfix-owej'', czytaj "do zmiennej "rownanie" zastosuj funkcję "FullForm" *)

rownanie //FullForm

(* Jak widzimy, podwójny znak równości to Equal. Posiadając taką wiedzę możemy elegancko wydobyć lewą i prawą stronę równania. Po prostu zamieniamy operatorem "./" (kropka slash, ReplaceAll czyli zamień wszystko) nagłówek Equal na List (lista). *)

lewaprawa=rownanie /.Equal->List

(* Zapamiętujemy pod nazwą lhs lewą stronę równania (ang. left-hand-side ), a pod nazwą rhs prawą stronę (ang. right-hand-side). W tym celu użyto operatora %, czyli poprzedniego wyniku. Średnik na końcu blokuje wyświetlanie na ekranie, bo i tak wiemy co się pojawi. *)

{lhs,rhs}=%;

(* Teraz narysujemy obie strony równania. W tym celu musimy chwilowo określić numeryczną wartość parametru p, bez tego komputer nic nie narysuje ani nie obliczy. Wygodnie można to zrobić poleceniem Block. Oprócz tego określam kolory linii na Red i BLue, aby było jasne która linia za co odpowiada: czerwnona  - lewa strona, niebieska - prawastrona. Nas interesuje głównie punkt przecięcia. *)

Block[{p=4/3},Plot[{lhs,rhs},{x,0,4},PlotStyle->{Red,Blue}]]

(* Punkt przecięcia możemy odczytać klikając na wykresie prawym przyciskiem myszy i wybierając polecenie "GetCoordinates". Następnie klikamy na miejsce przecięcia, nasiskamy Ctr+C i w następnej komórce (ang. Cell) Ctrl+V *)

{{2.756890007541163,1.6415505207734196`}}

(* Powyższa procedura z zasadzie kończy zadanie dla p=4/3. To co robimy dalej, to różne poziomy automatyzacji powyższego procesu. NA początek spróbojemy narysować kilka wykresów dla różnych wartości p. NA początek tabela wykresów dla p zaczynając od p=1 co 1/2 kończąc na p=4. *)

Table[Plot[{lhs,rhs},{x,0,4},PlotStyle->{Red,Blue}],{p,1,4,1/2}]

(* Aby wyświetlić powyższe wykresy równocześnie można użyć np: Show: *)

Show[%]

(* Analogicznie jak poprzednio, możemy przeklikać myszą wszystkie punkty i skopiować współrzędne. Nazwijmy to ,,rozwiązaniem 1'', rozw1.  *)

rozw1={{2.6288998620699893`,1.4623600453562577`},{2.7976141447365372`,1.7270944108141286`},{2.983781629058245,1.9752828784308818`},{3.1117717745294193`,2.2400172438887527`},{3.2339441861155405`,2.5378434050288554`},{3.344481129931555,2.8687613618511936`},{3.4666535415176756`,3.1996793186735317`}};

(* Aby dostać tablę wartości funkcji, zwracającej rozwiązanie (x-ową, czyli pierwszą współrzędną przecięcia) można użyć wielu sposobów. Poniżej jeden z nich. NAjpierw generujemy tabelę wartości parametru p, identyczną jak ta użyta do wygenerowania wykresów. NAstepnie, z listy rozw1 wybieramy wszystkie pierwsze elementy. Składnię [[1;;,1]] czytamy: wybierz elementy od 1 do (";;") końca (nic), ale potem z każdego wybierz pierwszy (,1) *)

{Table[p,{p,1,4,1/2}],rozw1[[1;;,1]]}

(* Teraz chcemy narysować wykres naszej funkcji, ale tabela powyżej ma zły kształt. Musimy zamienić wiersze z kolumnami, co robi funkcja Transpose *)

rozw1DAT=Transpose[%]

ListPlot[rozw1DAT,Joined->True,PlotMarkers->"OpenMarkers",AxesLabel->{"p","x"},PlotLabel->"Rozwiązanie równania  w zależności od parametru p"]

(* Teraz wejdziemy na wyższy poziom automatyzacji. Zamieniamy Table na Manipulate, równocześnie podmieniając nazwę parametru z p->P, co jest trickiem wymuszającym zmianę wartości p na wykresach wewnątrz Manipulate  *)

Manipulate[Plot[{lhs/.p->P,rhs},{x,0,4},PlotStyle->{Red,Blue}],{P,1,4,1/2}]

(* Zamiast wyznaczać przecięcie ręcznie, użyjemy FindRoot *)

Manipulate[FindRoot[Evaluate@rownanie/.p->P,{x,3}],{P,1,4,1/2}]

(* Teraz połączymy oba aplety, aby upewnić się, że FindRoot wyznacza prawidłowo punkt przecięcia.  *)

Manipulate[Plot[{lhs/.p->P,rhs},{x,0,4},PlotStyle->{Red,Blue},GridLines->{{x/.FindRoot[Evaluate@rownanie/.p->P,{x,3}]},None}],{P,1,4,1/2}]

(* Teraz wykonamy w pełni automatyczne rozwiązywanie równania *)

Table[{P,x/.FindRoot[Evaluate@rownanie/.p->P,{x,3}]},{P,1,4,1/2}]

(* Ponieważ komputer pracuje za nas, możemy dowolnie zmieniać zakres i liczbę punktów *)

rozw2DAT=Table[{P,x/.FindRoot[Evaluate@rownanie/.p->P,{x,3}]},{P,0,32,1/4}]

ListPlot[rozw2DAT,Joined->False,PlotMarkers->".",AxesLabel->{"p","x"},PlotLabel->"Rozwiązanie równania  w zależności od parametru p"]

(* Z wyprodukowanej tabeli liczb możemy utworzyć funkcję interpolującą *)

f=rozw2DAT//Interpolation

(* Którą dalej możemy używać prawie jak wbudowanych f. typu sinus czy logarytm *)

Plot[{f[p],2+Log[1+p],2+Sqrt[p]},{p,0,32}]

(* Inne, potencjalnie użyteczne sposoby rozwiązania zadania (bez komentarza) *)

g=Table[{x Log[Gamma[x]],x},{x,2,6,1/4}]//N//Interpolation

Plot[{f[p],g[p]},{p,0,12}]

Solve[E^(p/x)==Gamma[x],p]

Plot[{f[p],InverseFunction[# Log[Gamma[#]]&][p]},{p,0,12}]