JM

====== Lab. 1 (2017) ====== - Organizacja ćwiczeń (regulaminy, kartkówki, zadania dodatkowe, zaliczenia, konsultacje, kontakt). - MATLAB - wprowadzenie. - Reprezentacja liczb w formacie zmiennoprzecinkowym - standard IEEE 754. - Dokładność maszynowa - epsilon. ====== Lab. 2 (2017) ====== - Przybliżanie wartości pochodnej skończonym ilorazem różnicowym - badanie zależności błędu przybliżenia od przyrostu argumentu. - Możliwości (ręcznego lub symbolicznego) omijania "niebezpiecznych" operacji arytmetycznych na przykładzie obliczania wartości wyrażenia ''f = sqrt(x^2 +1) - x'' dla x od 1 do 10^20. W miarę możliwości - przykład VPA - variable-precision arithmetic. - Automatyczne różniczkowanie - przykład zastosowania liczb dualnych w języku C++.{{:dydaktyka:aniss:2017:cw2:lab2_stud.cpp|Plik pomocniczy w C++}} - Automatyczne różniczkowanie - konstrukcja grafu obliczania pochodnej w przód. - Różniczkowanie wstecz - przykład konstrukcji grafu. Przydatne funkcje MATLABa (osiągalne także w środowisku ''octave'': -''x = logspace(0,20,21);'' % ''x'' jest wektorem liczb - wyrazów ciągu geometrycznego: 1, 10, 100 ,..., 10^20 -''y = x.^2; '' % ''y'' jest wektorem kwadratów elementów wektora ''x''. -''figure(2), plot(x,y); '' % rysuje wykres nr 2 - funkcję ''y(x)''. Przybornik Symbolic Math Toolbox (z MuPAD) jest dostępny w MATLABie zainstalowanym w lab. 315 (C2) oraz 503 i każdym innym laboratorium UCI (C1). ====== Lab. 3 (2017) ====== - Kartkówka {{:dydaktyka:aniss:2017:pytania1.pdf|Przykładowe pytania}} - Tworzenie modeli ciągłych w SIMULINKU. - Rozwiązywanie zagadnienia początkowego w MATLABie. ====== Lab. 4 (2017) ====== - Model systemu dynamicznego w MATLABie i SIMULINKu. - Model pieca ogrzewającego pokój, z którego ciepło przepływa przez okna na zewnątrz. - Model drapieżnik - ofiara, np. dzięcioły - korniki. Przyrost liczby korników zależy od liczby korników i liczby dzięciołów. Liczba zjadanych korników jest proporcjonalna do iloczynu liczby korników i dzięciołów. Stopa przyrostu liczby dzięciołów zależy od ilości zjadanych korników. - Model balonu unoszonego ciepłym powietrzem. Należy uwzględnić zależność gęstości powietrza od temperatury, proporcjonalną zależność temperatury od wysokości lotu, chłodzenie gazu, dynamikę ruchu w pionie. - Model lodówki w izolowanym pomieszczeniu. Lodówka pobiera energię z sieci energetycznej, która w końcowym efekcie jest zamieniana na ciepło. - Model układu drabinkowego RCRCRC. - Model guza nowotworowego. Komórki zewnętrzne guza mnożą się, wewnętrzne się uspokajają. - Symulacja ruchu ciała o masie m po rzucie ukośnym, przy założeniu oporu powietrza proporcjonalnego do prędkości. - Symulacja ruchu samochodu przy parkowaniu równoległym. - Symulacja ruchu samochodu w pionie przy przejeździe przez próg zwalniający - przy założeniu liniowej/nieliniowej i symetrycznej/niesymetrycznej charakterystyki amortyzatorów. ====== Lab. 5 (2017) ====== Badanie błędu numerycznych metod rozwiązywania zagadnienia początkowego - Budowa diagramu dla równania ''d2x/dt2 + 0.5 dx/dt +3 = 0'' z warunkiem początkowym ''x(0)=1, dx/dt(0)=0''. - Symulacja metodą Eulera z krokiem ''h=0.2'' - Krok symulacji (opcja ''Fixed step size'' w menu okienka modelu: Simulation / Model Configuration Parameters ) ustalamy wpisując do okienka opcji nazwę zmiennej (np. ''h'') utworzoną w przestrzeni roboczej MATLABa (widoczną w SIMULINKu). - Analiza stabilności absolutnej. - Przeniesienie wyników symulacji do przestrzeni roboczej MATLABa (blokiem ''simout''). - Rozwiązanie analityczne ww. równania: - pierwiastki równania charakterystycznego mają część rzeczywistą ''a=-0.25'' i urojoną ''b=sqrt(47)/4'', - postać rozwiązania ''x=exp(a*t) .* (C1*cos(b*t) + C2*sin(b*t))'', - stałe ''C1'' i ''C2'' obliczamy z warunków początkowych i otrzymujemy: ''C1 = 1'', ''C2 = -a/b''. - W Edytorze MATLABa piszemy funkcję (np. ''function x=an_solv(t),'' która oblicza wartość rozwiązania analitycznego w chwilach czasowych ''t''. - Podobnie piszemy funkcję, która rysuje wykres rozwiązania numerycznego i analitycznego zaznaczając obliczone wartości odpowiednio np. kółkami i krzyżykami. - Uzupełniamy ww. funkcję o obliczanie średniego błędu rozwiązania numerycznego - z zastosowaniem 2 lub 3 norm (wyniki umieszczamy w wektorze np. ''err=[e1,e2];'', który jest zwracany przez funkcję np. ''function err=blad(h)''). - Praca w grupach: - wybór metody w opcji okienka modelu: Simulation / Model Configuration Parameters / Solver options = Fixed-step, Solver = jeden z: Euler, Heun, Bogacki-Shampine, Runge-Kutta. - obliczenie błędu średniego dla kroku ''h=1/20, 1/40, 1/80,...,1/2560'' (wyniki proszę wpisywać do tablicy, np. ''wyniki=zeros(9,2); h=0.1; for i=1:9, sim('//nazwa modelu//'); wyniki(i,:)=blad(h); h = h/2; end,'' - analiza wyników: - czy błąd maleje gdy krok maleje? - ile razy (w przybliżeniu)? - obliczenie krotności redukcji błędu przy dwukrotnym zmniejszeniu długości kroku: Wynik (błędy) z i-tego wiersza tablicy ''wyniki'' należy podzielić przez błędy otrzymane w następnym wierszu, czyli z krokiem2 razy krótszym. Można to zrobić tak: ''wyniki(1:8,:)./wyniki(2:9,:)'' //[zapis tab(3:7,:) oznacza tablicę utworzoną z wierszy od 3 do 7 i wszystkich kolumn tablicy tab].// - wyznaczenie relacji między wzrostem dokładności a wzrostem kosztu obliczeniowego. - Porównanie przykładowych wyników otrzymanych dla różnych metod. Tabela średniego błędu kwadratowego: krok Euler Heun B - S RK4 1 3.5805e-04 2.6577e-06 8.7422e-09 7.5712e-11 1/2 1.7769e-04 6.5511e-07 1.0766e-09 4.6300e-12 1/4 8.8513e-05 1.6263e-07 1.3358e-10 2.8627e-13 1/8 4.4173e-05 4.0514e-08 1.6635e-11 1.7785e-14 1/16 2.2066e-05 1.0111e-08 2.0755e-12 1.0896e-15 1/32 1.3962e-05 3.3028e-09 3.4762e-13 1.3585e-16 Tabela krotności zmniejszania błędu krok Euler Heun B - S RK4 1 -> 1/2 2.0150 4.0568 8.1199 16.353 1/2 -> 1/4 2.0075 4.0283 8.0599 16.173 1/4 -> 1/8 2.0038 4.0141 8.0300 16.096 1/8 -> 1/16 2.0019 4.0070 8.0148 16.323 1/16-> 1/32 1.5804 3.0613 5.9707 8.020 - Szacowanie błędu rozwiązania numerycznego w przypadku braku rozwiązania dokładnego - zastosowanie ekstrapolacji Richardsona. Należy: - wybrać metodę i krok oraz przeprowadzić symulację, - sprawdzić, czy krok leży w zakresie regularnego spadku błędu przy zmniejszaniu kroku, - na podstawie uzyskanych rozwiązań przybliżonych oszacować błąd i wyznaczyć rozwiązanie pozbawione głównej części dotychczasowego błędu, - porównać oszacowany błąd z błędem obliczonym za pomocą rozwiązania analitycznego. - Porównanie błędu i kosztu obliczeniowego algorytmu stałokrokowego i adaptacyjnego, w tym algorytmy FSAL. ====== Lab. 6 (2017) ====== Kartkówka - wg zagadnień i przykładowych zadań wysłanych do systemu Dziekanat. Ćwiczenie nr 6: **Wstęp do kalibracji modelu** **Zadanie**: Zbudować model umożliwiający symulację wypływu wody ze zbiornika. Dopuszczalne błędy: - względny błąd symulacji bieżącego poziomu wody ''rel_err_max'' = 1e-3, - błąd bezwzględny ''abs_err_max'' = 2[mm]. **Parametry zbiornika** - Geometria - średnica zbiornika ''D'' = 218 mm, - średnica poziomej rurki wypływowej ''d'' = 9 mm, - oś rurki wypływowej jest na wysokości ''h_r'' = 15 mm nad dnem zbiornika. - Na dnie zbiornika stoi (na większej podstawie) stożek ścięty (nie jest wypełniany wodą): - średnica dolna ''DS'' = 90 mm, - średnica górna ''ds'' = 50 mm, - wysokość ''h_s'' = 55 mm. - Początkowa wysokość lustra wody (liczona od dna zbiornika) ''h0'' = 145 mm. **Model teoretyczny** zgodny z równaniami Bernoulliego. Energia potencjalna zgromadzonej wody jest zamieniana na energię kinetyczną wody wypływającej. Otrzymujemy równanie różniczkowe: ''dh/dt = -k sqrt(h)'', gdzie ''k=(d/D)^2*sqrt(2*g)''. Model w SIMULINKu. **Dane pomiarowe:** Czas ''t'' (w sekundach), poziom wody mierzony od dna naczynia ''h'' (w metrach) t = [0 7.05 14.24 21.95 30.00 38.32 47.08 56.46 66.89 78.72 90.65 104.58 113.64 123.07 135.75]; h = [0.145 0.135 0.125 0.115 0.105 0.095 0.085 0.075 0.065 0.055 0.045 0.035 0.030 0.025 0.020]; **Porównanie wyników symulacji z danymi pomiarowymi**. - wykres danych, - sprawdzenie spełnienia założonych wymagań, - potencjalne źródła rozbieżności: - "driver" symulacji - jak zweryfikować? - model (parametry lub struktura) - analiza założeń upraszczających. **Ilościowa ocena błędu symulacji** Rozwiązanie analityczne dla ''h(0)=h0'': - ''h(t) = (sqrt(h0) - k*t/2)^2'', albo: - ''t(h) = 2*(sqrt(h0) - sqrt(h))/k''. Tę relację zapisujemy w pliku typu ''function'' (skorzystamy na następnym laboratorium). **Plan następnego ćwiczenia** - Kalibracja modelu - znajdowanie optymalnych wartości parametrów modelu. - Aproksymacja danych pomiarowych w przypadku braku opisu teoretycznego. ====== Lab. 7 (2017) ====== Kartkówka - (z zagadnień wymienionych w przesłanych materiałach - 18.12.2017): * aproksymacja danych pomiarowych - aproksymacja a interpolacja, miary błędu aproksymacji, kiedy zagadnienie aproksymacji sprowadza się do zagadnienia liniowego, * numeryczne metody znajdowania minimum funkcji - przykłady metod bezgradientowych oraz gradientowych 1, i 2, rzędu, różnice między metodą Newtona a metodami quasi-newtonowskimi, ogólny schemat metod: minimalizacji wzdłuż współrzędnych (coordinate descent method), najszybszego spadku, Newtona i metody "trust region". **Kalibracja modelu teoretycznego na przykładzie wypływu wody ze zbiornika**. - Wybór funkcji aproksymującej: * ''h(t) = (sqrt(h0) - k*t/2)^2'' * ''t(h) = 2*(sqrt(h0) - sqrt(h))/k'' - Sprowadzenie problemu do zagadnienia liniowego. - Obliczenie wartości przebiegu symulowanego w punktach pomiaru. - Funkcja błędu. - Obliczenie optymalnych wartości parametrów modelu. - Walidacja modelu. % Kalibracja modelu wypływu wody format long e, format compact % Dane pomiarowe: t = [0 7.05 14.24 21.95 30.00 38.32 47.08 56.46 66.89 78.72 90.65 104.58 113.64 123.07 135.75]; h = [0.145 0.135 0.125 0.115 0.105 0.095 0.085 0.075 0.065 0.055 0.045 0.035 0.030 0.025 0.020]-0.015; h0 = 0.13; % początkowy poziom lustra wody. % Wybór danych do kalibracji (ponad stożkiem): ti = t(1:10); hi = h(1:10); % Obliczania optymalnego parametru a(=1/k): ds = sqrt(h0)-sqrt(hi); a = 0.5*ds*ti'/(ds*ds'); % uwaga: w*w' to iloczyn skalarny k = 1/a % Wykres pomiarów ('o') i symulacji ('-') ts = 0:0.1:140; hs = (sqrt(h0)-k*ts/2).^2; plot(t,h,'o',ts,hs), legend('pomiary','symulacja'); **Budowa modelu na podstawie danych pomiarowych - Aproksymacja dyskretna.** * Aproksymacja przebiegu czasowego stanu wody w rzece. * Dane: 69 pomiarów na posterunku w Stróży przepływu wody w Rabie od 17.11.2015 g. 10.00 do 20.11 g. 6.00: q = [2.40 2.56 2.56 2.56 2.72 2.88 3.20 4.04 5.56 6.20 9.30 14.70 ... 18.20 21.80 27.80 30.20 31.40 30.20 29.00 27.80 25.40 24.20 23.00 21.80 ... 21.80 20.60 19.40 19.40 18.20 17.00 17.00 15.90 14.70 14.70 13.50 13.50 ... 12.30 12.30 11.10 11.10 10.50 10.50 10.50 9.90 9.90 9.90 9.30 9.30 ... 9.30 8.70 8.70 8.70 8.70 8.70 8.10 8.10 8.10 8.10 7.34 6.96 ... 6.96 6.96 6.96 6.96 6.58 7.34 6.96 6.96 6.96]; * Aproksymacja opadania fali - od 18. pomiaru, ''q=q(18:end)'' Funkcja szukająca optymalnych wartości parametrów funkcji wykładniczej: function [estimates, model] = fitcurveQ0(xdata, ydata) % Call fminsearch with a random starting point. n = 2; % n - number of parameters, start_point = rand(1, n); model = @expfun; opt = optimset('MaxFunEvals',60000,'MaxIter',50000,'Display','final'); estimates = fminsearch(model, start_point,opt); % % expfun accepts curve parameters as inputs, and outputs sse, % the sum of squares error for A*exp(-lambda*xdata)-ydata, % and the FittedCurve. FMINSEARCH only needs sse, but we want % to plot the FittedCurve at the end. % function [sse, FittedCurve] = expfun(params) A1 = params(1); lambda1 = params(2); FittedCurve = A1 .* exp(-lambda1 * xdata); ErrorVector = FittedCurve - ydata; sse = sum(ErrorVector .^ 2); end end * Wywołanie funkcji aproksymującej: % Dane pomiarowe: q = [2.40 2.56 2.56 2.56 2.72 2.88 3.20 4.04 5.56 6.20 9.30 14.70 ... 18.20 21.80 27.80 30.20 31.40 30.20 29.00 27.80 25.40 24.20 23.00 21.80 ... 21.80 20.60 19.40 19.40 18.20 17.00 17.00 15.90 14.70 14.70 13.50 13.50 ... 12.30 12.30 11.10 11.10 10.50 10.50 10.50 9.90 9.90 9.90 9.30 9.30 ... 9.30 8.70 8.70 8.70 8.70 8.70 8.10 8.10 8.10 8.10 7.34 6.96 ... 6.96 6.96 6.96 6.96 6.58 7.34 6.96 6.96 6.96]; t = 1:52; q=q(18:end); % [estimates, model] = fitcurveQ0(t,q); % aproksymacja nieliniowa [sse, FittedCurve] = model(estimates); % wykres znalezionej aproksymacji plot(t, q, '*',t, FittedCurve, 'ro') % xlabel('t'), ylabel('q'), title(['Fitting to function ', func2str(model)]); legend('data', ['fit using ', func2str(model)]) * Aproksymacja funkcją z większą liczbą parametrów 2. **Aproksymacja wielomianowa - obliczanie współczynników wielomianu za pomocą macierzy Vandermonde'a**: format compact x = -1.5:0.1:1.5; yt = 3.5 - 2*x + 1.5*x.^2 - 2*x.^3 + 2*rand(1,length(x)); y = yt'; figure(1), plot(x,y,'o'), hold on, V1 = vander(x); V = V1(:,end:-1:1); xc = linspace(x(1),x(end),1000); for p = 1:6, A = V(:,1:p+1); [UU,SS,VV] = svd(A); a = inv(A'*A)*A'*y, ya = polyval(a(end:-1:1),x); d = norm(yt-ya), plot(xc,polyval(a(end:-1:1),xc),'b'); pause end hold off 3. **Aproksymacja wielomianowa - kryterium wyboru stopnia wielomianu**: clear, close all, format short, format compact t = 0:0.001:1; y = sin(2*pi*t); % zależność aproksymowana tz = 0:0.05:1; yz = sin(2*pi*tz) + (rand(1,21)-0.5); % wyniki pomiarów % tp = tz(1:2:end); yp = yz(1:2:end); % zbiór treningowy t_test = tz(2:2:end); y_test = yz(2:2:end); % zbiór testowy % figure(1), plot(t,y,'r',tp,yp,'o',t_test,y_test,'x'), hold on, % err_trening = zeros(1,10); err_test = zeros(1,10); % tablice błędów for n = 1:10 p = polyfit(tp,yp,n); % p - współczynniki wielomianu ypa = polyval(p,tp); % wartości funkcji aproksymującej dy = ypa-yp; err_trening(n) = sqrt((dy*dy')/11), % błąd w punktach treningowych % y_test_a = polyval(p,t_test); % wartości f.aproks. w punktach testowych dy_test = y_test_a-y_test; err_test(n) = sqrt((dy_test*dy_test')/10), % błąd w punktach testowych % plot(t,polyval(p,t)), pause, end hold off % figure(2), plot(1:10,err_trening,'ro-',1:10,err_test,'bx-'); xlabel('stopien wielomianu') legend('blad w p. treningowych','blad w p. testowych') ====== Terminy zaliczeń (2017/18) ====== Drugi termin: 31 stycznia 2018. Trzeci termin: 7 lutego 2018. (Godziny - do ustalenia)