background image

METODY NUMERYCZNE - LABORATORIUM 

 

 

PODSTAWY JĘZYKA 

 

MATLAB 

(MATrix LABoratory) 

Scilab

 

(Scientific laboratory) 

 

Gdańsk 1

4.10. 2009 

 
 

 
 

dr inż. Marcin Kujawa 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab

 

Literatura: 

 
Literatura dotycząca MATLAB'a po polsku: 
 

A.Zalewski R.Cegieła: Matlab - Obliczenia numeryczne i ich zastosowania, wyd. Nakom, Poznań 1996. 

B.Mrozek, Z.Mrozek: MATLAB uniwersalne środowisko obliczeń naukowo-technicznych, PLJ 1996. 

J.Brzózka, L.Dorobczyński: Programowanie w Matlab, wyd. Mikom 1998. 

B.Mrozek, Z.Mrozek: MATLAB 5.x, Simulink 2.x, wyd. PLJ 1998. 

Z.Wróbel, R.Koprowski: Przetwarzanie obrazu w programie MATLAB, wyd. Uniw. Śl., K-ce 2001. 

A.Kamińska, B.Pańczyk: Matlab - przykłady i zadania, wyd. Mikom 2002. 

M.Stachurski: Metody numeryczne w programie Matlab, wyd. Mikom 2003. 

W.Regel: Statystyka matematyczna w Matlab, wyd. Mikom 2003. 

W.Regel: Wykresy i obiekty graficzne w MATLAB, wyd. Mikom 2003. 

R.Jankowski, I.Lubowiecka, W.Witkowski: Podstawy programowania w języku Matlab, Gdańsk 2003. 

B.Mrozek, Z.Mrozek: MATLAB i Simulink. Poradnik użytkownika, wyd. Helion 2004. 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

Matlab w Internecie: 

http://www.mathworks.com/ 
 
Lutecki: Przykłady zastosowań MATLABA 
http://www.lutecki.republika.pl/matlab.html 
 
Ł.Przewoźnik: Metody numeryczne i MATLAB 
http://oen.dydaktyka.agh.edu.pl/dydaktyka/matematyka/c_metody_numeryczne/ 
 
R.Buczyniski, R.Kasztelanic: Kurs Matlaba 
http://www.igf.fuw.edu.pl/ZOI/Matlab/index.html 
 
 

 

Zestawienie aktualnych propozycji książkowych dostępnych w księgarniach - porównanie produktów 
 
http://www.nokaut.pl/szukaj/matlab.html 

 
 
 

 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab

 

Literatura dotycząca SciLaba: 
 

A.Brozi: Scilab w przykładach, wyd. Nakom, 2007. 

 

C. T. Lachowicz: Matlab, Scilab, Maxima. Opis i przykłady 
zastosowań
, wyd. Politechniki Opolskiej, 2005. 

 
Scilab w Internecie: 
 
Strona domowa Scilaba: 
http://www.scilab.org/ 
 
Wikipedia Scilaba: 
http://wiki.scilab.org/ 
 
Dokumentacja do Scilaba: 
http://www.scilab.org/product/ 
http://wiki.scilab.org/howto/scilab_documentation_kit 
 
Komendy Scilaba: 
http://www.scilab.org/product/man/ 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

PODSTAWY 
 

1)  Wydawanie poleceń  - komendy można wpisywać z klawiatury - są wykonywane natychmiast po 

naciśnięciu klawisza „Enter”). Przed wykonaniem komendę można edytować przemieszczając kursor 
za pomocą strzałek w prawo i w lewo. Naciśnięcie klawisza „Enter” powoduje wykonanie całego 
polecenia) – niezależnie od aktualnej pozycji kursora w linii polecenia. 

 
 
2)  Skrypty - za pomocą dowolnego edytora tekstowego można utworzyć plik („czysty” plik tekstowy) 

o dowolnej nazwie i rozszerzeniu (tzw. skrypt) zawierający wiele poleceń, które zachowują się jak 
program - są  interpretowane w miarę  wczytywania kolejnych poleceń  z pliku przez Scilab. 
Zwyczajowo plikom zawierającym skrypty nadaje się rozszerzenie „.sce” jednak tak na prawdę jest 
to „nazwa zwyczajowa”, która przydaje się  głównie systemowi Windows (do rozpoznania jakim 
programem należy otworzyć  plik z takim właśnie rozszerzeniem). Podobnie jest z rozszerzeniem 
„.sci”, które zwyczajowo przypisywane jest plikom zawierającym

 

funkcje. 

 
 

3)  Edytor - Scilab wyposażony jest w edytor SciPad, który rozpoznaje składnię poleceń Scilaba i 

odpowiednio koloruje słowa kluczowe. Można w nim pisać programy i funkcje oraz od razu 
przesyłać je do uruchomienia w Scilabie. Uruchamiamy go z menu (pozycja „Editor”). Edytor ten 
rozpoznaje wybrany w Scilabie katalog roboczy i domyślnie właśnie w nim szuka plików do 
otwarcia. 

 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
4)  Historia - Scilab pamięta pewną  liczbę  wcześniej wydanych komend. Naciśnięcie strzałki 

skierowanej w górę  powoduje pojawienie się  ostatnio wydanej komendy, kolejne jej naciśnięcia 
przywołują wcześniejsze komendy. Strzałka w dół też działa w tym kontekście. 

 
 

 

5)  Dopasowanie - Plik „scilab.ini” - początkowo nie istnieje, można go jednak stworzyć w aktualnym 

(w chwili uruchamiania Scilaba) katalogu. Bezpośrednio po zainstalowaniu Scilaba katalogiem tym 
jest \Documents and Settings\nazwa użytkownika. W pliku tym można umieścić wszelkie komendy 
jakie chcielibyśmy aby Scilab wykonywał bezpośredni po uruchomieniu. W szczególności może 
znaleźć się tam polecenie chdir ('ścieżka') wskazujące na katalog, w którym trzymamy pliki Scilaba. 
Jednak definiowanie w tym pliku zmiennych, które chcielibyśmy by były zawsze dostępne, mija się z 
celem bowiem ich „żywotność” sięga nie dalej niż do najbliższego użycia polecenia clear. 

 
 
 
 
 
 
 
 
 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

Po uruchomieniu Scilaba w oknie wprowadzenia ukarze się: 
 
        ___________________________________________         
                       scilab-5.1.1 
 
                 Consortium Scilab (DIGITEO) 
               Copyright (c) 1989-2009 (INRIA) 
               Copyright (c) 1989-2007 (ENPC) 
        ___________________________________________         
   
Startup execution: 
  loading initial environment 
 
-->   gdzie --> jest znakiem zachęty. 

Dla przykładu komenda: 
-->A=[1 1 1;2 4 8;3 9 27] 
 
da na wyjściu: 
A  = 
 
    1.    1.    1.    
    2.    4.    8.    
    3.    9.    27. 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

