%% Działania w Matlabie na wyrażeniach symbolicznych cz. 2
% Opracował dr. inż Robert Jakubowski, Politechnika Rzeszowska, KSiSL
% Table of Contents
% Różniczkowanie
% Całkowanie
%     Obliczanie pola powierzchni z wykorzystaniem całki oznaczonej
%     Obliczanie pola powierzchni z wykorzystaniem podwójnej całki oznaczonej
%         Wynacznie pola łuku koła oraz pola koła o zadanym promieniu
%     Obliczanie objętości z wykorzystaniem całki oznaczonej podójnej
%         Ćwiczenie z wyznaczania objętości
%         Ćwiczenie z wykorzystaniem całki potrójnej do obliczeń objętości
% Równania różniczkowe
%     Przykłady zastosowania równań różniczkowych do układów mechanicznych
%         Ćwiczenia z zakresu analizy układu kinematycznego masy drgającej
%% Różniczkowanie
% Obliczenia pochodnych z wyrażeń algebraicznych można realizować funkcją diff. Przykładowe  obliczenia wartości pierwszej pochodnej z wyrażenia:
% y=3*x^2 +2*x-e^x

syms x
y=3*x^2+2*x-exp(x)
dydx=diff(y)

% Wyznaczenie drugiej pochodnej z wyrażenia y można zrealizować na dwa sosoby. Pierwszym jest wyznaczenie wartości pierwszej pochodnej, co już zrobiliśmy, a następnie wyznaczenie wartości pochodnej z pochodnej czyli:

d2ydx=diff(dydx)

% Drugim sposobem jest wyznaczenie od razu wartości drugiej pochodnej z wyrażenia y w następujący sposób

d2ydx_2=diff(y,2)

% Porównując otrzymane obydwoma metodami wyniki są one takie same.
% Gdybyśmy chcieli wyznaczyć trzecią pochodną z wyrażenia y to:

d3ydx=diff(y,3)

% Chcąc obliczyć wartość pochodnej w punkcie, można najpierw wyznaczyć określoną pochodną, a następnie w miejsce zmiennej podstawić określoną wartość wykorzystując funkcje subs i double

d2ydx_5=double(subs(d2ydx,x,5))

% W funkcji diff można też wskazać zmienną, po której będziemy relizować różniczkowanie. Jest to istotne wtedy, gdy mamy więcej zmiennych w wyrażeniu. W takim wypadku na drugiej pozycji wpisujemy tą zmienną po której chcemy różniczkować.

syms x y a
fxy=a*x^3+2*x^2*y-y^3-2*a*y

% Pierwsza i druga pochodna funkcji po zmiennej x będzie:

dfdx=diff(fxy,x)
dfdx2=diff(fxy,x,2)

% Zrobienie tego samego dla zmiennej y:

dfdy=diff(fxy,y)
dfdy2=diff(fxy,y,2)

%
%% Całkowanie
% Całkowaanie wyrażeń algebraicznych z wykorzystaniem obliczeń symbolicznych realizuje się wykorzystując funkcję int. Przykładowe obliczenie wartości całki można wykonać jak przedstawiono poniżej:
% \int \left(\sin \left(2x\right)-x^2 \right)\textrm{dx}

clear variables
syms x
y=sin(2*x)+x^2
cal_y=int(y)

% Gdy chcemy zdefiniować zmienną po której odbywa się całkowanie wtedy tą zmienną wstawia się na drugiej pozycji. Przyjżyjmy się wynikom całkowania, gdy jako zmiennną całkowania wyrażenia y wstawi się x lub z

syms z
cal_yx=int(y,x)
cal_yz=int(y,z)

% Obliczanie całki określonej w granicy od x1 do x2, gdzie x1 jest początkiem, a x2 końcem przedziału całkowania można wykonać podajac granice całkowania na kolejnych pozycjach. Chcąc wyliczyć przedstawioną poniżej całkę:i
% P=\int_0^2 \left(\sin \left(2x\right)-x^2 \right)\textrm{dx}
% w Matlebie będzie to zapisane:

P=int(y,x,0,2)
Pval=double(P)

% Wykorzystano funkcję double do przeliczenia zmiennej symbolicznej na zmienną numeryczną
% Jak powszechnie wiadomo całkę oznaczoną można wykorzystać do obliczeń pola powierzchni pod wykresem, albo objętości figury przestrzennnej.
%% Obliczanie pola powierzchni z wykorzystaniem całki oznaczonej
% Policzmy pole powierzchni pod krzywą zakreśloną przez funkcję sin(x) w przedziale od 0 do \frac{\pi }{2} .
% Zacznijmy od narysowania funkcji sin(x) w przedziale od 0 do \pi

fplot(sin(x),[0,pi])

% Pokażmy obszar który chcemy policzyć

hold on
plot([0,pi/2,pi/2],[0,0,sin(pi/2)],'r-')
fplot(sin(x),[0,pi/2],'r-')
hold off

% Wyznaczmy wartość pola pod krzywą wykorzystyjąc całkę oznaczoną, czyli:
% P=\int_0^{\frac{\pi }{2}} \sin \left(x\right)\textrm{dx}

P=double(int(sin(x),x,0,pi/2))

% Pole całego pokazanego przebiegu sinusoidy
% P=\int_0^{\pi } \sin \left(x\right)\textrm{dx}

P=double(int(sin(x),x,0,pi))

%
%% Obliczanie pola powierzchni z wykorzystaniem podwójnej całki oznaczonej
%% Wynacznie pola łuku koła oraz pola koła o zadanym promieniu
% Jeżeli w czasie całkowania wyznacza się zmianę w obydwu kierunkach x i y, to do obliczania pola wykorzystuje się całkę podwójną. Rozpatrzmy to na przykładzie obliczeń pola całego koła, lub jego części
% Rozwiązanie zadania można wykonać wykorzystyjąc opis położenia punktów x i y jako zależności od promienia r i kąta \alpha, jak pokzano na rysunku.

clear variables
syms x y
fimplicit(4^2==x^2+y^2)
axis('equal','on')
grid on
hold on
%plot max radius
plot([0,4*cos(3*pi/4)],[0 4*sin(3*pi/4)])
text(-1.6,2,'R',"Rotation",-45)
%plot radius
plot([0,4*cos(pi/4)],[0 4*sin(pi/4)])
plot([0,2],[0,2],'m')
text(1,1.3,'r',"Rotation",45)
%plot arc of ungle
xpom=cos(pi/4:-pi/4/10:0);
ypom=sin(pi/4:-pi/4/10:0);
plot(xpom,ypom,'r-')
text(1.2*cos(pi/10),1.2*sin(pi/10),'alpha')

%plot x y line
plot([2,2,0],[0,2,2],'k')
text(2.1,0.,'y=r cos(alpha)','Rotation',90)
text(0.,2.35,'x=r sin(alpha)')
hold off

%
% Zgodnie z przedstawionym rysunkiem współrzędne x i y dla zadanego kąta i długości promienia są wyrażone zależnością:
% \begin{array}{l}
% x=r*\cos \left(\alpha \right)\\
% y=r*s\textrm{in}\left(\alpha \right)
% \end{array}

syms x y r alpha 
x=r*cos(alpha)
y=r*sin(alpha)

% Oczywiście maksymalna wartość promienia zdefiniowana jest jako R.
% Definiując zmienną Prom jako długość promienia dla danej wartości r i \alpha
% Pole zakreślone przez promień (Prom) dla kąta \alpha można policzyć jako całkę podwójną ze zmiany promienia Prom względem jego długości od 0 do R i kąta od 0 do \alpha_k czyli
% P_l =\int_0^R \int_0^{\alpha_k } \textrm{Prom}\;\;d\alpha \;\textrm{dr}\;
% gdzie Prom można zdefiniować jako
% \textrm{Prom}=\sqrt{x^2 +y^2 }
% Ze względu na tok obliczeń, celowo wprowadzono zmienną Prom różną od r, wynika z tąd, że Prom jest funkcją r i \alpha.
% Załużmy, że obliczenia będę prowadzone dla 1/4 wycinka koła o promieniu R=4. Stąd zadanie sprowadza sie do rozwiązania.
% P_l =\int_0^4 \int_0^{\frac{\pi }{2}} \mathrm{Prom}\;d\alpha \;\mathrm{dr}\;

