Spis treści

Zajęcia 1

n = input('Podaj dlugosc ');
disp(n); % tylko wartosc zmiennej
display(n); % nazwa zmiennej i wartosc
 
% alokujemy pamiec tworzac tablice skladajaca sie z samych zer
% n,1 = wymiary
% podajac jeden argument = nxn
% mozna wiecej dodawac wymiarow po przecinku
tab = zeros(n,1);
 
% tab = 1:1:n stworzy wektor od 1 z krokiem 1 do n
 
 %i ma sie zmieniac od 1, z krokiem 1, do n (wykona sie dla n)
for i = 1:1:n
 
    % tablice sa indeksowane od 1
    tab(i)=i;
end
 
% brak { }. Po wszystkim zawsze end
 
% if
% else
% end
 
% wypisujemy tablice
disp(tab);
 
i = 1
while i <= n
    tab(i) = i;
    i = i + 1;
end
 
% jezeli na koncu nie dodamy srednika wypisze na ekranie np efekt
 
% tablice od -pi z krokiem 0.1 do pi
x = -pi:0.1:pi;
y = sin(x);
plot(x, y);
 
figure %otwiera nowe okno
plot(y);
 
% mnozymy y x y
%y = y * y; % nie mozna mnozyc przes siebie dwoch takich samych wektorow
% y = y * y.' % transpozycja [wynikiem jest skalar] wektor poziomy razy pionowy
y = y.' * y % transpozycja .' [wynikiem jest tablica] wektor pionowy razy poziomy
surf(y);
 
% Ax = b
% x =   A^{-1} * b
% x =   A \ b       <=> A^{-1} * b
%       b / A       <=> b * A^{-1}

Zajęcia 2

Trójkąt Pascala

n = input('podaj wymiar: ');
T = cell(n, 1);
 
for i = 1:1:n
    T{i} = ones(1, i);
 
    for j = 2:1:i-1
        T{i}(j) = T{i-1}(j-1) + T{i-1}(j);
    end
end
 
for i = 1:1:n
    disp(T{i});
end

Całka metodą prostokątów

s = 0.1;
x = -1:s:1;
y = x.^2;
 
c = 0;
 
% end = length(y)
%for i = 1 : end - 1
for i = 1 : length(y) - 1
    c = c + y(i) * s;
end
 
display(c)
 
% lub
d = sum(s * y(1: end - 1))

Całka metodą trapezów

s = 0.1;
x = -1:s:1;
y = x.^2;
 
% lub
c = sum( (y(1:end-1) + y(2:end)) * s/2 );
 
display(c)
 
% lub
c = trapz(x, y);
 
display(c)

Całka metodą Monte Carlo

Pp = 2;
Nf = 0;
N = 1000000;
x = 2 * rand(N, 1) - 1;
y = rand(N, 1);
 
for i = 1:1:N
    if y(i) <= x(i)^2
        Nf = Nf + 1;
    end
end
 
c = (Nf / N) * Pp;
 
display(c)
 
% lub
w = x.^2 > y;
%nnz = number of non zero elements
Nf = nnz(w);
 
c = (Nf / N) * Pp;
 
display(c)

Zajęcia 3

Macierze i wektory

w = rand(1,10);
% co trzeci element wektora zamieniamy na odwrotny
w(1:3:end) = -w(1:3:end);
disp(w);
 
% odwroc elementy mniejsze od 0.5
w = rand(1,10);
w(w < 0.5) = -w(w < 0.5);
disp(w);
 
% & - dla tablic, wektorów
% && - dla skalarow
% odwracamy elementy wieksze od 0.3 i mniejsze od 0.6
w = rand(1,10);
w(w > 0.3 & w < 0.6) = -w(w > 0.3 & w < .6)
disp(w);
 
%w = w.*.2; % mnozenie tablicowy, kazdy element mnozymy * 0.2 (0 mozna
%pominac) | 0.2 = .2
% randi - liczba calkowita od 0,1
w = randi(10,1,5);
 
%mean - srednia
% srednia bez wartosci najwiekszej i najmniejszej
%sr = mean(w(w > min(w) & w < max(w)));
%disp(sr);
 
