“Jeden obraz wart jest więcej niż tysiąc słów”
czyli trochę o “wizualizacji”

??*Plot*

(* Choćby z wyniku tej jednej instrukcji jest jasne, jak wiele pojęć w programie Mathematica jest związanych z tworzeniem rysunków *)

Clear[z,x,y];

Rysowanie linii na płaszczyźnie

Najprostsza sytuacja występuje, gdy krzywa jest wykresem funkcji: y = f(x)

MySin[z_]:=(E^(I*z)-E^(-I*z))/(2*I)

Plot[{Sin[x],MySin[x]},{x,0,4},PlotLegends->"Expressions"]

Krzywa płaska może być zadana równaniami parametrycznymi: x=x(t), y=y(t).
Dokładnie tak opisujemy ruch punktu materialnego na płaszczyźnie, a parametr t oznacza czas.

ParametricPlot[{Cos[t],Sin[t]},{t,0,11*Pi/8}]

Do rysowania krzywych możemy używać nie tylko współrzędnych kartezjańskich - czasem wygodniej jest użyć na przykład współrzędnych biegunowych. Wtedy odpowiednikiem związku y=f(x) jest ρ=ρ(ϕ)

(* Te krzywe występują w problemie Keplera, gdzie rozważany jest układ dwóch obiektów izolowanych od reszty wszechcświata, które oddziałują siłami opisanymi przez prawo powszechnego ciążenia. W zależności od warunków początkowych mamy różne krzywe stożkowe: *)

(* elipsa *)

PolarPlot[1/(1-(1/2)*Cos[ϕ]),{ϕ,0,2*Pi}]

(* parabola *)

PolarPlot[1/(1-Cos[ϕ]),{ϕ,0.2,2*Pi-0.2}]

(* hiperbola *)

PolarPlot[1/(1-2Cos[ϕ]),{ϕ,0,2*Pi}]

Rysowanie krzywych zadanych w sposób uwikłany: F(x,y) = 0.

ContourPlot[(Sqrt[(((-1+x))^2)+(((-1+y))^2)])-2==0,{x,-3/2,4},{y,-3/2,4},PlotRange->{{-3/2,4},{-3/2,4}},Axes->True,AxesLabel->{"x","y"},AspectRatio->Automatic]

Rysowanie obszarów zadanych warunkami postaci F(x,y)  ≤ 0.

RegionPlot[x^2+3*y^4-4≤ 0,{x,-5/2,5/2},{y,-5/2,5/2},Axes->True,AxesLabel->{"x","y"},AspectRatio->1]

(* W powyższym przykładzie wystąpił jeden warunek, ale to można zmienić i narzucić więcej warunków na poszukiwany obszar *)

RegionPlot[x^2+3*y^4-4≤ 0&&x+4*y≤0,{x,-5/2,5/2},{y,-5/2,5/2},Axes->True,AxesLabel->{"x","y"},AspectRatio->1]

Rysowanie zbioru punktów na płaszczyźnie

pkty={{-1,0},{0,-1},{0,1},{1,0}, {0,0}};

ListPlot[pkty,AspectRatio->1,AxesLabel->{"x","y"},PlotStyle->{Red,PointSize[Large]}]

Ten sam efekt przy pomocy liczb zespolonych

liczbyzespolone=Table[pkty[[i]][[1]]+I*pkty[[i]][[2]],{i,1,Length[pkty]}]

ComplexListPlot[liczbyzespolone,AspectRatio->1,AxesLabel->{"x","y"},PlotStyle->{Red,PointSize[Large]}]

Ten sam zbiór punktów można wykorzystać do narysowania łamanej

ListLinePlot[pkty,AspectRatio->1,AxesLabel->{"x","y"},PlotStyle->{Red,Dashed,Thickness[0.01]}]

ClearAll["Global`*"]

Rysowanie krzywych i powierzchni w przestrzeni

Krzywa w przestrzeni może być zadana równaniami parametrycznymi: x=x(t), y=y(t), z=z(t).
Tak opisujemy ruch punktu materialnego w przestrzeni, a parametr t oznacza czas.

plot1=ParametricPlot3D[{t,2t,5Sin[t]},{t,-10,20},PlotStyle->{Black,Thick},AxesLabel->{"x","y","z"},PlotRange->{{-10,20},{-15,15},{-6,6}}]

Powierzchnie możemy zadawać na bardzo wiele sposobów. Odpowiednikiem y = f(x)  jest zależność z=g(x,y)

plot2=Plot3D[(x+2*y)/10,{x,-10,20},{y,-4,10},AxesLabel->{"x","y","z"}]