Działania na wektorach i macierzach 
 
Podstawowe polecenia 

Definicja macierzy

>> A=[1,3,6;2,7,8;0,3,9] 

Odpowiedź: 

A = 
     1     3     6 
     2     7     8 
     0     3     9 

Wielkość macierzy

>> size(A) 

Odpowiedź: 

ans = 
     3     3 

Transpozycja macierzy

>> A' 

Odpowiedź: 

ans = 
     1     2     0 
     3     7     3 
     6     8     9 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
 
Odwołanie do pojedynczych kolumn czy wierszy macierzy

>> A(:,3) 

Odpowiedź: 

ans = 
     6 
     8 
     9   

 

 

 

(trzecia kolumna macierzy A) 

 
>> A(1,:) 

Odpowiedź: 

ans = 
     1     3     6 

 

 

(pierwszy wiersz macierzy A)  

 
Działania na wybranych kolumnach i wierszach 
np. dodawanie: 

>> A(1,:)+A(3,:) 

Odpowiedź: 

ans = 
     1     6    15 

 
 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
Dodawanie macierzy

>> B=[3,4,5;6,7,2;8,1,0]; 
>> C=A+B 

Odpowiedź: 

C = 
     4     7    11 
     8    14    10 
     8     4     9 

Odejmowanie macierzy

>> C=A-B 

Odpowiedź: 

C = 
    -2    -1     1 
    -4     0     6 
    -8     2     9 

Mnożenie macierzy

>> C=A*B 

Odpowiedź: 

C = 
    69    31    11 
   112    65    24 
    90    30     6 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

Wspólne funkcje MATLABA i Scilaba: http://www.scilab.org/product/dic-mat-sci/M2SCI_doc.htm 

 

Funkcje wykorzystywane przy działaniach macierzowych 

 

Symbol Wyjaśnienie 
inv (wspólna) 

obliczanie macierzy odwrotnej 

det (wspólna) 

obliczanie wyznacznika macierzy 

rank (wspólna) 

obliczanie rzędu macierzy 

(liczba niezależnych wierszy i 

kolumn)

 

cond (wspólna) 

wskaźnik uwarunkowania macierzy 

(oszacowanie z jaką 

(maksymalnie) dokładnością (do ilu miejsc po przecinku) możemy 
podać wynik)

 

eye(n,n) (wspólna) 

definicja macierzy jednostkowej o wymiarach n na n 

trace (wspólna) 

obliczanie śladu macierzy 

(suma elementów na głównej 

przekątnej macierzy)

 

zeros(n,m) (wspólna) 

definicja macierzy wypełnionej zerami 

expm (wspólna) 

potęgowanie macierzy 

eig (mtlb_eig lub spec) 

 

obliczanie wartości własnych macierzy 

lu (wspólna) 

dekompozycja macierzy na macierz trójkątną górną i 
dolną 

chol (wspólna) 

rozkład Cholewskiego 

qr (wspólna) 

Dekompozycja macierzy na ortogonalną i trójkątną 
górną 

\ (wspólna) 

znacznik używany do rozwiązania liniowych równań 
algebraicznych 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
 
Funkcje wykorzystywane przy analizie danych 
 
Symbol Wyjaśnienie 
min (max) (wspólna) 

zwraca najmniejszy (największy) element wektora 
lub gdy jest to macierz najmniejszy (największy) 
element każdej kolumny 

sum (wspólna) 

zwraca sumę elementów wektora lub gdy jest to 
macierz sumy elementów każdej kolumny 

std (mtlb_std) 

Zwraca odchylenie standardowe wektora lub gdy jest 
to macierz odchylenia standardowe elementów 
każdej kolumny 

sort (wspólna) 

sortuje wektor w porządku wartości rosnących lub 
gdy jest to macierz sortuje poszczególne elementy 
każdej kolumny 

mean (mtlb_mean) 

zwraca średnia arytmetyczną elementów wektora lub 
gdy jest to macierz średnie arytmetyczne elementów 
każdej kolumny 

 
 
 
 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
Funkcje wykorzystywane przy analizie równań wyższych rzędów 
 
Symbol Wyjaśnienie 
poly oblicza 

współczynnik wielomianu 

charakterystycznego 

roots (wspólna) 

zwraca pierwiastki równania 

polyval zwraca 

wartość wielomianu 

polyfit aproksymacja 

wielomianem 

 
Funkcje wykorzystywane przy analizie nieliniowych równań algebraicznych 
 
Symbol Wyjaśnienie 
fmin szuka 

najmniejszej 

wartości funkcji 

fzero 

szuka miejsca zerowego funkcji 

 
PRZYKŁADY ZASTOSOWANIA WYBRANYCH FUNKCJI 
 
>> x=fzero('sin',10); % szuka miejsca zerowego funkcji sinus w okolicach 10 
>> m=fminbnd('sin',10,11); % szuka najmniejszej wartości funkcji sinus w przedziale (10,11) 
>> p=[3,-2,4,1]; % definicja wielomianu p=3*x^3-2*x^2+4*x+1 
>> r=roots(p); % zwraca pierwiastki równania 
>> w=polyval(p,2); % zwraca wartość wielomianu 3*2^3-2*2^2+4*2+1 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
Funkcje matematyczne 
 
Symbol Wyjaśnienie 
sin(x) (wspólna) 

sinus 

cos(x) (wspólna) 

cosinus 

tan(x) (wspólna) 

tangens 

asin(x) (wspólna) 

arcus sinus 

acos(x) (wspólna) 

arcus cosinus 

atan(x) (wspólna) 

arcus tangens 

sinh(x) (wspólna) 

sinus hiperboliczny 

cosh(x) (wspólna) 

cosinus hiperboliczny 

tanh(x) (wspólna) 

tangens hiperboliczny 

asinh(x) (wspólna) 

arcus sinus hiperboliczny 

acosh(x) (wspólna) 

arcus cosinus hiperboliczny 

atanh(x) (wspólna) 

arcus tangens hiperboliczny 

sqrt(x) (wspólna) 

pierwiastek kwadratowy 

exp(x) (wspólna) 

e do x 

log(x) (wspólna) 

logarytm naturalny 

log2(x) (wspólna) 

logarytm przy podstawie z 2 

log10(x) (wspólna) 

logarytm przy podstawie z 10 

 
 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

 
 
Funkcje związane z obliczeniami w dziedzinie liczb zespolonych 
 
Symbol Wyjaśnienie 
abs(x) (wspólna) 

macierz modułów elementów macierzy x 

angle(x) (wspólna) 

macierz argumentów macierzy x 

real(x) (wspólna) 

macierz części rzeczywistych elementów macierzy x 

imag(x) (wspólna) 

macierz części urojonych elementów macierzy x 

conj(x) (wspólna) 

macierz o elementach sprzężonych z elementami 
macierzy x 

 
Funkcje dodatkowe 
 
Symbol Wyjaśnienie 
round(x) (wspólna) 

