METODY NUMERYCZNE - LABORATORIUM

PODSTAWY JĘZYKA

MATLAB

(MATrix LABoratory)

i

Scilab

(Scientific laboratory)

Gdańsk 1

4.10. 2009

dr inż. Marcin Kujawa

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.

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

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/

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.

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.

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.

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

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

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

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

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

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

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

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

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() ?

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

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

2

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

3

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

4

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

5

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

6

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

7

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

8

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

9

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

MATLAB

Metody Numeryczne

•

1. Interpolacja i aproksymacja

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

11

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

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

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

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

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

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

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

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

2

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

3

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

4

Metody Numeryczne

•

2.Wielomiany, miejsca zerowe

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

5

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

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

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

2

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

3

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

4

Metody Numeryczne

•

3. Całkowanie numeryczne

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

5

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

Metody Numeryczne

•

3. Całkowanie numeryczne

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

6

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!

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

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

2

Metody Numeryczne

•

4. Układy równań liniowych

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

3

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 n

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.

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

4

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

5

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

6

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ć x 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

7

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

8

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

P

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

9

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

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ślona) w 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

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

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

MATLAB & Scilab -

DODATEK

WEKTORY I WARTOŚCI WŁASNE MACIERZY

Jeżeli pewien wektor x zostanie poddany przekształceniu opisanemu przez
macierz A w wyniku którego otrzymamy wektor y: y=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

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

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 A
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

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

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

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

2

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

3

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

4

Scilab

Wynik pracy programu

Metody Numeryczne

•

5. Równania różniczkowe

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

5

Metody Numeryczne

•

5. Równania różniczkowe

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

6

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

⎡ ⎤ ⎡ ⎤

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

=

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

⎢ ⎥ ⎢ ⎥

⎣ ⎦ ⎣ ⎦

.

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

7

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

8

Scilab

9

Metody Numeryczne

•

5. Równania różniczkowe

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

Wynik pracy programu

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.

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

ε

+

+

−

=

⋅

.

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

ε

−

−

⋅

≤

.

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ą.