(* Składamy na jednym rysunku powierzchnię i linię *)

Show[plot1,plot2]

(* Zapytajmy, w jakich punktach linia przebija wybrany fragment płaszczyzny *)

points=Reduce[z==(x+2*y)/10&&x==t&&y==2t&&z==5Sin[t]&&-10≤ x≤ 20&&-4≤ y≤ 10,{x,y,z,t},Reals]

points[[1]][[2]][[2]]

P1=({x,2x,(1/10) (x+2 y)}/.{y->2x})/.{x->0}

P2=({x,2x,(1/10) (x+2 y)}/.{y->2x})/.{x->points[[1]][[2]][[2]]}

N[P2]

(* Rysujemy punkty P1 i P2 *)

plot3=ListPointPlot3D[{P1,P2},PlotRange->{{-10,20},{-15,15},{-6,6}},PlotStyle->{Red,PointSize[0.02]}]

(* Składamy na jednym rysunku powierzchnię, linię i pojedyncze punkty *)

Show[plot1,plot2,plot3]

ClearAll["Global`*"]

Nie zawsze rysowanie powierzchni przy pomocy zależności z=g(x,y) jest wygodne. Widać to choćby dla sfery, którą trzeba poskładać z górnej i dolnej części

r=4;Plot3D[Sqrt[r^2-x^2-y^2],{x,-r,r},{y,-r,r},AxesLabel->{"x","y","f(x,y)"}]

r=4;Plot3D[-Sqrt[r^2-x^2-y^2],{x,-r,r},{y,-r,r},AxesLabel->{"x","y","f(x,y)"}]

r=4;Plot3D[{-Sqrt[r^2-x^2-y^2],Sqrt[r^2-x^2-y^2]},{x,-r,r},{y,-r,r},AspectRatio->1,PlotRange->{{-r-0.1,r+0.1},{-r-0.1,r+0.1},{-r-0.1,r+0.1}}]

Parametryzacja powierzchni w przestrzeni wymaga podania x, y i z jako funkcji dwóch parametrów. Przydaje się znajomość współrzędnych krzywoliniowych

(* część powierzchni bocznej walca *)

r=4;ParametricPlot3D[{r*Cos[ϕ],r*Sin[ϕ],zp},{zp,-3,3},{ϕ,Pi/3,Pi},AspectRatio->1,AxesLabel->{"x","y","z"}]

ClearAll["Global`*"]

(* fragment sfery *)

r=4;ParametricPlot3D[{r*Sin[θ]*Cos[ϕ],r*Sin[θ]*Sin[ϕ],r*Cos[θ]},{θ,0,Pi/3},{ϕ,Pi/3,2Pi},AspectRatio->Automatic,AxesLabel->{"x","y","z"}]

(* cała sfera *)

r=4;ParametricPlot3D[{r*Sin[θ]*Cos[ϕ],r*Sin[θ]*Sin[ϕ],r*Cos[θ]},{θ,0,Pi},{ϕ,0,2*Pi},AspectRatio->Automatic,AxesLabel->{"x","y","z"}]

ClearAll["Global`*"]

Dla powierzchni mamy odpowiednik ContourPlot, który pozwala korzystać z przepisu F(x,y,z) = 0

(* i jeszcze jeden sposób narysowania sfery ... *)

ContourPlot3D[x^2+y^2+z^2-4==0,{x,-2,2},{y,-2,2},{z,-2,2},AxesLabel->{"x","y","z"}]

(* bardziej skomplikowana powierzchnia *)

ContourPlot3D[Abs[z-1]==x^2+y^2,{x,-2,2},{y,-2,2},{z,0,2},AxesLabel->{"x","y","z"}]

(* kula jako przykład obszaru w przestrzeni *)

RegionPlot3D[x^2+y^2+z^2≤4,{x,-2,2},{y,-2,2},{z,-2,2},Axes->True,AxesLabel->{"x","y","z"},AspectRatio->1]

(* część kuli *)

RegionPlot3D[x^2+y^2+z^2≤4&&z>Sqrt[x^2+y^2],{x,-2,2},{y,-2,2},{z,0,2},Axes->True,AxesLabel->{"x","y","z"},AspectRatio->1]

ClearAll["Global`*"]

Jeszcze inny sposób rysowania powierzchni, explicite we współrzędnych sferycznych. Jet to odpowiednik PolarPlot we współrzędnych biegunowych na płaszczyźnie.

r=1;SphericalPlot3D[r,{θ,0,Pi},{ϕ,0,2*Pi},AspectRatio->Automatic,AxesLabel->{"x","y","z"}]