%Wyznaczamy zmienną Prom
Prom=sqrt(x^2+y^2)
%Obliczamy całkę po kącie alpha ze zmennej Prom
P_p_alpha=int(Prom,alpha,[0,pi/2])
%Obliczamy całkę po zmiennej r ze zmiennej P_p_alpha 
P_p=int(P_p_alpha,r,[0,4])
P_p_val=double(subs(P_p))

% Dla sprawdzenia wyznaczmy wartość pola ćwierci koła z klasycznej metody

Pp=pi*4^2/4

% Wynik obliczony za pomocą całki podwójnej jest zgodny z wynikiem otrzymanym ze wzoru na pole powierzchni 1/4 wycinka kołowego
% Policzmy teraz pole całego koła, czyli całkę podwujną w postaci:
% P_l =\int_0^4 \int_0^{2\pi } \textrm{Prom}\;d\alpha \;\textrm{dr}\;

%Obliczamy całkę po kącie alpha ze zmennej Prom
P_p_alpha=int(Prom,alpha,[0,2*pi])
%Obliczamy całkę po zmiennej R ze zmiennej P_p_alpha 
P_p=int(P_p_alpha,r,[0,4])
P_p_val=double(subs(P_p))
%Dla sprawdzenia wyznaczmy wartość pola ćwierci koła z klasycznej metody
Pp=pi*4^2

% Podany sposób można wykorzystać do wyznaczenia zależności ogólnej na pole koła. W tym celu musielibyśmy rozwiązać całkę podwójną gdzie granica całkowania po r jest zmienną symboliczną R:
% P_k =\int_0^R \int_0^{2\pi } \textrm{Prom}\;d\alpha \;\textrm{dr}\;
% Ponieważ całkę wewnętrzną mamy policzonną pod zmienną P_p_alpha, wiec zostaje tylko policzyć całkę zewnętrzną

syms R
P_k=int(P_p_alpha,r,[0,R])

% Wynik obliczeń dał nam dobrze znany wzór na pole koła w funkcji promienia R.
%
%% Obliczanie objętości z wykorzystaniem całki oznaczonej podójnej
% Zgodnie z definicją objętość bryły będzie wyrażona za pomocą całki podwójnej, czyli
% \int_{\textrm{y1}}^{\textrm{y2}} \int_{\textrm{x1}}^{\textrm{x2}} f\left(x,y\right)\textrm{dxdy}
% Taką metodą można obliczać pola figur gdzie zmienna z jest okreśłona jako funkcja zależna od dwóch pozostałych.
% Załóżmy, że mamy powieżchnię opisaną równaniem:
% z={3-x}^6 -y^2y
% Policzmy objętość jaką będzie stanowić bryła ograniczona od dołu płaszcyzną XY od góry powierzchnią z oraz na osi x i y przedziałem -1 do 1. Matematycznie było by to zapisane
%  \int_{-1}^1 \int_{-1}^1 z\left(x,y\right)\textrm{dxdy}
% Zapiszmy to w postaci wyrażenia symbolicznego i przedstawmy na wykresie

clear variables
syms x y
z(x,y)=3-x.^6-y.^2

% Policzmy wartości z1=z(x,y) dla x i y  w przedziale od -1 do 1.
% Przedstawienie graficzne obliczanej objętości
% Wykorzystajmy funkcję meshgrid do przygotowania tablicy z wartościami x i y i następnie te tablice podstwimy do funkcji z


[x1,y1]=meshgrid(-1:0.05:1);
[z1]=double(z(x1,y1));
z1(size(z1))=0;

% Wykorzystując funkcję meshz narysujemy interesujący nas profil, a przez dodanie punktu 0,0,0 narysujemy wycinek zaczynajcy się od płaszczyzny XY