zaokrągla elementy macierzy x do najbliższej liczby 
całkowitej 

rem(x,y) (wspólna) 

oblicza resztę z dzielenia odpowiadających sobie 
elementów macierzy x i y 

gcd(a,b)  

oblicza największy wspólny dzielnik liczb a i b 

lcm(a,b) oblicza 

najmniejszą wspólną wielokrotną liczb a i b 

 
 

background image

METODY NUMERYCZNE - LABORATORIUM 

 

MATLAB & Scilab 

ĆWICZENIA 
 

•  Zdefiniować macierz A wymiaru n następująco (sprawdź szczegóły dotyczące funkcji diag za pomocą 

Help): 

A  = 
   - 2.    0.    0.    0.    0.   
     0.  - 1.    0.    0.    0.   
     0.    0.    0.    0.    0.   
     0.    0.    0.    1.    0.   
     0.    0.    0.    0.    2.  
 

•  Niech A będzie macierzą kwadratową A=[1,2,3;4,5,6;7,8,9]; czym będzie diag(diag(A)) ? 
 
•  Co stanie się gdy wydamy polecenie clear i czym się różni ono od polecania clear A ? (sprawdź 

szczegóły dotyczące funkcji clear za pomocą Help). 

 

•  Zdefiniuj macierz z jedynkami na głównej przekątnej o wymiarach 3 na 3. 

 

•  Narysuj wykres funkcji f(x)=x

3

+2x

2

-3 w przedziale przy x od -4 do 4. 

 

•  Do czego służy funkcja clf() ? 

 
 

background image

 

METODY NUMERYCZNE 

 

Temat 1 

INTERPOLACJA I APROKSYMACJA 

 

 
Ćwiczenia opracowano na podstawie: 

1)  Zalwewski A., Cegieła R. „MATLAB-obliczenia numeryczne i ich 

zastosowanie” Nakom, Poznań 1997. 

2)  Jankowski R., Lubowiecka I., Witkowski W. „Podstawy programowania w 

języku Matlab” Politechnika Gdańska, WIL, 2002 

3)  http://www.mathworks.com/access/helpdesk.html 
4)  Brozi A. Scilab w przykładach” Nakom, Poznań 2007 

 
 
 
 

background image

MATLAB & Scilab

 

INTERPOLACJA 

DEF.: 

Interpolacja to wyszukiwanie analitycznego wzoru funkcji na podstawie 
częściowej tablicy jej wartości, zwanych węzłami interpolacji. 
 

INTREPOLACJA FUNKCJI JEDNEJ ZMIENNEJ 

Standardowe procedury MATLAB-a realizują interpolację za pomocą 
następujących metod: 

- interpolacja wielomianami 
- interpolacja za pomocą funkcji sklejanych 

 
1) interpolacja wielomianami pierwszego i trzeciego stopnia 
Powszechnie znanym wielomianem jest wielomian interpolacyjny Lagrange’e 
postaci: 

1

1

1

1

1

1

1

(

)...(

)(

)...(

)

( )

(

)...(

)(

)...(

)

N

i

i

n

N

i

i

i

i

i

i

i

i

n

x x

x x

x x

x x

W x

y

x

x

x

x

x

x

x

x

+

=

+

=

 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

Podstawową wadą takich interpolacji jest skłonność do występowania silnych 
oscylacji między węzłami interpolacji dla wielomianów interpolacyjnych 
wysokich stopni. 
 
2) interpolacja za pomocą funkcji sklejanych 
Metoda ta sprowadza się do podziału przedziału, w którym dokonywana jest 
interpolacja, na kilka podprzedziałów oraz dokonywania interpolacji na 
kolejnych obszarach interpolowanej funkcji za pomocą wielomianów niskiego 
rzędu przy jednoczesnej gwarancji „gładkiego” przejścia między kolejnymi 
przedziałami. 
 
 

W dziedzinie obliczeń numerycznych interpolacja spełnia zwykle role 

pomocniczą 

np. dla konstrukcji metod całkowania i różniczkowania numerycznego 

lub w metodach optymalizacji. 

 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB 

 

FUNKCJA BIBLIOTECZNA 
 
Interpolacja funkcji jednej zmiennej w punktach określonych wektorem x

i

 

 

yi=interp1(x,y,xi,’metoda’)

 

gdzie: 
x,y - węzły interpolacji (elementy wektora x muszą tworzyć ciąg rosnący) 
’metoda’ - opcjonalny parametr pozwalający na wybór metody interpolacji 
 

’linear’ - interpolacja funkcja łamaną, 
’cubic’ - interpolacja wielomianami trzeciego stopnia (odległości 
między elementami x muszą być równe). 
’spline’ - interpolacja funkcjami sklejanymi trzeciego stopnia, 
 

 
 
 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB 

 

PRZYKŁAD  
 
Interpolacja liniowa 
x=0:20; y=sin(x)+sin(2*x); 
xi=0:.2:20; 
figure(1) 
yi=interp1(x,y,xi,’linear’); 
plot(x,y,’o’,xi,yi,xi,sin(xi)+sin(2*xi)) 
xlabel(‘x’) 
ylabel(‘y’) 
 

interpolacja funkcji y=sin(x)+sin(2x) 

łamaną, uzyskaną wskutek 

uwzględnienia 21 węzłów 

rozmieszczonych równomiernie w 

przedziale <0,20> 

 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB 

 

Interpolacja typu - cubic 
figure(2) 
yi=interp1(x,y,xi,’cubic’); 
plot(x,y,’o’,xi,yi,xi,sin(xi)+sin(2*xi)) 
xlabel(‘x’) 
ylabel(‘y’) 
 
 
 

interpolacja funkcji y=sin(x)+sin(2x) 

wielomianami trzeciego stopnia, przy 

uwzględnieniu 21 węzłów 

rozmieszczonych równomiernie w 

przedziale <0,20> 

 
 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

 

MATLAB 

 

Interpolacja typu - spline 
figure(3) 
yi=interp1(x,y,xi,’spline’); 
plot(x,y,’o’,xi,yi,xi,sin(xi)+sin(2*xi)) 
xlabel(‘x’) 
ylabel(‘y’) 
 
 
 
 

interpolacja funkcji y=sin(x)+sin(2x) 

funkcjami sklejanymi, przy 

uwzględnieniu 21 węzłów 

rozmieszczonych równomiernie w 

przedziale <0,20> 

 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab 

PRZYKŁAD  

 

Interpolacja liniowa 

 
x=0:20; y=sin(x)+sin(2*x); 

//węzły interpolacji 

f1=scf(1); 

//wykres 1

 

plot2d(x,y,style=-9); 
xi=0:.2:20; 
plot2d(xi,sin(xi)+sin(2*xi),style=5); 
yi=interpln([x;y],xi); 
plot2d(xi,yi,style=2); 
xtitle

 

('Interpolacja liniowa', 'x', 'y') 

 
 
 
 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

 