SphericalPlot3D[Sin[θ],{θ,0,Pi},{ϕ,0,2*Pi},AspectRatio->Automatic,AxesLabel->{"x","y","z"}]

p=SphericalPlot3D[Sin[θ]*(2-Cos[ϕ]),{θ,0,Pi},{ϕ,0,2*Pi},AspectRatio->Automatic,AxesLabel->{"x","y","z"}]

Ważne: tworzone rysunki można eksportować do plików w różnych formatach graficznych

SetDirectory[NotebookDirectory[]];

Export["p.jpg",p]

Export["p.gif",p]

Mathematica® ma wbudowane pewne obiekty graficzne typu punkt, łamana, łuk, okrąg, wielobok, koło, ... . Mogą one występować w wersji 2D lub 3D. Ponieważ można je animować, więc z powodzeniem mogą służyć do modelowania prostych układów fizycznych. 
Typowy sposób postępowania:
obiekt → Graphics [{opcje, obiekt }] → Show → Manipulate → Animate

Przykład 1: ruch kulki w jednym wymiarze z zależnością czasową opisaną funkcją Sin[t]

ClearAll["Global`*"]

podloga1=Line[{{-1,-2},{1,-2}}];

gpodloga1=Graphics[{Thick,podloga1}];

x1[t_]:=0;y1[t_]:=Sin[t];

kulka1[t_]:=Disk[{x1[t],y1[t]},1]

gkulka1[t_]:=Graphics[kulka1[t]];

Manipulate[Show[{gpodloga1,gkulka1[t]},PlotRange->{{-2,2},{-3,2}}],{t,0,100}]

Animate[Show[{gpodloga1,gkulka1[t]},PlotRange->{{-2,2},{-3,2}}],{t,0,100}]

ClearAll["Global`*"]

h=20; (* wysokość środka kulki *)

r=1;

g=10; (* m/s^2 *)

x[t_]:=0;y[t_]:=h-(1/2)*g*t^2;

podloga=Line[{{-1,0},{1,0}}]

gpodloga=Graphics[{Thick,podloga}];

kulka[t_]:=Disk[{x[t],y[t]},r]

gkulka[t_]:=Graphics[kulka[t]];

tmax=Sqrt[2*(h-r)/g];

Animate[Show[{gpodloga,gkulka[t]},PlotRange->{{-2,2},{-2,h+2}}],{t,0,tmax}]

vmax=g*tmax;

(* t in (0,tmax)
y(t)=h-(1/2)*g*t^2 ;
t in (tmax2*tmax)
y(t)= vmax*(t-tmax)-(1/2)*g*(t-tmax)^2;
*)

ClearAll["Global`*"]

Przykład 2: pionowy ruch kulki w jednym wymiarze w jednorodnym ziemskim polu grawitacyjnym: kulka spada z wysokości h i odbija się sprężyście od podłogi.

h=20;r=1;g=10;tmax=Sqrt[2*(h-r)/g];

N[tmax]

f[t_]:=
Module[{t1=t,t2=tmax},
While[t1>2*t2,t1=t1-2*t2];Piecewise[{{h-(1/2)*g*t1^2,t1≤t2}},r+g*t2*(t1-t2)-(1/2)*g*(t1-t2)^2]]

p1=Plot[h-(1/2)*g*t^2,{t,0,tmax},PlotRange->{{0,5},{0,20}},PlotStyle->Brown]

p2=Plot[r+g*tmax*(t-tmax)-(1/2)*g*(t-tmax)^2,{t,tmax,2*tmax},PlotRange->{{0,5},{0,20}},PlotStyle->Brown]

p3=Plot[f[t],{t,0,5}]

Show[p1,p2,p3]

x[t_]:=0;y[t_]:=f[t];

podloga=Line[{{-1,0},{1,0}}]

gpodloga=Graphics[{Thick,podloga}];

kulka[t_]:=Disk[{x[t],y[t]},r]

gkulka[t_]:=Graphics[kulka[t]];

tmax=Sqrt[2*(h-r)/g];

Animate[Show[{gpodloga,gkulka[t]},PlotRange->{{-2,2},{-2,h+2}}],{t,0,10*tmax}]

ClearAll["Global`*"]

Przykład 3: żonglowanie kilkoma kolorowymi kulkami

R=10;r=2;

ϕ[t_]:=t;

kolor={Blue,Red,Green,Black,Brown,Yellow,Pink,Purple,Orange};

n=Length[kolor];

x[t_,i_]:=R*Cos[ϕ[t]+(i-1)/n*2*Pi]

y[t_,i_]:=R*Sin[ϕ[t]+(i-1)/n*2*Pi]

