Matlab Przykłady praktyczne Równania różniczkowe

1. MOiPP Wykład 5 • Matlab • Przykłady praktyczne • Równania różniczkowe

2. Wykresy funkcji symbolicznych - ezplot syms x f=x^2 fp=diff(f) subplot(2,1,1) ezplot(f, [-2,2]) subplot(2,1,2) ezplot(fp, [-2,2])

3. Długość krzywej (np. paraboli) clc clear syms x disp('parabola') f=x^2; disp(f) d=int(sqrt(1+diff(f)^2), x, -1, 1); %całka ozn. disp('całkowita długość') disp (subs(d)) %obliczenie długości parabola x^2 d = 5^(1/2)-1/2*log(-2+5^(1/2)) całkowita długość 2.9579 3

4. Długość krzywej paraboli – równania parametryczne clc, clear symst %równania parametryczne paraboli x=t/2 y=t^2/4 l=int(sqrt(diff(x)^2+diff(y)^2),t,-2, 2) L=subs(l) %oblicz wartość liczbową t=-2:0.1:2; x=subs(x); y=subs(y); plot(x,y) UWAGA: inne granice całkowania! (t) x = 1/2*t y = 1/4*t^2 l = 5^(1/2)-1/2*log(-2+5^(1/2)) L = 2.9579

5. Można też liczyć w sposób przybliżony metodą siecznych Metoda: Sumujemy odcinki prostej pomiędzy punktami węzłowymi x=-1:0.1:1; s=0; for k=1:length(x)-1 y(k)=x(k).^2; y(k+1)=x(k+1).^2; s=s+sqrt((y(k+1)-y(k)).^2+(x(k+1)-x(k)).^2); end disp(s) Długość odcinka wg prawa Pitagorasa s = 2.9564 wynik przybliżony

6. Długość f(x)=sin(x) w przedziale (0, pi) t=-1:0.1:pi; s=0; for k=1:length(t)-1 s=s+sqrt((sin(t(k+1))-sin(t(k)))^2+(t(k+1)-t(k))^2); end disp(s) s = 3.7609

7. Zaspa Iloczyn dwóch paraboli

8. clear %wykres zaspy a=10;b=20;h=0.5; [x,y]=meshgrid(-a/2:0.1:a/2, -b/2:0.1:b/2); z = h*(a^2-4*x.^2).*(b^2-4*y.^2)/b^2/a^2; mesh(x,y,z) %obliczenie objętosci syms y x h a b zsym= h*(a^2-4*x^2)*(b^2-4*y^2)/b^2/a^2; p=int(int(zsym,x,-a/2,a/2),y,-b/2,b/2) a=10;b=20;h=0.5; p0=subs(p)

10. Wyznaczanie ekstremum funkcji – badanie pochodnych clc,clear syms x f=x^3-5*x+2 ezplot(f,[-3, 3]) %rysujemy wykres funkcji df=diff(f,x);% 1-sza pochodna ext=subs(solve(df,x)) %znajdujemy zera pochodnej (ekstrema) d2f=diff(f,x,2); % 2-ga pochodna for k=1:length(ext) x=ext(k); % wstawiamy ekstremum do x – teraz x jest liczbą inf=subs(d2f); % wartość 2-giej pochodnej w punkcie extremum if imag(inf)==0 % jeśli wartość nie jest zespolona if inf<0 fprintf('w punkcie %f jest maksimum\n', x) else fprintf('w punkcie %f jest minimum\n', x) end end end;

11. Wyniki badania extremum f = x^3-5*x+2 ext = 1.2910 -1.2910 w punkcie 1.290994 jest minimum w punkcie -1.290994 jest maksimum

12. Rozwiązywanie równań różniczkowych • funkcja dsolve() . Funkcja ta oblicza symbolicznie rozwiązania równań różniczkowych zwyczajnych. Równania są określane przez symboliczne wyrażenia zawierające literę D do oznaczenia stopnia. Symbole D2, D3... DN, odnoszą się do drugiej, trzeciej,..., n-tej pochodnej. D2y jest zatem odpowiednikiem symbolicznym Zmienna niezależna domyślna to t.

13. UWAGI Nazwy zmiennych symbolicznych nie powinny zawierać D. Zmienną niezależną domyślną t można zmienić i podać jako drugi argument. Warunki początkowe mogą być określone przez dodatkowe równania. Jeśli nie określono warunków początkowych, rozwiązania zawierają stałe całkowania: C1, C2, itp.

14. Przykłady dsolve('Dy=1+y^2') %zmienna t domyślna • ans = • tan(t+C1) dsolve('Dy=1+y^2','x')%zmienna x ustalona ans = tan(x+C1) parametr symboliczny dsolve('Dx = -a*x') ans = C1*exp(-a*t)

15. Po wstawieniu warunków początkowych: f = dsolve('Dy=1+y^2','y(0)=1') • f = • tan(t+1/4*pi) Uwaga: y jest w obszarze roboczym MATLAB, ale t nie jest, a zatem polecenie diff(y,t) zwraca błąd. Aby umieścić t w obszarze roboczym należy: symst pochodna=diff(y,t) pochodna= 1+tan(t+1/4*pi)^2

16. Przykład 4 warunki początkowe u = dsolve('D3u=u','u(0)=1','Du(0)=-1','D2u(0) = pi','x') • D3u reprezentuje d3u/dx3 • D2u(0) odpowiada u"(0) Wynik u = 1/3*pi*exp(x)-1/3*(pi+1)*3^(1/2)*exp(-1/2*x)*sin(1/2*3^(1/2)*x)+ (1-1/3*pi)*exp(-1/2*x)*cos(1/2*3^(1/2)*x)