Scilab 

 
 

Interpolacja typu – spline 

 
f2=scf(2);

 //wykres 2

 

plot2d(x,y,style=-9); 
plot2d(xi,sin(xi)+sin(2*xi),style=5); 
d=splin(x,y); 
yi=interp(xi,x,y,d); 
plot2d(xi,yi,style=2); 
xtitle

 

('Interpolacja typu - spline ', 'x', 'y') 

 
 
 
 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

 

MATLAB

 

INTREPOLACJA FUNKCJI DWÓCH ZMIENNYCH 

 
FUNKCJA BIBLIOTECZNA 
 
Interpolacja funkcji dwóch zmiennych w punktach określonych macierzą [X, Y] 

 

yi=interp2(x,y,z,xi,yi,zi,’metoda’)

 

PRZYKŁAD 
 
[X,Y] = meshgrid(-3:.25:3); 
Z = peaks(X,Y);

%dowolna funkcja dwóch zmiennych 

[XI,YI] = meshgrid(-3:.125:3); 
ZI = interp2(X,Y,Z,XI,YI,’linear’); 
mesh(X,Y,Z), hold, mesh(XI,YI,ZI+25) 
hold off 
axis([-3 3 -3 3 -5 30]) 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

10 

background image

 

MATLAB

 

 
 

 

 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

11 

background image

MATLAB

 

 

APROKSYMACJA 

DEF.: 

Aproksymacja to dowolne przybliżenie pewnej wartości lub pewnego zbioru 
wartości bądź ciągu wartości (np. funkcji) według zadanej dokładności. 
 

Aproksymacja wielomianem 

 
FUNKCJA BIBLIOTECZNA 
 

f=polyfit (x,y,r)

 

gdzie: 
x,y - wektor danych 
r - stopień wielomianu 
Funkcja ta dla wektorów danych x i y generuje współczynniki wielomianu 
stopnia  r  przybliżającego najlepiej w sensie średniokwadratowym zależność 
miedzy serią danych x i y. 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

12 

background image

MATLAB

 

PRZYKŁAD 
x=[6.1 8 10.3 11.3 12.3 13.2 14.2 16.1];

%zbiór rzędnych

 

y=[61.3 50 46.4 46.4 40.8 31.6 29.9 22];

%zbiór wartości

 

n=1;

%stopień wielomianu (dla n=1 funkcja liniowa) 

 
f=polyfit(x,y,n);

%tablica współczynników wielomianu

 

 
p=min(x)

%początek zbioru rzędnych xn

 

k=max(x)

%koniec zbioru rzędnych xn

 

dx=(k-p)/50; 

%krok podziału zbioru xn

 

xn=(p:dx:k); 

%zbiór rzędnych xn

 

for

 i=1:length(xn) 

    yn(i)=polyval(f,xn(i));

%zbiór wartości funkcji

 

end 

figure(1) 
plot(x,y,

'o'

,xn,yn) 

 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

13 

background image

MATLAB

 

n=3;

%stopień wielomianu

 

f=polyfit(x,y,n); 

for

 i=1:length(xn) 

    yn(i)=polyval(f,xn(i));

%zbiór wartości wielomianu

 

end 

figure(2) 
plot(x,y,

'o'

,xn,yn) 

 
 
 
 
 
 
 
 
 
 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

14 

background image

Scilab 

PRZYKŁAD 

 
 

//------------------------------------------ 
//Aproksymacja liniowa 
//f=a0+a1*x - równanie funkcji aproksymującej 
//a0=ymean-a1*xmean - stała a0 
//a1=(n*c1-c2)/(n*c3-c4) - stała a1 
 
//------------------------------------------ 
//Dane 
//------------------------------------------ 

x=[6.1 8 10.3 11.3 12.3 13.2 14.2 16.1]; 
y=[61.3 50 46.4 46.4 40.8 31.6 29.9 22]; 
 
 
 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

15 

background image

Scilab 

//Start programu 
//------------------------------------------ 

n=size(x,'*'); 

//liczba prób

 

xmean=mean(x); 

//średnia z x 

ymean=mean(y); 

//średnia z y 

//czynnik c1 

for i=1:n do 
  c(i)=x(i)*y(i); 
end 
c1=sum(c); 

//czynnik c2 

c2=sum(x)*sum(y); 

//czynnik c3 

for i=1:n do 
  d(i)=x(i)^2; 
end 
c3=sum(d); 

 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

16 

background image

Scilab 

17 

Metody Numeryczne 

 1. Interpolacja i aproksymacja  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

//czynnik c4 

c4=(sum(x))^2; 
a0=ymean-a1*xmean; 
a1=(n*c1-c2)/(n*c3-c4); 

//------------------------------------------ 
//Funkcja aproksymująca 
//------------------------------------------ 

p=min(x) 

//początek zbioru rzędnych xn

 

k=max(x) 

//koniec zbioru rzędnych xn 

dx=(k-p)/50; 

//krok podziału zbioru xn 

xn=(p:dx:k); 

//zbiór rzędnych xn

 

plot2d(x,y,style=-9); 
plot2d(xn,a0+a1*xn,style=5); 

//------------------------------------------ 
//Koniec programu 

background image

 

METODY NUMERYCZNE 

 

Temat 2 

WIELOMIANY, MIEJSCA ZEROWE 

 

 
Ćwiczenia opracowano na podstawie: 

1)  Zalwewski A., Cegieła R. „MATLAB-obliczenia numeryczne i ich 

zastosowanie” Nakom, Poznań 1997. 

2)  Jankowski R., Lubowiecka I., Witkowski W. „Podstawy programowania w 

języku Matlab” Politechnika Gdańska, WIL, 2002 

3)  http://www.mathworks.com/access/helpdesk.html 
4)  Brozi A. Scilab w przykładach” Nakom, Poznań 2007 

 
 
 
 

background image

MATLAB

 

 
 
 

MATLAB posiada kilka funkcji bibliotecznych znajdujących zera funkcji jednej 
zmiennej oraz pozwalających wygodnie operować na wielomianach. 
 
 
FUNKCJA BIBLIOTECZNA 
fzero - 

oblicza miejsce zerowe funkcji o podanej nazwie zaczynając obliczenia z 

punktu startowego x0, który winien być początkowym przybliżeniem wartości 
szukanego miejsca zerowego. 
 

z=fzero(‘F’,x0,esp) 

gdzie: 
esp

 - żądana dokładność, 

 
 
 

Metody Numeryczne 

 2.Wielomiany, miejsca zerowe  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

PRZYKŁAD

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 MATLAB

 

Tworzymy skrypt funkcji w pliku o nazwie fun.m 

function

 [Y]=fun(X) 

%funkcja 

Y=sin(X)-X*cos(X) 

% koniec fun.m 

Obliczamy wartość miejsca zerowego w oknie komend 
 

z=fzero('fun',4.71239,10^-3) 
z=

 

4.4934 