k[t_,i_]:=Disk[{x[t,i],y[t,i]},r]

gk[t_,i_]:=Graphics[{kolor[[i]],k[t,i]}]

gkulki[t_]:=Table[gk[t,i],{i,1,n}]

a=Animate[Show[gkulki[t],PlotRange->{{-R-2*r,R+2*r},{-R-2*r,R+2*r}}],{t,0,10}]

SetDirectory[NotebookDirectory[]];

Export["kolorowe_kulki.avi",a,"FrameRate"->3];

ClearAll["Global`*"]

Prosta konstrukcja “ruchomej” sprężyny. Konstrukcja jest zamknięta w module, by można ją było bezpiecznie używać w dowolnym otoczeniu.

gsprezyna[Pp0_,Pk0_,skala0_]:=Module[{Pp=Pp0,Pk=Pk0,skala=skala0,h,t,a,tab1,d,theta,i,tab5,tab6,Sprezyna6},
h[t_,a_,skala_]:=If[t≥a/18&&t≤17*a/18,skala*Sin[18*Pi/a*(t-a/18)],0];
tab1[a_,skala_]:=Simplify[Table[{t,h[t,a,skala]},{t,0,a,a/36}],a>0];
d[Pp_,Pk_]:=Sqrt[(Pk-Pp).(Pk-Pp)];
theta[Pp_,Pk_]:=ArcTan[Pk[[1]]-Pp[[1]],Pk[[2]]-Pp[[2]]];
tab5[Pp_,Pk_,skala_]:=Table[RotationMatrix[theta[Pp,Pk]].tab1[d[Pp,Pk],skala][[i]],{i,1,37}];
tab6[Pp_,Pk_,skala_]:=Transpose[{Transpose[tab5[Pp,Pk,skala] ][[1]]+Pp[[1]],Transpose[tab5[Pp,Pk,skala]][[2]]+Pp[[2]]}];
Sprezyna6[Pp_,Pk_,skala_]:=Line[tab6[Pp,Pk,skala]];Graphics[Sprezyna6[Pp,Pk,skala]]
]

gsprezyna[{1,1},{10,10},1]

gsprezyna[{1,11},{10,0},1/2]

ClearAll["Global`*"]

Rozważymy teraz prosty układ fizyczny: punktowa masa 𝓂 jest zawieszona na sprężynie o stałej sprężystości k i waha się w ustalonej płaszczyźnie pionowej

Elementy tego układu bardzo łatwo jest narysować. Masę reprezentuje czarne kółko, a sprężynę udaje zmieniająca swoją długość i orientację w przestrzeni łamana. Dodatkowe  elementy, na przykład “sufit” z punktem zawieszenia też łatwo można dodać. 
Ale to nie wystarcza, by stworzyć realistyczny opis zachowania tego układu. Potrzebujemy praw fizycznych !
Te z kolei przybierają postać układów równań różniczkowych i dlatego potrzebujemy funkcji DSolve i NDSolve !

Skąd wziąć równania różniczkowe rządzące ruchem układu ?
Jeśli siły działające w układzie są łatwe do określenia, to wystarczy II zasada dynamiki Newtona.
I tak byłoby w przypadku naszej masy wahającej się na (idealnej) sprężynie w jednorodnym ziemskim polu grawitacyjnym.
Jednak istnieje wiele zagadnień, gdzie siły są bardzo skomplikowane, na przykład w związku z ograniczenia narzuconymi na ruch (więzy). Wtedy warto skorzystać z formalizmu Lagrange’a, w którym kluczową rolę odgrywa energia kinetyczna i potencjalna układu oraz zbudowana z nich funkcja zwana lagrangianem !
Ten formalizm służy do znajdowania równań różniczkowych opisujących zachowanie układu
i łatwo może być zaimplementowany w programie Mathematica. 
Formalizm Lagrange’a stanowi część tzw. mechaniki teoretycznej i bardzo skrótowo przedstawię jego cechy w specjalnie przygotowanej prezentacji.
Mam nadzieję, że po tym wprowadzeniu przykład z tego notatnika i wiele innych przykładów dostępnych na stronie przedmiotu staną się zrozumiałe i przekonają Państwa o wielkim znaczeniu i wielkiej użyteczności tego podejścia do mechaniki.

x[t]=l[t]*Sin[ϕ[t]];

y[t]=-l[t]*Cos[ϕ[t]];

x'[t]=D[x[t],t];

y'[t]=D[y[t],t];

(* Energia kinetyczna punktowej masy *)

Ekin=FullSimplify[(1/2)*m*(x'[t]^2+y'[t]^2)]