% srednia z zabezpieczeniem odrzucenia np. dwoch maksymalnych elementow
%sr = (sum(w) - min(w) - max(w)) / length(w) - 2;
%disp(sr);
 
disp(w);
 
% wyswietl ile jest liczb parzystych, a ile nieparzystych
a = [ length(w(mod(w,2)==0)), nnz(mod(w,2))];
disp(a);
 
% randn - liczba losowa o rozkladzie normalnym o wartosci oczekiwanej od 0,1
w = randn(1,10);
 
% dodanie wszystkich ujemnych, dodanie wszystkich dodatnich
a = [ sum(w(w<0)), sum(w(w>0))];
disp(a);
% macierz 3x5
n = randi(10,3,5);
disp(n);
 
% dodalismy numery wierszy do macierzy
%n = [(1:3).' n];
%disp(n);
 
% dodalismy numery kolumn do macierzy
%n = [0:5; n];
%n = [0:size(n,2)-1 ; n];
%disp(n);
 
% ; - nowy wiersz w macierzy
% [ a,b ] = macierz pozioma
% [ a ; b ] = macierz pionowa
 
%n = [0:5 ; [(1:3).' n]];
% lub 
%n = [0:size(n,2) ; [(1:size(n,1)).' n]];
%disp(n);
 
% usuwamy kolumne macierzy w ktorej wystepuje wartosc mniejsza od 3
%n(:,i) - i-ta kolumna macierzy
%disp(n);
 
%n(:, min(n) < 3) = [];
% lub
%n(:, ~(min(n) > 3)) = [];
% lub
%n = n(:, min(n) > 3);
 
% usuwamy wiersz macierzy w ktorej wystepuje wartosc mniejse od 6
%n(min(n.').' < 3, :) = [];
%disp(n);
 
% pomijamy drugi argument, szukaj w wierszach(2) lub w kolumnach(1) min(n, [], 2)
%n(min(n, [], 2) < 3, :) = [];
%disp(n);
 
% znajdz wszystkie wartosci minimalne i zamien na NaN;
% n(min(n, [], 2) == min(min(n)), min(n) == min(min(n))) = NaN;
% disp(n);

Sortowanie

% macierz 3x5
n = randi(10,3,5);
disp(n);
 
%[a b] = sort(n)
% a = posortowana macierz
% b = kolejnosc jak on to zmienia
 
% sortuje tylko pierwsza kolumne (i zwraca do a, n jest po staremu)
[ a b ] = sort(n(:, 1));
% lub ~pomija argument [tylko w nowych wersjach]
%[ ~, b ] = sort(n(:, 1));
 
% sortujemy tylko pierwsza kolumne
%n(:,1) = n(b,1);
 
%sortujemy wszystkie kolumny zgodnie z kluczem pierwszej kolumny
n = n(b,:);
 
% sortujemy pierwszy wiersz
[a b ] = sort(n(1,:))
% sortujemy kolumny wg klucza z poprzedniego sortowania
n = n(:, b)
 
disp(n);

Zajęcia 4

Funkcja całki z wielomianu

calka.m

function wy = calka(w, a, b)
    if (nargin == 1)
        wy = polyint(w);
    elseif (nargin == 2)
        wy = polyint(w, a);
    else
        tmp = polyint(w);
        wy = polyval(tmp, b) - polyval(tmp, a);
    end
end

main.m

p = [ 1 0 0 ];
%x^4 + 5x^3 - 2x^2 + 0x + 3
 
%polyint(p); %calka z wielomianu
%polyint(p, k); %warunek couchiego drugi artument, wartosc wyrazu wolnego
 
n = calka(p, 0, 1);
disp(n);

Miejsca zerowe i ekstrema na wykresie

% zamknij wszystkie otwarte wykresy
close all;
 
%p = [ 1 -2 -1 2 4 ];
p = conv(conv(conv([1 1], [1 -1]), [1, 0]), [1 2]);
x = -2:0.01:2;
 
y = polyval(p, x);
dp = polyder(p); % wielomian ktory jest pochodna wielomianu p
 
%plot(x, y);
%hold on;
%plot(r, zeros(size(r)), 'rx');
% mozna tak jak wyzej lub jeden plot na dole
 
%roots - funkcja zwraca pierwiastki wielomianu
r = roots(p);
g = roots(dp);
%plot(r, zeros(size(r)), 'rx', x, y);
 
% polyval(p,g) - wartosc wielomianu p w punktach g
plot(g, polyval(p, g), 'X', r, zeros(size(r)), '*', x, y);

Ekstrema funkcji

% zamknij wszystkie otwarte wykresy
close all;
 
% szukamy elementow ktore nie sa ekstremami, ich druga pochodna musi byc = 0
p = [ 1 0 0 0 ];
x = -2:0.01:2;
 
y = polyval(p, x);
z = roots(p);
 
dp = polyder(p);
rdp = roots(dp);
 
dp2 = polyder(dp);
rdp2 = roots(dp2);
 
for i=1:length(rdp2)
   rdp( rdp == rdp2(i) ) = []; 
end
 
disp(rdp);
disp(rdp2);
 
plot(rdp, polyval(p, rdp), 'X', z, zeros(size(z)), 'r*', x, y);

Praca domowa

p1 = [1 0 0];
p2 = 1;
 
% funkcja kwadratowa
% funkcja stala
 
% napisac program ktory dostaje dwa wieomiany p1 i p2 i zwraca wspolrzedne
% punktow przeciecia
 
% roznica dwoch wielomianow, i liczymy miejsca zerowe tej roznicy
% wielomianow

Zajęcia 5

Kolokwium

Zajęcia 6

GUI

function pushbutton1_Callback(hObject, eventdata, handles)
   wiersze = str2double(get(handles.edit1,'String'));
   kolumny = str2double(get(handles.edit2,'String'));
   A = ones(wiersze, kolumny, 3);
 
   n = str2double(get(handles.edit3,'String'));
   for i = 1:1:n
    x = randi(wiersze);
    y = randi(kolumny);
    rgb = rand(1,3);
    A(x,y,:) = rgb;
   end
   set(handles.figure1,'UserData', A);
   imshow(A);
 
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
 
 
% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
    set(handles.pushbutton2,'UserData', 1);
    A = get(handles.figure1,'UserData');
    B = get(handles.figure1,'UserData');
    n = str2double(get(handles.edit1,'String'));
    m = str2double(get(handles.edit2,'String'));
    todo = true;
    while todo && get(handles.pushbutton2,'UserData')==1
        todo = false;
        for i = 2:1:n-1
           for j = 2:1:m-1
               if(all(B(i,j,:)==1))
                   if(~all(A(i-1,j,:)==1))
                       B(i,j,:) = A(i-1,j,:) ;
                       todo = true;
                   end
                   if(~all(A(i,j-1,:)==1))
                       B(i,j,:) = A(i,j-1,:);
                       todo = true;
                   end
                   if(~all(A(i,j+1,:)==1))
                       B(i,j,:) = A(i,j+1,:);
                       todo = true;
                   end
                   if(~all(A(i+1,j,:)==1))
                       B(i,j,:) = A(i+1,j,:);
                       todo = true;
                   end
               end
           end
        end
        set(handles.figure1,'UserData', A);
        A = B;
        imshow(A);
        pause(0.03);
    end
 
    set(handles.pushbutton2,'UserData', 0);

Zajęcia 7

fun_celu

function y=fun_celu(X,r, a, b)
    N = size(X,1)/2;
    X = reshape(X, N, 2);
    suma = 0;
    for i = 1:N
        if X(i,1)<=a+r || X(i,1)>=b-r || X(i,2)<=a+r || X(i,2)>=b-r
            suma = suma +1000;
        end
        for j = i+1:N
            odl = norm(X(i,:)-X(j,:));
            if odl<2*r
                suma = suma +(2*r - odl);
            end
        end
    end
    y = suma;
end

main

a = 0;
b = 10;
N = 10;
r = 1.5;
k = 0:0.01:2*pi;
X = rand(N,2)*(b-a-2*r)+a+r;
 
for i = 1:N
    x = r * cos(k) + X(i,1);
    y = r * sin(k) + X(i,2);
    fill(x,y,'r','FaceAlpha',0.5);
    hold on;
end
set(gca,'Xlim',[a,b]);
set(gca,'Ylim',[a,b]);
axis equal;
 
X = reshape(X, 2*N, 1);
[x_opt y_opt] = fminsearch(@(x)fun_celu(x,r,a,b), X);
 
X = reshape(x_opt, N, 2);
figure;
for i = 1:N
    x = r * cos(k) + X(i,1);
    y = r * sin(k) + X(i,2);
    fill(x,y,'g','FaceAlpha',0.5);
    hold on;
end
set(gca,'Xlim',[a,b]);
set(gca,'Ylim',[a,b]);
axis equal;

Zajęcia 8

1

Obiekt zawieszony na sprężynie

m = 1;
b = 0.4;
k = 0.5;
 
dx = @(t, x)[x(2);...
    -(k*x(1)+b*x(2))/m];
 
%x(0) = 1
%x'(0) = 0
[t,x] = ode45(dx, 0:0.1:100, [1 0]);
 
plot(t, x(:,1));

2

fun_celu

function ts = fun_celu(m, X)
    dt = 0.05;
    dx = @(t, x)[x(2);-(X(1)*x(1)+X(2)*x(2))/m];
 
    [t,x] = ode45(dx, 0:0.1:100, [1 0]);
    pol = x(:,1);
 
    for i=size(pol):-1:0
        if pol(i) > dt || pol(i) < -dt
            ts = t(i);
            break;
        end
    end
end

main

m = 1;
b = 0.4;
k = 0.5;
 
X = [k;b];
 
[x_opt y_opt] = fmincon(@(X)fun_celu(m, X),[0.5;0.4],[],[],[],[],[0;0],[5;5]);
 
display(y_opt);
display(x_opt);
 
%znalezc czas stabilizacji zeby k,b -> ts -> 0
%do optymalizacji potrzebujemy funkcji celu

Zajęcia 9

Ciężarek z zajęć nr 8

simulink
 
plot(x);

Zawieszenie samochodu

Zajęcia 10

Wariant 1

main

n = 10;
m = 1;
b1 = 0.1;
b2 = 0.5;
F = 1;
 
%a1 = x1;
%a2 = x1';
%a3 = x2;
%a4 = x2';
 
da = @(t, a)[a(2); (sila(t)-b1*a(2)-b2*(a(2)-a(4)))/n;...
             a(4); b2*(a(2)-a(4))/m];
 
%                  [0 10] - od 0 do 10 sekund, dla takiego czasu dostajemy rozwiazanie         
%                          [0 - poczatkowe polozenie wozka, 0 - pocztkowa predkosc wozka, 0 - poczatkowe polozenie klocka, 0 - poczatkowa predkosc klocka]
[t, a] = ode45(da, [0 10], [0 0 0 0]);
%                          [a1 a2 a3 a4]
 
plot(t, a);

sila

function [ f ] = sila( t )
    if t < 5
        f = 1;
    else
        f = -1;
    end
end

Wariant 2

main

M = 10;
m = 1;
b1 = 0.5;
b2 = 0.5;
T = 10;
F = 1;
 
 
 
[F, odl]= fmincon(@FC,2*rand(3,1)-1,[],[],[],[],[-2;-2;-2],[2;2;2]);
 
 
 display(F);
 display(odl);
 
 da = @(t,a)[a(2);(sila(t,F)-b1*a(2)-b2*(a(2)-a(4)))/M ;...
     a(4);b2*(a(2)-a(4))/m];
 
 [t,a] = ode45(da,[0 30], [0 0 0 0]);
 
 display(a);
 plot(t,a);

sila

 
function f=sila(t,F)
    if t<10
        f=F(1);
    elseif t<20 && t>=10
        f=F(2);
    else
        f=F(3);
    end
end

fc

function fc = FC(F)
M = 10;
m = 1;
b1 = 0.5;
b2 = 0.5;
 
da = @(t,a)[a(2);(sila(t,F)-b1*a(2)-b2*(a(2)-a(4)))/M ;...
    a(4);b2*(a(2)-a(4))/m];
[t,a] = ode45(da,[0 30], [0 0 0 0]);
 
fc = (a(end,1)-5)^2+ 10*a(end,2)^2; %ograniczenia na przejechana droge a(end,1)-5)^2 i na to zeby predkosc na koncu byla 0
 
if max(abs(a(:,1)-a(:,1))) >0.1     % ograniczenie zeby klocek nie spadl
    fc = 1000;
end
end

Zajęcia 11

Calka

main

p = [ 1 0 0 ];
 
% liczenie analityczne
% calka               , 1 - przedzial calkowania (od 0 do 1)
c = polyval(polyint(p), 1) - polyint(polyder(p), 0);
disp(c);
 
%liczenie calki numerycznie
% c1 = CALKA[0 1] x^2 dx
c1 = quad(@(x)polyval(p,x), 0, 1);
disp(c1);
 
%liczenie calki numerycznie ktorej funkcja podcalkowa jest w fun.m
p = [ 1 0 ];
L = quad(@(x)fun(p, x), 0, 1);
disp(L);

fun.m

Całka służąca do liczenia długości odcinka

function y  = fun( p, x )
    dp = polyder(p);
    y = sqrt(1 + polyval(dp, x).^2);  
end

Calka 2

x = 0:10;
y = rand(size(x));
arr = zeros(length(x)-1, 2);
 
% zwraca wspolczynniki wielomianu
for i = 1:length(x)-1
    %dopasowanie wielomianu
    p = polyfit(x, y, i);
    arr(i, 1) = quad(@(x)fun(p, x), 0, 1);
    arr(i, 2) = i;
 
    figure(i);
    plot(p);
end
 
bar(arr(:, 1));

Funkcja sklejalna

main

%dopasowanie funkcji sklejanej
% s = spline(x, y);
% y = ppval(s, x);
 
x = 0:10;
y = x.^2;
 
s = spline(x ,y);
sum = 0;
for i = 1:s.pieces
    p = s.coefs(i, :);   
 
    sum = sum + quad(@(x)fun(p, x), 0, 1);
end
 
disp(sum);

fun.m

Całka służąca do liczenia długości odcinka

function y  = fun( p, x )
    dp = polyder(p);
    y = sqrt(1 + polyval(dp, x).^2);  
end

Zajęcia 12

T = 1;
fs = 8192;
t = 0:1/fs:T;
n = length(t);
 
y1 = 2 * sin(2 * pi * t * 200);
 
y2 = sin (2 * pi * t * 600);
 
plot(t(1:200), y1(1:200), t(1:200), y2(1:200));
 
y = y1 + y2;
figure;
plot(t(1:200), y(1:200));
 
% widmo sygnalu
% fft - fast fourier transform
w = abs(fft(y))/n;
% f - wektor czestotoliwosci
f = (0 : n/2) * fs / n;
 
bar (f, w(1:length(f)));
 
% filtrowanie sygnalu
% 1szy arg - rzad filtru
% 2gi  arg - czestotoliwosc odciecia (musi byc z zakresu [0 1]
% 3gi  arg - 'low' - przepuszczaj czestotliwosci ponizej
%            'high' - przepuszczaj czestotliwosci powyzej
[l,m] = butter(12, 2 * 400 / fs, 'low');
y3 = filter(l, m, y);
figure;
plot(t(1:200), y3(1:200), t(1:200), y1(1:200));
 
% rysowanie charakterystyki filtra
[h, theta] = freqz(l, m);
figure;
plot(theta * fs/(2 * pi), abs(h));
 
% dodamy szum do y1
y1 = y1 + randn(size(y1));
figure;
plot(t(1:200), y1(1:200));
 
% widmo sygnalu zaszumionego
s = abs(fft(y1))/n;
% find(s > 0.5)
figure;
plot(f, s(1:length(f)));
 
% probujemy odfiltrowac szum
% dwa filtry (dolnoprzepustowy [wszystko ponizej 220] + gornoprzepustowy [wszystko powyzej 180])
% [l1, m1] = butter(12, 2 * 180 / h, 'high');
% [l2, m2] = butter(12, 2 * 220 / h, 'low');
% y4 = filter(l1, m1, y1);
% y5 = filter(l2, m2, y4);
% figure;
% plot(t(1:200), y5(1:200));
 
% symulacja klawiatury telefonicznej (klawisz 5)
f1 = 1336;
f2 = 770;
y = sin(2 * pi * t * f1) + sin ( 2 * pi * t * f2);
sound(y * fs);

Zajęcia 13

Korelacja

T = 1;
fs = 8192;
t = 0:1/fs:T;
y1 = sin(2 * pi * t * 200);
y2 = sin (2 * pi * t * 500);
y = y1 + y2;
 
%wspolczynnik korelacji
% c = corr(y.', y.');
% disp(c);
% 
% x = 0:0.01:2*pi;
% y1 = sin(x);
% y2 = -sin(x);
% c1 = corr(y1.', y2.');
% disp(c1);
% 
% x = -1:0.01:1;
% y = x.^2;
% c2 = corr(x.', y.');
% disp(c2);
 
% for i=1:200
%    c(i) = corr(y(1:1000).', y(i:i+999).');
% end
% plot((0:199)/fs, c);
Dokładność aproksymacji wielomianowej
% x = -1:0.1:1;
% y = x.^2 + 0.1 * randn(size(x));
% plot(x, y, 'r+');
% p = polyfit(x,y,2);
% 
% x1 = -1:0.01:1;
% y1 = polyval(p, x1);
% hold on;
% plot (x1, y1);
 
x = -1:0.1:1;
y = x.^2 + 0.1 * randn(size(x));
plot(x, y, 'r+');
[p S] = polyfit(x,y,2);
 
x1 = -1:0.01:1;
y1 = polyval(p, x1);
hold on;
plot (x1, y1);
[Y D] = polyconf(p, x, S);
plot(x, Y-D, 'r--');
plot(x, Y+D, 'g--');

Funkcje statystyczne

x = [1 1 1 1 1 1 100];
 
%minimum
v_min = min(x);
v_max = max(x);
%srednia
v_mean = mean(x);
%srednia z x ale odrzucamy 20% pomiarow
%(10 najmniejszych wartosci, 20 najwiekszych)
v_trimmean = trimmean(x, 20);
 
%srednia geometryczna
v_geomean = geomean(x);
 
%srednia harmoniczna
v_harmmean = harmmean(x);
 
%mediana
v_median = median(x);
 
%moda (najczesciej wystepujaca wartosc)
v_mode = mode(x);
 
%zakres
v_range = range(x);
 
%odchylenie standardowe
v_std = std(x);
 
%wariancja
v_var = var(x);
 
% roznica pomiedzy 1szym a 3tym kwantylem x
v_iqr = iqr(x);

Spadek z wysokości

%spadek z wysokosci
close all;
g = 9.81;
h0 = [100 200 300 400 500];
n = length(h0);
m = 20;
t = zeros(m, n);
v = zeros(m, n);
 
for i=1:n
   for j=1:m
      t(j,i) = sqrt(2 * h0(i) / g ) + 0.1 * randn();
      v(j,i) = sqrt(2 * h0(i) * g ) + randn();
   end
end
 
% wykres ramka-wasy
boxplot(t, h0);
figure;
boxplot(v, h0);
 
h1 = [];
for i = 1:n
    h1 = [h1; ones(m,1) * h0(i)];
end
t1 = reshape(t, n * m, 1);
v1 = reshape(v, n * m, 1);
 
figure;
gscatter(v1, t1, h1, 'rbgky', 'xo+sx');
 
t2 = t(:, 1);
v2 = v(:, 1);
figure;
gscatter(v2, t2);
 
% wyznaczenie najbardziej odstajacych, wyznacza kontury wewnatrz ktorych sa
% wszystkie punkty
k = convhull(v2, t2);
hold on;
plot(v2(k), t2(k));
hold off;
 
figure;
 
gscatter(v2, t2);
% wygenerowanie trojkatow
T = delaunay(v2, t2);
hold on;
triplot(T, v2, t2);
hold off;
 
[m1 s1] = normfit(v2);