FUNKCJA BIBLIOTECZNA 
poly

 - wyznacza współczynniki a(1), a(2),…, a(n+1) wielomianu 

1

( , )

(1)

(2)

(

1)

n

n

W x a

a

x

a

x

a n

=

+

+

+

 

o pierwiastkach podanych w wektorze r. 

a=poly(r) 

 
Współczynniki są uszeregowane według malejących potęg zmiennej x

Metody Numeryczne 

 2.Wielomiany, miejsca zerowe  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

PRZYKŁAD

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 MATLAB

 

W oknie komend wpisujemy 
r=[1 2 3 4 5] 
a=poly(r) 
 
a=[1   -15    85  -225   274  -120] 
 
FUNKCJA BIBLIOTECZNA 
roots

 - funkcja oblicza pierwiastki wielomianu 

 o postaci, jak podano w 

opisie funkcji poly 

( , )

W x a

 

r=roots(a) 

PRZYKŁAD 
W oknie komend wpisujemy 
a=[1   -15    85  -225   274  -120] 
r=roots(a) 
 
r=[5 4 3 2 1] 

Metody Numeryczne 

 2.Wielomiany, miejsca zerowe  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Metody Numeryczne 

 2.Wielomiany, miejsca zerowe  

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

Scilab

 

 
inv_coeff() 

- definicja wielomianu 

a=[2 3 4 5]; 
p=inv_coeff(a) 
      

           2      3 

    2 + 3x + 4x + 5x 
 

degree() 

- funkcja określa stopień wielomianu 

c=degree(p) 
c=3 

 
roots

 - funkcja oblicza pierwiastki wielomianu 

 
a=[1   -15    85  -225   274  -120]; 
r=roots(a)’ 
 
r=1 2 3 4 5 

background image

 

METODY NUMERYCZNE 

 

Temat 3 

CAŁKOWANIE NUMERYCZNE 

 

 
Ćwiczenia opracowano na podstawie: 

1)  Zalwewski A., Cegieła R. „MATLAB-obliczenia numeryczne i ich 

zastosowanie” Nakom, Poznań 1997. 

2)  Jankowski R., Lubowiecka I., Witkowski W. „Podstawy programowania w 

języku Matlab” Politechnika Gdańska, WIL, 2002 

3)  http://www.mathworks.com/access/helpdesk.html 
4)  Brozi A. Scilab w przykładach” Nakom, Poznań 2007 

 

 
 

 

background image

MATLAB & Scilab

 

CAŁKOWANIE NUMERYCZNE 

 
Metody całkowania numerycznego, zwane kwadraturami sprowadzają się do 
przybliżenia funkcji podcałkowej  f lub jej kolejnych fragmentów na danym 
przedziale <a, b> za pomocą innej funkcji, dla której wartości całki jest 
określona analitycznie. 
 
Stosuje się w trakcie analizy kwadratury: 
- złożone, 
- adaptacyjne. 
 
Idea kwadratury złożonej polega na podziale przedziału <a, b> na pewną liczbę 
podprzedziałów (np. ze względu na dużą zmienność funkcji w przedziale 
całkowania) i obliczenie wartości całki na każdym z tych przedziałów, a 
następnie zsumowanie tych wartości. 
 
 

Metody Numeryczne 

 3. Całkowanie numeryczne 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

 
Techniki adaptacyjne
 polegają na wstępnym podziale przedziału całkowania na 
dwie części równej długości i obliczenie wartości całek w każdym z nich za 
pomocą kwadratury. Przedział, w którym nie osiągnięto wystarczającej 
dokładności obliczeń (porównuje się wartość funkcji w danym przedziale z 
wartością funkcji po podziale na dwie części) ponownie jest dzielony na dwie 
części i tak do czasu aż osiągniemy zadana zbieżność. 
 
 
 
 
 
 
 
 
 
 

 

Metody Numeryczne 

 3. Całkowanie numeryczne 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB

 

 
Standardowe procedury MATLAB-a pozwalające na numeryczne obliczenie 
wartości całki oznaczonej funkcji jednej zmiennej to: 
- quad - kwadratura adaptacyjna oparta o interpolacje wielomianem 2 stopnia, 
- quadl - kwadratura adaptacyjna oparta o kwadraturę Lobatto. 
 
FUNKCJE BIBLIOTECZNE 
 

Q1=quad(f,a,b,tol) 

Ql=quadl(f,a,b,tol) 

gdzie: 
f - łańcuch zawierający nazwę funkcji umieszczonej w osobnym skrypcie, 
a, b - przedział całkowania, 
tol - wymagana tolerancja względna (domyślnie 10^-3 - przy braku parametru). 
 
 
 

 

Metody Numeryczne 

 3. Całkowanie numeryczne 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Metody Numeryczne 

 3. Całkowanie numeryczne 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

MATLAB

 

PRZYKŁAD - Całkowanie numeryczne 

2

0

sin( )

x dx

π

 

 

Tworzymy skrypt funkcji podcałkowej 
 
function [Y]=quadfun(X) 
%funkcja podcałkowa 
Y=sin(X.*X); 
% koniec quadfun.m 
 
Obliczamy wartość całki 
 

Q1=quad(‘quadfun’,0,pi)=

 

0.7727 

lub 

Ql=quadl(‘quadfun’,0,pi)=

 

0.7727 

background image

Metody Numeryczne 

 3. Całkowanie numeryczne 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

Scilab 

 
 
PRZYKŁAD 
- Całkowanie numeryczne 
 

Kwadratura standardowa 
x0=0;x1=%pi; 
X=integrate('sin(x.*x)','x',x0,x1) 
X=

 

0.7726517 

 
Kwadratura adaptacyjna oparta o metodę trapezów (przeznaczona do danych 
eksperymentalnych) 
x=0:0.1:%pi; 

- gęstość podziału

 

X=inttrap(x,sin(x.*x)) 
X=

 

0.7803791 

 

Należy zbadać zbieżność rozwiązania w zależności od gęstości podziału! 

background image

 

METODY NUMERYCZNE 

 

Temat 4 

UKŁADY RÓWNAŃ LINIOWYCH 

I DEKOMPOZYCJA MACIERZY 

 

 
Ćwiczenia opracowano na podstawie: 

1)  Zalwewski A., Cegieła R. „MATLAB-obliczenia numeryczne i ich 

zastosowanie” Nakom, Poznań 1997. 

2)  Jankowski R., Lubowiecka I., Witkowski W. „Podstawy programowania w 

języku Matlab” Politechnika Gdańska, WIL, 2002 

3)  http://www.mathworks.com/access/helpdesk.html 
4)  Brzóska J., Dorobczyński L. „Programowanie w MATLABIE” „MIKOM” 

Warszawa 1998 

5)  Brozi A. Scilab w przykładach” Nakom, Poznań 2007 

background image

WŁASNOŚCI MACIERZY

  

 

 

 

 

 

 

 

 

 

MATLAB & Scilab

 