(* Energia potencjalna zależy od y i rozciągnięcia sprężyny *)

Epot=  FullSimplify[m*g*y[t] +(1/2)*k*(l[t]-l0)^2]

(* Funkcja Lagrange'a układu *)

L=FullSimplify[Ekin-Epot]

Etot=FullSimplify[Ekin+Epot]

d1=D[L,ϕ'[t]];

d2=D[L,ϕ[t]];

d3=D[d1,t];

d4=D[L,l'[t]];

d5=D[L,l[t]];

d6=D[d4,t];

s=Solve[{d3==d2,d6==d5},{ϕ''[t],l''[t]}];

FullSimplify[s]

eq1=Flatten[FullSimplify[s]/.Rule->Equal]

(* Parametry układu: długość i masa wahadła matematycznego oraz wartość przyspieszenia ziemskiego *)

k=10; (* N/m *)

l0=1/5 ; (* m *)

m=1/4; (* kg *)

g=981/100; (* m/s^2 *)

(* Warunki początkowe  *)

li=1/4; (* m *)

vi=0; (* m/s *)

(Subscript[ϕ, 0])=Pi/3;

(Subscript[ω, 0])=0;

tmax=20; (* s *)

WP={l[0]==li,l'[0]==vi,ϕ[0]==(Subscript[ϕ, 0]),ϕ'[0]==(Subscript[ω, 0])}

numerycznie=NDSolve[{eq1,WP},{l[t],ϕ[t]},{t,0,tmax},WorkingPrecision->24]

Flatten[numerycznie]

(* Wyciągam funkcję ϕ[t[ z numerycznego rozwiązania *)

l[t_]:=Evaluate[l[t]/.Flatten[numerycznie]]

ϕ[t_]:=Evaluate[ϕ[t]/.Flatten[numerycznie]]

(* Teraz mogę używać funkcji ϕ[t] tak, jak innych funkcji *)

l[1]

l'[2]

ϕ[3]

ϕ'[11]

x[t_]:=l[t]*Sin[ϕ[t]];

y[t_]:=-l[t]*Cos[ϕ[t]];

Plot[ϕ[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","ϕ(t) [rad]"}]

Plot[ϕ'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","ϕ'(t) [rad/s]"}]

Plot[l[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","l(t) [m]"}]

Plot[l'[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","l'(t) [m/s]"}]

Plot[Evaluate[Ekin],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Ekin(t) [J]"}]

Plot[Evaluate[Epot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Epot(t) [J]"}]

Plot[Evaluate[Etot],{t,0,tmax},PlotRange->All,AxesLabel->{"t [s]","Etot(t) [J]"}]

(* Tu używam modułu dla sprężyny *)

gsprezyna[Pp0_,Pk0_,skala0_]:=Module[{Pp=Pp0,Pk=Pk0,skala=skala0,h,t,a,tab1,d,theta,i,tab5,tab6,Sprezyna6},
h[t_,a_,skala_]:=If[t≥a/18&&t≤17*a/18,skala*Sin[18*Pi/a*(t-a/18)],0];
tab1[a_,skala_]:=Simplify[Table[{t,h[t,a,skala]},{t,0,a,a/36}],a>0];
d[Pp_,Pk_]:=Sqrt[(Pk-Pp).(Pk-Pp)];
theta[Pp_,Pk_]:=ArcTan[Pk[[1]]-Pp[[1]],Pk[[2]]-Pp[[2]]];
tab5[Pp_,Pk_,skala_]:=Table[RotationMatrix[theta[Pp,Pk]].tab1[d[Pp,Pk],skala][[i]],{i,1,37}];
tab6[Pp_,Pk_,skala_]:=Transpose[{Transpose[tab5[Pp,Pk,skala] ][[1]]+Pp[[1]],Transpose[tab5[Pp,Pk,skala]][[2]]+Pp[[2]]}];
Sprezyna6[Pp_,Pk_,skala_]:=Line[tab6[Pp,Pk,skala]];Graphics[Sprezyna6[Pp,Pk,skala]]
]

gs[t_,skala_]:=gsprezyna[{0,0},{x[t],y[t]},skala];

sufit=Rectangle[{-1/2,0},{1/2,1/20}];

masa[t_]=Disk[{x[t],y[t]},1/30];

gsufit=Graphics[sufit];

gmasa[t_]:=Graphics[masa[t]]

Animate[Show[gsufit,gs[t,1/30],gmasa[t],PlotRange->{{-6/5,6/5},{-6/5,1/5}},AspectRatio->Automatic],{t,0,tmax}]