JM

Ta strona jest tylko do odczytu. Możesz wyświetlić źródła tej strony ale nie możesz ich zmienić.
====== 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:

  [estimates, model] = fitcurveQ0(t,q);
  %
  plot(t, q, '*');  hold on
  [sse, FittedCurve] = model(estimates);
  plot(xdata, FittedCurve, 'r')
  % 
  xlabel('t')
  ylabel('q')
  title(['Fitting to function ', func2str(model)]);
  legend('data', ['fit using ', func2str(model)])
  hold off
  
 * 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
  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 = zeros(1,10); err_t = 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(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_t(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,'ro-',1:10,err_t,'bx-');
  xlabel('stopien wielomianu')
  legend('blad w p. treningowych','blad w p. testowych')
JM

Narzędzia użytkownika

Narzędzia witryny

Narzędzia strony