Rząd macierzy – największy stopień uzyskanej z niej podmacierzy kwadratowej 
nieosobliwej tzn. takiej której wyznacznik jest różny od zera. 
Do obliczenia rzędu macierzy służy funkcja rank. 
 

r=rank(A) 

 
Według twierdzenia Cronecera –Capelliego układ równań liniowych Ax=b 
można jednoznacznie rozwiązać gdy rząd macierzy głównej jest równy rzędowi 
macierzy dołączonej 
 
PRZYKŁAD 
A=[1 2 -1; 2 1 -2; 1 2 2]; b=[1 2 3]’; 
rA=rank(A)=3 
rAb=rank([A b])=3 
 
 
 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

MATLAB & Scilab

 

Wyznacznik macierzy (wykorzystywane we wzorach Cramera) 
 

1

2

1

2

det( )

sgn( )

...

n

p

p

np

p P

A

A

p a

a

a

=

=

⋅ ⋅

 

A - macierz kwadratowa n x n 
P
 – zbiór wszystkich permutacji liczb (możliwych uszeregowań tych liczb) 
od 1 do 
 
Dla każdej permutacji p obliczamy następujący iloczyn: 
 

1

2

1

2

sgn( )

...

n

p

p

np

p a

a

a

⋅ ⋅

 

 

1,  gdy   jest permutacją patrzystą

sgn( )

1,  gdy   jest permutacją niepatrzystą

p

p

p

= ⎨

 

 

i sumujemy dla kolejnych permutacji. 
 

background image

MATLAB & Scilab

 

 
Obliczanie wyznacznika z definicji jest niezmiernie czasochłonne, gdyż wymaga 
wyznaczenia n! permutacji, dlatego w praktyce nie wykonuje się go w ten 
sposób. 
Metoda obliczania wyznacznika polega na sprowadzeniu macierzy do postaci, 
w której wyznacznik jest łatwo obliczyć np. trójkątnej, gdzie wyznacznikiem 
macierzy jest iloczyn elementów na przekątnej. 
Do obliczenia wyznacznika macierzy służy funkcja det. 
 

d=det(A) 

 
PRZYKŁAD 
A=[1 2 10 -2; 3 -3 6 9; 2 4 -1 8; 4 5 6 10]; 
d=det(A) 
d=-351 
 
 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

UKŁADY RÓWNAŃ LINIOWYCH I DEKOMPOZYCJA MACIERZY 
W wielu sytuacjach obliczeniowych (wyznaczniki i rozwiązywanie układów 
równań) korzystne jest przedstawienie danej macierzy A w postaci iloczynu 
macierzy trójkątnych dolnej i górnej: A=LU 
Funkcja lu znajduje macierze L i U takie, że L*U=A 
Procedura 1 

 

     [L,U]=lu(A) 

PRZYKŁAD 
A=[1 2 1;0 4 2; 6 1 1]; 
[L,U]=lu(A) 
L= 
    0.1667    0.4583    1.0000 
             0    1.0000             0 
    1.0000              0            0 

macierz nie jest macierzą trójkątną dolną 

U= 
    6.0000    1.0000    1.0000 
             0    4.0000    2.0000 
             0              0   -0.0833 

dla macierzy zachodzi zależność A=LU 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

Procedura 

2      

 

[L,U,P]=lu(A) 

gdzie zachodzi własność LU=PA, a P to macierz permutacji (nieosobliwa 
zbudowana z zer i jedynek; P

-1

=P) 

[L,U,P]=lu(A) 
L = 
    1.0000              0             0 
             0     1.0000             0 
    0.1667     0.4583    1.0000 
U = 
    6.0000    1.0000    1.0000 
             0    4.0000    2.0000 
             0             0   -0.0833 
P = 
     0     0     1 
     0     1     0 
     1     0     0 
 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

ZASTOSOWANIE ROZKŁADU LU 
UKŁADY RÓWNAŃ LINIOWYCH 
Rozkład LU ułatwia rozwiązanie różnych zadań. Najprostszym jest rozwiązanie 
układów równań liniowych postaci:   Ax=b 
 
Jeżeli znamy rozkład macierzy: 

A=PLU, 

 

równanie wyjściowe wówczas ma postać: 
 

PLUx=b lub LUx=Pb    

(P

-1

=P). 

 

Aby wyznaczyć trzeba rozwiązać pomocniczy układ równań: 
 

Ly=Pb, 

a uzyskany y wstawić do równania 

Ux=y. 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

PRZYKŁAD - Rozwiązanie układu równań liniowych 
A=[10 2 1;0 40 2; 6 1 10]; 
b=[20 30 40]’; 

-

 

to musi być wektor kolumnowy 

 

%Sprawdzenie rzędu macierzy 

rank(A)=3 
rank([A b])=3 

%Rozwiązanie układu równań 

x=A\b 

x = 

1.5808 
0.6004 
2.9915 

Funkcja ta rozwiązują układ równań dwukrotnie szybciej niż polecenie: 

x=inv(A)*b 

 

inv(A) - wyznacza odwrotność macierzy A 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

ZASTOSOWANIE ROZKŁADU LU 
DO ODWRACANIA MACIERZY 
A=PLU   A

-1

=(PLU)

-1

=U

-1

L

-1

P

-1

=U

-1

L

-1

 
PRZYKŁAD 
A=[1 2 1;0 4 2; 6 1 1]; 
inv(A)= 
    1.0000   -0.5000   -0.0000 
    6.0000   -2.5000   -1.0000 
  -12.0000    5.5000    2.0000 
[L,U,P]=lu(A); 
A=P*L*U 
invA=inv(U)*inv(L)*P 
invA = 
    1.0000   -0.5000   -0.0000 
    6.0000   -2.5000   -1.0000 
  -12.0000    5.5000    2.0000 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

MATLAB & Scilab

 

ZASTOSOWANIE ROZKŁADU LU 
DO WYZNACZANIA WARTOŚCI WYZNACZNIKA MACIERZY 
 
PRZYKŁAD 
A=[1 2 1;0 4 2; 6 1 1]; 
d=det(A)=2 
 
[L,U,P]=lu(A) 
wyznacznik A jest równy iloczynowi elementów na głównej przekątnej U 

det(A)=det(U)= 

 

1

n

ii

i

u

=

a=diag(U)’ 
 
d=abs(a(1)*a(2)*a(3))=2 
 
 
 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

10 

background image

MATLAB & Scilab

 

ROZKŁAD CHOLEWSKY’EGO 
 
