1 Pakiety numeryczne Równania różniczkowe Łukasz Sztangret Katedra Informatyki Stosowanej i Modelowania

2 Równania różniczkowe zwyczajne Do rozwiązywania równań różniczkowych zwyczajnych I rzędu lub układów równań różniczkowych zwyczajnych I rzędu można wykorzystać jedną z poniższych funkcji Matlaba: ode23 ode45 ode113 ode15s ode23s ode23t ode23tb Składnia: [t,y] = ode45(dy,T,y0); Uchwyt do funkcji zwracającej pochodną szukanej funkcji y (dla jednego równania) lub wektor pochodnych (dla układu równań) Przedział czasu Warunek początkowy Wektor czasu Wektor zawierający wartości szukanej funkcji y

3 Równanie różniczkowe I rzędu Wykonujemy przekształcenie aby otrzymać wzór na y’ Rozwiązanie w Matlabie: dy=@(t,y)1-y; [t,y]=ode45(dy,[0 1],2); plot(t,y); Definiujemy uchwyt do funkcji zwracającej pochodną szukanej funkcji y Wywołujemy funkcję ode45 podając dy, przedział czasu oraz warunek początkowy. Rozwiązaniem analitycznym jest funkcja

4 Układ równań różniczkowych I rzędu Rozwiązanie w Matlabie: dy=@(t,y)[y(1)*log(y(2));... -y(2)]; [t,y]=ode45(dy,[0 10],[exp(-2), exp(2)]); plot(t,y); Rozwiązaniem analitycznym są funkcje Definiujemy uchwyt do funkcji zwracającej wektor zawierający pochodne szukanych funkcji y

5 Równanie różniczkowe II rzędu Równanie należy przekształcić na układ równań I rzędu wprowadzając nowe zmienne: Podstawiając otrzymujemy: Obliczamy pochodne y 1 oraz y 2 : Rozwiązanie w Matlabie: dy=@(t,y)[y(2);... 1+t-y(1)+2*y(2)]; [t,y]=ode45(dy,[0 2],[0 0]); plot(t,y(:,1)) Rozwiązaniem analitycznym jest funkcja

6 Równanie różniczkowe I rzędu Wykonujemy przekształcenie aby otrzymać wzór na y’ Rozwiązanie w Matlabie: dy=@(t,y)-2*y; [t,y]=ode45(dy,0.5:-0.5:0,exp(-1)); [t,y]=ode45(dy,[0,1],y(end)); plot(t,y); Rozwiązaniem analitycznym jest funkcja Rozwiązujemy równanie różniczkowe dla czasu od 0.5 do 0.

7 Atraktor Lorenza Układ trzech równań różniczkowych: s=10; r=28; b=8/3; dy=@(t,y)[s*(y(2)-y(1));... -y(1)*y(3)+r*y(1)-y(2);... y(1)*y(2)-b*y(3)]; [t y]=ode45(dy,[0 100],[0 0.5 1]); plot3(y(:,1),y(:,2),y(:,3)) view(20,20)

8 Model ciężarka Wprowadzamy nowe zmienne Obliczamy pochodne x 1 i x 2 Rozwiązanie w Matlabie: m=5; b=1.5; k=1; F=1; dx=@(t,x)[x(2);... (F-k*x(1)-b*x(2))/m]; [t,x]=ode45(dx,[0 50],[0 0]); plot(t,x(:,1));

9 Model ciężarka m=5; b=1.5; k=1; F=1; dx=@(t,x)[x(2);... (F-k*x(1)-b*x(2))/m]; [t,x]=ode45(dx,[0 50],[0 0]); plot(t,x(:,1)); m=5; b=1.5; k=1; dx=@(t,x)[x(2);... (sila(t)-k*x(1)-b*x(2))/m]; [t,x]=ode45(dx,[0 50],[0 0]); plot(t,x(:,1)); function F=sila(t) if t

10 Pakiet SIMULINK Simulink jest interaktywnym pakietem przeznaczonym do modelowania, symulacji i analizy układów dynamicznych.

11 Sources

12 Sinks

13 Math Operations

14 Model 1 Sine Wave Amplitude:1 Bias: 0 Frequency:2 Phase:0 Ramp Slope:0.1 Start time:0 Initial output:0

15 Model 2 – równanie różniczkowe

16 Model 3 – układ równań różniczkowych

17 Model 4 – równanie różniczkowe II rzędu

18 Model ciężarka Prędkość to całka z przyspieszenia Położenie to całka z prędkości Przyspieszenie jest równe (F-kx-bx’)/m

19 Model wahadła na wózku

20

21

22

23 M=1; m=1; b1=1; b2=4; L=3; F=5; g=9.81; x_0=0; a_0=0; x1_0=0; x1_max=5; a1_0=0; T=50; dt=1e-1; time=0:dt:T; U=[time.',... [ones(101,1);... zeros(100,1);... -1*ones(100,1);... zeros(200,1)]]; sim('model.mdl'); figure('Color',[1 1 1]) subplot(4,2,1) plot(t,U(:,2)) title('Sterowanie') subplot(4,2,3) plot(t,x); title('Polozenie wozka') subplot(4,2,5) plot(t,x1); title('Predkosc wozka') subplot(4,2,7) plot(t,x2); title('Przyspieszenie wozka') subplot(4,2,4) plot(t,a); title('Wychylenie wahadla') subplot(4,2,6) plot(t,a1); title('Predkosc katowa wahadla') subplot(4,2,8) plot(t,a2); title('Przyspieszenie katowe wahadla')

24 Model wahadła na wózku

25 Model ciężarka – SIMULINK/Simscape

26 Prezentacja udostępniona na licencji Creative Commons: Uznanie autorstwa, Na tych samych warunkach 3.0. Pewne prawa zastrzeżone na rzecz autorów. Zezwala się na dowolne wykorzystywanie treści pod warunkiem wskazania autorów jako właścicieli praw do prezentacji oraz zachowania niniejszej informacji licencyjnej tak długo, jak tylko na utwory zależne będzie udzielana taka sama licencja. Tekst licencji dostępny jest na stronie: http://creativecommons.org/licenses/by-sa/3.0/deed.pl