meshz(x1,y1,z1)
xlabel('x')
ylabel('y')
zlabel('z')
axis([-1.5 1.5 -1.5 1.5 0 3])

% Wyznaczanie objętości
% Policzymy oznaczoną całkę podwójną rozbijając ją na dwie całki wyznaczane kolejno względem zmiennej x i y z wyrażena z

V_x=int(z,x,-1,1)
V=double(int(V_x,y,-1,1))

%% Ćwiczenie z wyznaczania objętości
% Policz ile będzie wynosiła objętość bryły gdy granica całkowania dla y zmieni się na przedział od 0 do 1.
% \int_0^1 \int_{-1}^1 z\left(x,y\right)\textrm{dxdy}
%  Czy wynik będzie połową objędtości wyliczonej powyżej?



%
%% Ćwiczenie z wykorzystaniem całki potrójnej do obliczeń objętości
% Na podstawie przykładu obliczeń pola koła policzyć objętość kuli jako całkę potrójną zaleną od promienia r, kąta alpha na płaszczyżnie XY i kąta theta dla określenia odchylenia promienia w przestrzeni  względem płaszczyzny XY. Wzór na objętość kuli w tym wypadku będzie przedstawiony jako:
% P_l =\int_0^4 \int_0^{2\pi } \int_0^{2\pi } \textrm{Pom}\;d\alpha \;d\theta \;\textrm{dr}\;
% Podpowiedź: do promienia (zmienna Prom) wyznaczonego na płaszczyźnie należy dodać składową określającą zmienę względem osi z zależną od r i kąta \theta.



%
%% Równania różniczkowe
% Wprowadźmy zapis równania różniczkowego w wyrażeniu symbolicznym. Potrzebne jest do tego zdefiniowanie funkcji symbolicznej np. y(x). Definiuje się ja podobnie jak zmienną poprzez funkcję syms

clear variables
syms y(x)

% Teraz mozna zdefiniować rówananie rózniczkowe zmiennej x w postaci
% \frac{\textrm{dy}}{\textrm{dx}}=-2*x+4*y\left(x\right)
% i zapisać to jako zmienną r_rx

r_rx=diff(y,x)==-2*x^2-4*y

% Do roziązywania równania rózniczkowego wykorzystuje się funkcję dsolve

y_roz=dsolve(r_rx)

% Oczywiście takie rozwiązanie spowoduje pojawienie się stałych całkowania, które można wyznaczyć wprwadzając warunki początkowe. Np podając wartość początkową  y(x) dla x=-2, np y(0)=-2

w_p=y(0)==-2

% Po wprowadzeniu warunków początkowych rozwiąznie równania sprowadza sie do wywałoania funkcji dsolve z parametrami dsolve(r_rx,w_p), czyli

y_roz_p=dsolve(r_rx,w_p)

% Teraz możemy przedstawić rozwiązanie, czyli przebieg fynkcji y(x) np dla x od 0 do 1.5

fplot(y_roz_p,[0 1.5],'b')

% Sprwdżmy jak będzie wyglądała funkcją będąca różniczką funkcji y w tym przedziale i narysujmy ją na wykresie

dydx=diff(y_roz_p,x)
hold on
fplot(dydx,[0,1.5],'r--')
grid on

% Poszukajmy extremum funkcji y dla prezentowanego przedziału, czyli znajdźmy warość x dla której pochodzna równa się 0
% czyli \frac{\mathrm{dy}}{\mathrm{dx}}=0. W tym celu, założymy że x jest w interesującym nas przedziale 0<x<1.5.

assume(x>0)
assumeAlso(x<1.5)
x0=double(solve(dydx==0,x))
legend(['y(x) ';'dy/dx';'ymax '])

% Wyznaczmy wartość y od x0 i przedstawny go na wykresie

y0=double(subs(y_roz_p,x,x0))
plot(x0,y0,'mo')
axis([0 1.5 -2 3])
hold off