Obok znajdującego różnorodne zastosowanie rozkładu macierzy A=LU należy 
wspomnieć o rozkładzie Cholewsky’ego zwanym także rozkładem 
Banachiewicza. Polega on na zapisaniu macierzy A  (dodatnio określonej - 
zawsze odwracalna i jej odwrotność jest również dodatnio określonaw postaci 
iloczynu LL

T

, gdzie L jest macierzą trójkątna dolną. 

Rozkład ten realizuje funkcja: 

L=chol(A) 

PRZYKŁAD 
A=[10 2 1; 2 40 2; 6 1 10]; 
L=chol(A) 
 
L = 
    3.1623    0.6325    0.3162 
             0    6.2929    0.2860 
             0             0    3.1334 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

11 

background image

MATLAB & Scilab

 

 
ROZKŁAD NA MACIERZ ORTOGONALNĄ I TÓJKĄTNĄ GÓRNĄ 
 
Funkcja: 

[Q,R]=qr(A) 

A=QR 
gdzie: 
R - macierz trójkątna górna 
Q - macierz ortogonalna 
 
PRZYKŁAD 
A=[10 2 1; 2 40 2; 6 1 10]; 
[Q,R]=qr(A) 
Q = 
   -0.8452    0.1427   -0.5151 
   -0.1690   -0.9856    0.0043 
   -0.5071    0.0907    0.8571 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

12 

background image

MATLAB & Scilab

 

R = 
  -11.8322    -8.9586   -6.2541 
              0  -39.0480   -0.9212 
              0              0    8.0646 
 
Sprawdzenie ortogonalności macierzy Q 
 

Q

T

Q=QQ

T

=I 

 

Q’*Q= 
    1.0000    0.0000    0.0000 
    0.0000    1.0000    0.0000 
    0.0000    0.0000    1.0000 
Q*Q’= 
    1.0000    0.0000    0.0000 
    0.0000    1.0000    0.0000 
    0.0000    0.0000    1.0000 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

13 

background image

MATLAB & Scilab - 

DODATEK

 

WEKTORY I WARTOŚCI WŁASNE MACIERZY 
 
Jeżeli pewien wektor x zostanie poddany przekształceniu opisanemu przez 
macierz w wyniku którego otrzymamy wektor yy=A*x. Jeżeli y=lambda*x 
to wtedy x jest wektorem własnym macierzy A, zaś  lambda jest jej wartością 
własną. Spełnione są związki: 
 

det(A-lambda*I)=0 

oraz 

(A-lambda*I)*X=0 

 
Numeryczne wyznaczanie wartości własnych i wektorów własnych nie jest 
zadaniem  łatwym, a dostępne algorytmy ich uzyskania nie są niezawodne 
zwłaszcza gdy mamy do czynienia z wielokrotnymi wartościami własnymi 
macierzy lub wartościami w postaci liczb zespolonych. 
 
 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

14 

background image

MATLAB - 

DODATEK

 

 
 
PRZYKŁAD
 
Obliczyć wartości własne macierzy 
A=[3 -1 1; -1 5 -1; 1 -1 3]; 

 

Współczynniki wielomianu charakterystycznego według malejących potęg 
p=poly(A) 
p = 
1.0000  -11.0000   36.0000  -36.0000 
Wartości własne macierzy 
lambda=roots(p)’ 
lambda = 
6.0000    3.0000    2.0000 
 
 
 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

15 

background image

MATLAB & Scilab - 

DODATEK

 

 
Wartości własne można też wyznaczyć przy pomocy polecenia 
 

lambda=eig(A) 

lub 

[X,lambda]=eig(A) 

gdzie: 
X - macierz wektorów własnych macierzy 
lambda - 
macierz diagonalana, której główna przekątna zawiera wartości własne 
 
PRZYKŁAD
 
A=[3 -1 1; -1 5 -1; 1 -1 3]; 
lambda=eig(A)’ 

% MATLAB 

lambda=spec(A)’; 

// Scilab 

 

lambda=mtlb_eig(A)’ 

// Scilab 

funkcja niedostęna już od wersji 5.1.2

 

lambda= 
2.0000    3.0000    6.0000 

 

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

16 

background image

Metody Numeryczne 

 4. Układy równań liniowych 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

17 

MATLAB - 

DODATEK

 

[X,lambda]=eig(A) 
 
X = 
    0.7071   -0.5774    0.4082 
   -0.0000   -0.5774   -0.8165 
   -0.7071   -0.5774    0.4082    

 

 

 

 

wektory własne macierzy A

 

lambda = 
    2.0000         0         0 
         0    3.0000         0 
         0         0    6.0000  

wartości własne macierzy A na głównej przekątnej 

 
 
 
 
 
 
 

 

background image

 

METODY NUMERYCZNE 

 

Temat 5 

RÓWNANIA RÓŻNICZKOWE 

 

 
Ćwiczenia opracowano na podstawie: 

1)  Zalwewski A., Cegieła R. „MATLAB-obliczenia numeryczne i ich 

zastosowanie” Nakom, Poznań 1997. 

2)  Jankowski R., Lubowiecka I., Witkowski W. „Podstawy programowania w 

języku Matlab” Politechnika Gdańska, WIL, 2002 

3)  http://www.mathworks.com/access/helpdesk.html 
4)  Brozi A. Scilab w przykładach” Nakom, Poznań 2007 

 

 
 

 

background image

Scilab

 

Rozwiązywanie równań różniczkowych zwyczajnych 

Równania pierwszego rzędu 
Weźmy pod uwagę równanie opisujące ruch jednostajny z prędkością v wzdłuż 
osi x 

dx

v

dt

=

Rozwiązaniem jest funkcja: 

0

( )

x t

v t x

= ⋅ +

FUNKCJA BIBLIOTECZNA 

wynik=ode(x0,t0,t,f) 

gdzie: 
x0 i t0 - warunki początkowe 
t - chwila czasu odpowiadająca rozwiązaniu 
f - funkcja określająca równanie różniczkowe 
Procedura służy do rozwiązywania równań różniczkowych postaci: 

( , )

dx

f t x

dt

=

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab

 

PRZYKŁAD 
 

//Program-różniczka funkcji jednej zmiennej 
//czyszcenie pamięci 

clc() 
clear 
xdel(winsid()) 

// 

v=1; 

//prędkość początkowa 

x0=0; 

//położenie początkowe

 

t0=0; 

//czas początkowy

 

tk= 10; 

//czas końcowy 

t=linspace(t0,tk); 

//deklaracja wektora czasu 

//deklaracja funkcji f(t,x)  

function wynik=f(t,x) 
   wynik=v; 
endfunction 

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab

 

// 

x=ode(x0,t0,t,f); 

// 
//wykres położenia od czasu 

xinit() 

//inicjalizacja draiverów grafiki

 