17. Układ równań różniczkowych Funkcja dsolve rozwiązuje także układ równań różniczkowych zwyczajnych kilku zmiennych, z warunkami początkowymi lub bez.

18. Przykład • Dwa równania liniowe, pierwszego rzędu: • S = dsolve('Df = 3*f+4*g', 'Dg = -4*f+3*g') • S = • f: [1x1 sym] • g: [1x1 sym] • Rozwiązania obliczane są zwracane w strukturze S. • Można określić wartości f i g, wpisując: f = S.f f = exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) g = S.g g = exp(3*t)*(C1*cos(4*t)-C2*sin(4*t))

19. Jeśli chcemy uzyskać f i g bezpośrednio, oraz uwzględnić także warunki początkowe, wpisujemy: • [f,g] = dsolve('Df=3*f+4*g, Dg =-4*f+3*g', 'f(0) = 0, g(0) = 1') • f = • exp(3*t)*sin(4*t) • g = • exp(3*t)*cos(4*t)

20. Przykład 3 • Równanie różniczkowe drugiego stopnia z dwoma warunkami początkowymi: f = dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0', 'x') f = 4/3*cos(x)-1/3*cos(2*x) f2=simplify(f) %uproszczenie f2 = 4/3*cos(x)-2/3*cos(x)^2+1/3

21. Sprawdzenie rozwiązania s= diff(diff(f2))-cos(2*x)+f2 s2=simplify(s) s = -4/3*sin(x)^2+2/3*cos(x)^2-cos(2*x)+1/3 s2 = 0

22. Jeszcze jeden przykład składni w SymbolicMathToolbox. warunek pocz.: y(0) = 1 • y = dsolve('Dy+4*y = exp(-t)', 'y(0) =1') • spr=diff(y,t)+4*y %sprawdzenie rozwiązania • spr1 =simplify(spr) %uproszczenie wyrażenia • t=0%sprawdzenie warunku początkowego • wp=subs(y) • y = • (1/3*exp(3*t)+2/3)*exp(-4*t) • spr = • exp(-4*t)*exp(3*t) • spr1 = • exp(-t) • t = • 0 • wp = • 1

23. Podsumowanie Podstawowe funkcje SymbolicTool syms– inicjacja zmiennych symbolicznych solve – symboliczne rozwiązywanie równań subs – podstawianie danych do wyrażeń symbolicznych limit – wyznaczanie granic funkcji (limes) diff - wyznaczanie pochodnych funkcji int - wyznaczanie całek nieoznaczonych i oznaczonych dsolve – rozwiązywanie równań różniczkowych simplify – upraszczanie wyrażeń symbolicznych ezplot– wykres funkcji symbolicznej

24. Przykłady zastosowań równań różniczkowych 1. Rzut ukośny pod kątem 0z prędkością początkową v0 (pomijamy opór powietrza) warunki początkowe ax(t) = = 0 vx(0) =x'(0) = v0 cos 0 x(0)=0 zapis w Matlabie: = -g ay(t) = vy(0) =y'(0) = v0 sin 0 y(0)=0 x1 = dsolve('D2x = 0','Dx(0) = v0*cos(alfa0)','x(0)=0','t') y1 = dsolve('D2y = -g','Dy(0) = v0*sin(alfa0)','y(0)=0','t')

25. clear,clc syms a t0 t v0 v g alfa0 x1 = dsolve('D2x = 0','Dx(0) = v0*cos(alfa0)','x(0)=0','t') y1 = dsolve('D2y = -g','Dy(0) = v0*sin(alfa0)','y(0)=0','t') alfa0=pi/4; v0=10; g=9.81; %-------------------------położenie y(x) % wstawienie danych do wzorów xu1=subs(x1) yu1=subs(y1) % -----------------------wektory dla punktów czasu t=0:0.01:2; xp1=subs(xu1) yp1=subs(yu1) subplot(3,1,1) plot(xp1,yp1) title('Położenie y(x)') cd.

26. %---------------------------składowe prędkości vx(t) vy(t) symst vx=diff(xu1) vy=diff(yu1) t=0:0.01:2; %wektor punktów czasu vx=subs(vx) %wstawienie danych vx1=ones(1,length(t))*vx; % vx=const - wektor 1*vx vy1=subs(vy); % wektor vy dla punktów czasu subplot(3,1,2) plot(t,vx1,t,vy1) title('Składowe prędkoscivx(t) vy(t)') % ----------------------------------wyznaczenie s(t) syms t s=int(sqrt(diff(xu1,t)^2+diff(yu1,t)^2),t,0,t) %całka oznaczona t=0:0.01:2; st=subs(s); subplot(3,1,3) plot(t,st) title('s(t)')

28. 2. Oscylator harmoniczny bez tłumienia • gdzie: • y(t) – położenie ciała, – częstość drgań, • 3. Oscylator z tłumieniem – współczynnik tłumienia.

29. clc syms x f=dsolve('D2y+1*y=0','y(0)=4','Dy(0)=0','x') t=0 subs(f) ezplot(f,[0 10*pi, -6 6],1) f2=dsolve('D2y+0.1*Dy+1*y=0','y(0)=4','Dy(0)=0','x') ezplot(f2,[0, 10*pi, -5 5],2)