%
%% Przykłady zastosowania równań różniczkowych do układów mechanicznych
% Przykład rozwiązania dla masy połączonej z podłóżem za pomocą tłumika i sprężyny.
%
% Mamy układ przedstawiony jak na rysunku. Masa m jest połączona z podłożem za pomocą sprężyny o współczynniku sprężystości k (N/m) i tłumika drgań o współczynniku tłumienia c (N s/m). Określmy przemieszceznia masy m względem osi y, gdy zadziała na niego siła F skierowana zgodnie z osią y
% Układ opisuje równanie:
% m\frac{d^2 y}{{\textrm{dt}}^2 }+c\frac{\textrm{dy}}{\textrm{dt}}+\textrm{ky}=F\left(t\right)
% W stanie równowagi, gdy na układ nie działa żadna siła równie to ma postać:
% m\frac{d^2 y}{{\textrm{dt}}^2 }+c\frac{\textrm{dy}}{\textrm{dt}}+\textrm{ky}=0
% Wprowadźmy zatem zmienne symboliczne odpowidajace parametrom równania

syms y(t) m c k F(t)

% Wprowadźmy ograniczenia na zmienne, zakładajace, że charakterystyki elementów, masa oraz przemieszczenie ma wartości dodatnie oraz że siła jest wielkością należacą do zbioru liczb rzeczywistych.

assume([k c m],'positive')
assume(y(t),'real')

% Zacznijmy od rozwiązania zadania gdy siła F=0. Wprowadźmy rónanie sił działających na układ.

RM=m*diff(y,t,2)+c*diff(y,t)+k*y==0

% Do rozwiązania zadania należy podać warunki początkowe, dotyczące punktu startu, poczatkowej prędkości czy przyśpieszenia. Załużmy, że w początkowym stanie położenie jest odchylone od stanu równowagi o -8cm y(0)=-0.08 oraz, że prędkość początkowa masy była równa 0, yp - jest zmiennną określającą początkowe położenie, vp - określa początkową prędkość, aby ją podać, najpierw należy wyznaczyć funkcję reprezentującą pochodną po y(t), stąd definicja warunków początkowych jest następująca.

yp=y(0)==-0.08;
Dy(t)=diff(y(t),t)
vp=Dy(0)==0;

% Rozwiążmy to zadanie

Y_m = dsolve(RM,[yp vp])

% Podstawmy wartości za zmienne. Załużmy, że masa m=70 kg, wspólczynnik sprężystości k=30 N/m, a współczynnik tłumienia c=15 Ns/m. Podstawmy wartości i narysujmy wykres zmian położenia masy od czasu t w przedziale od 0 do 150 s.
%

Y_m_val(t)=subs(Y_m,[m,c,k],[70,15,30])
ts=[0:0.5:100];
Y_val=double(Y_m_val(ts));
plot(ts,Y_val)
% fplot(Y_m_val,[0 100]),

% Dołużmy wykres zmian prędkości. W tym celu  zróżniczkujmy otrzymane wyrażenie na drogę Y_m_val(t). Następnie policzymy wartości prędkości dla zadanych wartości czas ts

V_y=diff(Y_m_val,t)
V_val=double(V_y(ts));
hold on
plot(ts,V_val,'g')
xlabel('t [s]')
ylabel('y [m], Vy [m/s]')
%fplot(V_y,[0 100],'g')
legend([' y ';' Vy'])
grid on
hold off

%% Ćwiczenia z zakresu analizy układu kinematycznego masy drgającej
% Sprawdzić jak zachowa się układ, gdy wymuszenie, czyli początkowe wychylenie masy od punktu równowagi będzie -4 cm. Wynik przedstawić graficznie.
% Podpowiedź: należy wykonać taką samą analizę jak poprzednio zmieniajjąc warunki początkowe.



% Wskazać jak zmieni się charakreystyka drgań masy gdy zmienimy wspólczynnik tłumienia. Zrobić ćwiczenie wstawiając trzy wartości współczynnika tłumienia np c=[10 15 20] Ns/m. Przedstawić wynik na wykresie.
%
% Zwiększyć wspólczynnik tłumienia do c=50 Ns/m



%