plot2d(t,x) 
xtitle('Zaleznosc polozenia od czasu w ruchu jednostajnym','t (czas)','x 
(polozenie)') 

//deklaracja zmiennych graficzych 

w=gca(); 
w.title.font_size=4; 
w.x_label.font_size=4; 
w.y_label.font_size=4; 
w.children.children.thickness=3; 
w.children.children.foreground=2;  
 
 

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab

 

Wynik pracy programu 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

Scilab

 

Równania drugiego rzędu 
W przypadku równań rzędu drugiego (lub n) zawsze możemy zapisać to 
równanie jako układ dwóch równań (lub n równań), jak poniżej: 

2

2

d x

a

dt

=

 

 

dx

v

dt

dv

a

dt

=

⎪⎪

=

⎪⎩

 

 

x

v

d

v

a

dt

⎡ ⎤ ⎡ ⎤

=

⎢ ⎥ ⎢ ⎥

⎣ ⎦ ⎣ ⎦

Układ dwóch równań drugiego rzędu 
Rzut ukośny „szkolny” (zmiana dwóch współrzędnych x i y) 
 

2

2

d r

a

dt

=

 

 

dr

v

dt

dv

a

dt

=

⎪⎪

=

⎪⎩

 

 

x

y

x

x

y

y

x

v

y

v

d

v

a

dt

v

a

⎡ ⎤ ⎡ ⎤

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

=

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

⎣ ⎦ ⎣ ⎦

 

background image

Scilab

 

PRZYKŁAD 

 
//Program - rzut szkolny 
// 

clc() 
clear 
xdel(winsid()) 

// 

deg=%pi/180; 

//jeden stopień w rad

 

g=[0;-9.81]; 

//[m/s^2] wektor przyśpieszenia ziemskiego 

r0=[0;0]; 

//[m] wektor początkowy punktu rzutu

 

kat=[40:5:50]*deg; 

//trzy wyjściowe kąty wyrzutu

 

wv0=50; //[m/s] 

//prędkość początkowa rzutu 

v0=[cos(kat);sin(kat)]*wv0; 
t0=0;tk=10;t=linspace(t0,tk); 

 
 

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab

 

//deklaracja funkcji 

function wynik=szkolny(t,Y) 
  wynik=zeros(4,1); 

//incjacja wektora w postaci zer

 

  wynik(1:2)=Y(3:4); 

//piersze dwa elementy to 3 i 4 element argumentu Y 

(prędkości) 

  wynik(3:4)=g; 

//przyśpieszenia (3 równa zero, 4 równa -9.81)

 

endfunction 
x40=ode([r0;v0(:,1)],t0,t,szkolny); 
x45=ode([r0;v0(:,2)],t0,t,szkolny); 
x50=ode([r0;v0(:,3)],t0,t,szkolny); 
xinit() 
plot2d([x40(1,:)',x45(1,:)',x50(1,:)'],[x40(2,:)',x45(2,:)',x50(2,:)']) 
xtitle('Trajektoria rzutu szkolnego','x (polozenie) [m]','y (wysokosc) [m]') 
legend('40','45','50') 
w=gca();w.title.font_size=4;w.font_size=3; 
w.x_label.font_size=4; 
w.y_label.font_size=4; 

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

background image

Scilab

 

Metody Numeryczne 

 5. Równania różniczkowe 

• Opracował: Marcin Kujawa • Politechnika Gdańska • WILiŚ • Katedra Mechaniki Budowli i Mostów 

 

Wynik pracy programu 
 

background image

ZADANIE DOMOWE nr 1 

Dokonaj aproksymacji mając dany zbiór danych empirycznych przedstawionych w 
tabeli x i y 

x y 
2.0 12.0 
1.2 8.0 
14.8 76.4 
8.3 17.0 
8.4 21.3 
3.0 10.0 
4.8 12.5 
15.6 97.3 
11.5 25.0 
14.2 38.6 
14.0 47.3 
16.1 88.0 

MATLAB 
1) Aproksymuj dane prostą oraz wielomianami stopnia 2, 3 i 4 wykorzystując dostępne funkcje języka 
Matlab. 
2) Określ błędy korelacji i determinacji każdej z wykonanych aproksymacji. 
3) Napisz własną procedurę aproksymacyjną stopnia 2. 

Scilab 

1) Napisz własną procedurę aproksymacyjną stopnia 2. 
2) Określ błędy korelacji i determinacji każdej z wykonanych aproksymacji. 

background image

ZADANIE DOMOWE nr 2 

 
Napisz program poszukujący rozwiązania równania 

( )

x

f x

=

, implementując metodę 

iteracyjną Newtona-Raphsona postaci: 
 

1

( )

( )

i

i

i

i

f x

x

x

f x

+

= −

 
z dokładnością do błędu szacowanego zgodnie z formułą: 
 

1

1

100%

i

i

i

x

x

x

ε

+

+

=

background image

ZADANIE DOMOWE nr 3 

Napisz program poszukujący rozwiązania układu równań algebraicznych postaci 

{ } { }

[ ]

A X

C

=

, implementując metodę iteracyjną Gaussa-Sidela. 

 
Należy wyznaczyć niewiadome z równań: 

1

12

2

13

3

1

1

11

2

21

1

23

3

2

2

22

1

1

2

2

,

1

1

...

,

...

,

...

.

n

n

n

n

n

n

n

n n

n

n

nn

C

a X

a X

a X

X

a

C

a X

a X

a X

X

a

C

a X

a X

a

X

X

a

− −

=

− −

=

− −

=

#

 

Następnie przyjmując 

yznaczyć 

2

3

0

n

X

X

X

=

=

=

 w

1

1

zem z pozostałymi podstawić do równania 

nr 2 wyznaczając 

2

 i tak dla wszystkich niewiadomych do 

n

. W podobny sposób wykonujemy kolejne 

iteracje do momentu gdy błąd jest mniejszy od przyjętego zgodnie z formułą: 

11

C

X

a

=

 i ra

X

X

1

100%

j

j

i

i

j

i

x

x

x

ε

background image

ZADANIE DOMOWE nr 4 

 
Napisz program poszukujący rozwiązania równania 

( )

x

f x

=

, implementując metodę 

iteracyjną bisekcji (połowienia). 

Pierwiastek poszukujemy w miejscu gdzie funkcja zmienia znak, a więc: 

 

( ) ( ) 0

l

u

f x f x

<

 dla sąsiednich punktów  

l

x

   i   

u

x

  

Najpierw poszukujemy metodą przyrostową zmiany znaku funkcji i następnie metodą bisekcji poszukujemy 
pierwiastka funkcji. 

1. Dla 

odciętych 

l

x

  i  

u

x

, dla których 

( ) ( ) 0

l

u

f x f x

<

 obliczamy 

2

l

u

r

x

x

x

+

=

 

2. Następnie sprawdzamy: 

a) Jeżeli 

( ) ( ) 0

l

r

f x f x

<

 pierwiastek znajduje się w lewej części przedziału więc podstawiamy 

u

r

x

x

=

b) Jeżeli 

 pierwiastek znajduje się w prawej części przedziału więc podstawiamy 

( ) ( ) 0

l

r

f x f x

>

l

u

x

x

=

c) Jeżeli 

( ) ( ) 0

l

r

f x f x

=

     

r

x

 jest poszukiwanym pierwiastkiem. 

Działania powtarzamy, aż wyznaczymy pierwiastek z żądaną dokładnością. 


Document Outline