Metody numeryczne w nauce i technice.

MATLAB

Laboratorium

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Plan zajęć:

Wprowadzenie – informacje wstępne:

1. Wstęp do programowania w środowisku Matlab.
2. Macierze, wektory, układy równań liniowych.
3. Interpolacja i aproksymacja.
4. Przybliżone rozwiązywanie równań nieliniowych.
5. Całkowanie i różniczkowanie numeryczne.
6. Rozwiązywanie równań różniczkowych zwyczajnych.
7. Kolokwium zaliczeniowe.

1.0 Wprowadzenie do metod numerycznych.

Metody numeryczne

są działem matematyki, zajmującej się opracowywaniem metod

przybliżonego rozwiązywania zagadnień matematycznych, których rozwiązanie sposobami
analitycznymi (ścisłymi) jest trudne, albo wręcz niemożliwe.

1.1 Oszacowanie dokładności obliczeń i analiza błędów.

Błędy opisu matematycznego analizowanego zjawiska

Błędy początkowe

(błędy danych wejściowych)

Błędy obcięcia

(skończona liczba działań)

Błędy zaokrągleń

(niedokładna reprezentacja liczb rzeczywistych)

Zadania źle uwarunkowane

– niewielka względna zmiana danych powoduje duże względne

zmiany jego rozwiązania.

Algorytm jest

numerycznie poprawny

, jeżeli dla dokładnych danych wejściowych daje

rozwiązania będące rozwiązaniami zadania dla niewiele zaburzonych tych danych.

(nie dotyczy deterministycznych układów chaotycznych☺)

1.2 Podstawowe komendy Matlaba

%

- komentarz, objaśnienia

clc

- czyści okno poleceń Command Window

clf

- czyści okno graficzne Figure

clear x

- usuwa zmienną x z pamięci

clear all

- usuwa wszystkie zmienne z pamięci

who,whos

- wyświetla zdefiniowane zmienne

lenght(x)

- zwraca długość wektora

size(M)

- zwraca rozmiar macierzy, ilość wierszy i kolumn

t=t1:Dt:t2

- definiuje wektor o poczatku t1, krok Dt, końcu t2.

M’

- transpozycja macierzy M

A*B

- iloczyn macierzy A i B

sum(x)

- suma elementów wektora x

prod(x)

- iloczyn elementów wektora x

det(A)

- wyznacznik macierzy kwadratowej A

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

2

Laboratorium z metod numerycznych – Informatyka Stosowana I .

lu

- rozkład macierzy na macierze trójkątne

[L,U]-lu(A)

- L*U=A

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

- L*U=PA

-

L

– macierz trójkątna dolna

-

U

– macierz trójkątna górna

-

P

– macierz permutacji

eye(n)

- macierz jednostkowa o wymiarze (n x n)

inv(A)

- macierz odwrotna

eig(A)

-

wektor

wartości własnych macierzy kwadratowej A

[V,D]=eig(A)

-

D

- diagonalna macierz wartości własnych A

-

V

– macierz, której kolumny odpowiadają wektorom własnym

spełniającym A*V=V*D

ones(m,n)

- macierz jedynek (m x n)

zeros(m,n)

- macierz zer (m x n)

2.0 Macierze, wektory, układy równań liniowych.

M-plik funkcyjny –

nazwa funkcji.m

Kod:

fun1.m

----------------------------

function y=fun1(x)

% opis fukcji

y=4*sin(3*x).*exp(-x/2);

============================

p2.m

----------------------------

clear all

t=1:0.2:10;
y=2*sin(2*t-1);

b=fun1(t);
plot(t,b,’xr’,t,y), grid on
pause(3)

grid off

Wynik przedstawia się następująco:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

3

Laboratorium z metod numerycznych – Informatyka Stosowana I .

2.1 Macierze i wektory

M1=[1

2

3

W1=[1

2

3]

4 5 6]

M1=

W1=

1 2 3

1 2 3

4 5 6

W2=[1

2

3]’

M2=[1 2 3;4 5 6;7 8 9]

W2=

1

M2=

2

1 2 3

3

4 5 6

7 8 9

W3=1:2:8

W3=

1 3 5 7

2.2 Rozwiązywanie układów równań liniowych

Układ równań liczbowych

a

11

*x

1

+a

12

*x

2

+…+a

1k

*x

k

+…+a

1n

*x

n

=b

1

a

21

*x

1

+a

22

*x

2

+…+a

2k

*x

k

+…+a

2n

*x

n

=b

2

a

k1

*x

1

+a

k2

*x

2

+…+a

kk

*x

k

+…+a

kx

+x

n

=b

k

a

n1

*x

1

+a

n2

*x

2

+…+a

nk

*x

k

+…+a

xn

*x

n

=b

n

Przykład liczbowy

( )

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

=

=

≠

=

−

=

−

+

+

−

=

−

+

−

+

−

=

+

+

−

)

(

)

(

*

0

)

det(

*

15

4

2

4

11

4

2

9

7

2

21

2

6

6

7

3

3

2

1

1

3

x

2

x

1

x

x

x

x

z

y

x

x

B

A

x

B

x

A

A

x

z

y

x

x

z

y

x

y

x

z

y

x

Zapis w Matlabie

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

4

Laboratorium z metod numerycznych – Informatyka Stosowana I .

2.3

Rozkład LU

Kod:

% Równanie

A * x = B

(1)

% Rozkład LU

P * A = L * U

% Równanie (1) obustronnie mnożymy przez P

P * A *x = P * B

%

Stosujemy rozkład LU

L * U *x = P * B

% Wprowadzamy nową zmienną

U * x = y

% Otrzymujemy równanie o macierzy trójkątnej

L * y = P * B

Przykład m2.m

Kod:

A=[1 2 3;2 3 4;3 4 6];
B=[7;8;9];

[L,U,P]=lu(A);
n=length(B);

PB=P *B;

% dla macierzy L

y(1)=PB(1)/L(1,1);

for k=2:n
y(k)=(PB(k)-L(k,1:k-1)*y(1:k-1))’)/L(k,k);

end

% dla macierzy U

x(n)=y(n)/U(n,n);
for k=1:n-1

x(n-k)=(y(n-k)-U(n-k,n-k+1:n)*x(n-k+1:n))’)/U(n-k,n-k);
end
x’

2.4 Rozwiązania w Matlabie

Macierz kwadratowa → A * x = B

det(A) ≠ 0

Kod:

x=A

-1 *

B

x=inv(A) * B

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

5

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Inne przypadki: układy nadokreślone i niedookreślone

x=A\B

%Matrix left division \ of Galois arrays

x=pinv(A)B

%Moore-Penrose pseudoodwrotna macierz

Przykład 1

Rozwiąż równanie

Lx –B

T

=2* Ex

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

−

−

=

6

5

0

9

2

1

4

7

5

A

[

]

30

20

10

−

=

B

L – jest macierzą trójkątna dolna po permutacji macierzy A
E – jest macierzą jednostkową o wymiarach macierzy A

Kod:

A=[5 -7 4; =3 2 9; 0 5 6];

B=[-10 20 30];
[L,U,P]=lu(A);

E=eye(size(A));
X=(L-2*E)\B’

Otrzymujemy wynik:

x=

10.0000

-20.0000
-27.2000

2.5

Zadania

1.

4x+3y+6=14+3x-2y-8;
6y+15+3z=-4x+3;
2x-25+2y=3x+3y-3x;

Oblicz y oraz s = x

2

+y

2

+z

2

2.

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

6

Laboratorium z metod numerycznych – Informatyka Stosowana I .

7

⎤

⎢

⎢

⎢

⎣

⎡

−

−

=

9

2

3

4

7

5

A

⎥⎦

6

5

9

⎥

⎥

[

]

30

20

10

−

=

B

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

=

3

2

1

x

x

x

x

U – jest macierzą trójkątna górna macierzy A
E – jest macierzą jednostkową o wymiarach macierzy A
D – jest wyznacznikiem macierzy A

Rozwiąż równanie

Ux-5B

T

=d* Ex

Rozwiązania

1.

Kod:

A=[1 3 8;4 6 0;5 -1 -3];

B=[8;-15;25];
x=A\B;

y=x(2)
s=x’*x

Wynik:

y= -6.7000
s= 103.9400

2.

Kod:

A=[5 -7 4;-3 2 9;0 5 6];
B=[-10 20 30];
[L,U]=lu(A);

d=det(A);
E=eye(size(A));

x=(U-d*E)\(5*B’)

Wynik:

x=
-0.1397

0.2740

0.4109

2.6 Układ nadmiarowy – ilośc równań większa od ilości niewiadomych.

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Dla pomiarów [t

1

,y

1

]

t

1

0.1 0.3 0.5 0.9 1.2 1.6 2.2 2.8

y

1

2.5 2.3 2.7 2.9 2.5 2.0 1.4 1.2

Wyznaczyc współczynniki C

1

równania y=f(t)

2

-3t

3

-t

2

1

e

*

C

e

*

C

C

y(t)

+

+

=

Kod:

t=[0.1 0.3 0.5 0.9 1.2 1.6 2.2 2.8]’
y=[2.5 2.3 2.7 2.9 2.5 2.0 1.4 1.2]’
plot(t,y.’o’)

% y=C1+C2*exp(-t)+C3*exp(-3*t.^2)

A=[ones(size(t)) exp(-t) exp(-3*t.^2)]

C=A\y
pause

%dalej

–

“enter”

w oknie polecen matlaba

clf
tt=[0:0.01:3]’;

AA=[ones(size(tt)) exp(-tt) exp(-3*tt.^2)];
yy=AA*C;
plot(tt,yy,t,y,’o’),grid on, hold on

2.7 Układ niedookreślony – ilość równań mniejsza niż ilość niewiadomych

% R*x=B

m5.m

% R=[R1 R2],

R=[6

8

7

3

5

% x=[x1’ x2’]’

3 5 4 1 3

% R1*x1+R*x2=B

2 5 3 7 2]

% x1*x1=B-R2*x2

[n,m]=size(R)

R1=R(:,1:n)

R2=R(:,n+:m)

B=[2 1 4]’

x2=[1 1]’

x1=R1\)B-R2*x2)
x=[x1’

x2’]’

blad=B-R*x

3.0 Interpolacja i aproksymacja.

Interpolacja

Zadanie znalezienie funkcji F(x) przechodzącej przez zadane punkty (x,y) nazywamy

I N T E R P O L A C J Ą

(po łacinie inter – między, polus – węzeł)

Punkty (x,y) nazywamy węzłami interpolacji.

Istnieje funkcja:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

8

Laboratorium z metod numerycznych – Informatyka Stosowana I .

]

[

),

(

0 n

x

x

x

x

f

y

∈

=

Należy wyznaczyć funkcję F(x), taką aby:

y

x

F

=

)

(

3.1 Metoda Funkcji bazowych

F(x) jest liniową kombinacją n=1 funkcji bazowych

∑

=

=

n

i

i

i

n

x

a

x

F

x

x

x

0

1

0

)

(

)

(

),

(

)

(

),

(

ϕ

ϕ

ϕ

ϕ

K

Wprowadzając macierz bazową

φ

[

]

)

(

)

(

)

(

)

(

1

0

x

x

x

x

n

ϕ

ϕ

ϕ

φ

K

=

i macierz współczynników A

[

]

T

n

a

a

a

A

K

1

0

=

A

x

x

F

∗

=

)

(

)

(

φ

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

∗

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

n

n

n

n

n

n

n

n

n

y

y

y

a

a

a

x

x

x

x

x

x

x

x

x

x

F

x

F

x

F

M

M

K

L

L

K

K

M

1

0

1

0

1

0

1

1

1

1

0

0

1

1

0

0

1

0

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

ϕ

ϕ

ϕ

ϕ

ϕ

ϕ

ϕ

ϕ

ϕ

Y

x

x

F

Y

A

∗

∗

=

=

∗

−1

)

(

)

(

φ

φ

φ

Jeżeli funkcje bazowe są wielomianami to

interpolacja wielomianowa

np. interpolacja

wielomianami

Langrange`a

,

jednomianami

n

n

x

x

=

)

(

ϕ

3.2 Wielomiany

2

7

3

6

4

)

(

)

(

2

4

5

0

+

+

−

+

=

∈

=

∑

−

−

−

x

x

x

x

x

p

N

i

x

p

x

p

n

i

i

n

i

n

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

9

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Funkcje MATLABA

p= [4 6 0 -3 7 2]

% zapis wielomianu o współczynnikach p

y=polyval(p,x)

% oblicza wartość wielomianu dla x

p=polyfit(x,y,n)

% znajduje współczynniki wielomianu n-tego stopnia
aproksymującego dane (x

i

y

i

)

Kod:

xi=[ 2 4 6 8 9];

yi=[-5 0 2 -1 7];
plot(xi,yi,’o’),grid on,pause(2)

n=length(xi);
p=polyfit(xi,yi,n-1);
x=1:0.01:9.2;

y=polyval(p,x);
plot(xi,yi,’o’,x,y),grid on, pause(2)

pp=polyfit(xi,yi,8);
yy=polyval(pp,x);

plot(x1,yi,’o’,x,y,x,yy,’r’),grid on
axis([1 10 -10 12])

%dopasowanie osi wykresu

Interpolacyjne funkcje MATLABA

Matlab posiada rozbudowane funkcje interpolacyjne dedykowane do odpowiednich danych:

• Funkcje dla danych jednowymiarowych

interp1

• Funkcje dla danych dwuwymiarowych

inrerp2

• Funkcje dla danych wielowymiarowych

interp3, intern

Kod:

y=interp1(xi,yi,x,

’metoda’

)

%wykonuje interpolację w punktach określonych wektorami x dla zadanych węzłów (x

i

y

i

).

Elementy wektora

xi

muszą tworzyć ciąg rosnący.

metoda –

wybór metody interpolacji

linear

–

interpolacja funkcją łamaną

spline

–

interpolacja funkcjami sklejanymi

cubic

–

interpolacja wielomianami 3-go stopnia

nearest

– interpolacja wartościa najbliższych punktów

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

10

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Kod:

xi=[ 2 4 6 8 9];

yi=[-5 0 2 -1 6];
plot(xi,yi,’o’),grid on, pause(2)

x=1:0.01:9.1;
y1=interp1(xi,yi,x,’linear’);
plot(x,y1,xi,yi,’o’),grid on, pause(2)

y2=interp1(xi,yi,’o’,’spline’);
plot(x,y2,xi,yi,’o’),grid on, pause(2)

y3=interp1(xi,yi,x,’cubic’);
plot(x,y3,xi,yi,’o’),grid on, pause(2)

y4=interp1(xi,yi,x,’nearest’);
plot(x,y4,xi,yi,’o’),grid on, pause(2)
plot(x,y1,x,y2,x,y3,x,y4,xi,yi,’o’)

Wynik:

3.3

Aproksymacja

Aproksymacja

Zadaniem znalezienie funkcji F(x) z żądaniem, aby jak najbardziej zbliżała się do funkcji f(x) w
pewnym określonym przedziale (a,b) nazywamy

A P R O K S Y M A C J Ą.

Kryterium Aproksymacj metodą najmniejszych kwadratów jest minimalizacja całki:

∫

−

=

b

a

dx

x

F

x

f

Q

2

))

(

)

(

(

Dla przypadku funkcji dyskretnej, minimalizacja sumy:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

11

Laboratorium z metod numerycznych – Informatyka Stosowana I .

[

]

∑

=

−

=

a

i

i

i

x

F

y

Q

1

2

)

(

Wyznaczanie współczynników prostej regresji

(

)

(

)

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

=

⎟

⎠

⎞

⎜

⎝

⎛ −

⎟

⎠

⎞

⎜

⎝

⎛ −

∗

⎟

⎠

⎞

⎜

⎝

⎛ −

=

∗

−

∗

−

∗

=

=

∗

−

∗

+

⎟

⎠

⎞

⎜

⎝

⎛

−

=

∗

−

∗

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

∗

−

+

∗

−

=

−

=

=

⎟

⎠

⎞

⎜

⎝

⎛

−

∗

+

=

−

+

∗

=

∂

∂

=

∗

−

+

∗

=

∂

∂

−

+

∗

=

−

=

∆

=

−

=

∆

+

∗

=

=

n

i

i

n

i

i

i

n

i

n

i

i

i

n

i

i

i

n

i

i

n

i

i

i

n

i

i

n

i

n

i

i

i

n

i

i

n

i

i

n

i

i

i

n

i

i

n

i

i

n

i

n

i

n

i

i

i

i

i

i

n

i

i

i

n

i

n

i

i

i

i

n

i

i

i

i

i

x

x

y

y

x

x

a

x

x

x

x

y

y

x

a

y

x

x

y

x

x

x

a

y

x

x

x

a

y

x

a

x

a

y

b

x

n

a

y

n

b

y

b

n

x

a

y

b

x

a

b

b

a

Q

x

y

b

x

a

a

b

a

Q

y

b

x

a

y

y

y

b

a

Q

y

y

y

b

x

a

y

n

i

y

x

1

2

__

1

__

__

1

1

__

2

1

__

1

1

1

__

1

1

__

2

1

1

1

__

__

2

__

__

1

1

1

1

1

1

1

2

2

1

2

0

0

1

1

0

2

2

)

,

(

0

2

)

,

(

)

(

)

(

)

,

(

..

1

)

(

∑

=1

Kod:

xi=0:0.5:10;

n=length(xi);

y1=2*xi+3;
y2=3*(

rand

(1,n)-rand(1,n));

yi=y1+y2;
plot(xi,yi,’*’),pause(2)

xs=

sum

(xi)/n;

% wartość średnia x

ys=sum(yi)/n;

% wartość średnia y

a=(xi*yi’-n*xs*ys)/(xi*xi’-n*xs^2);
b=ys-a*xs;
x=0:0.01:10;

y=a*x+b;
plot(xi,yi,’*’,x,y),grid on

parametry_prostej=[a b]’

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

12

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Wynik:

parametry_prostej =

1.9143
3.9445

3.4 Zadania

Zadanie 1

Dokonano następujących pomiarów:

x

-5

-3

0 2 5 6 8

y 10 40 80 90 130 145 160

1. Wyznacz prostą regresji.
2. Oblicz y(1) i y(7).

Kod:

x=[ -5 -3 0 2 5 6 8];

y=[10 40 80 90 130 145 160];

c=polyfit(x,y,1)
y1_7=polyval(c,[1,7])

Wynik:

y= 11.2889x + 70.0444
y(1)= 81.3333

y(7)=149.0667

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

13

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Zadanie 2.

Wyznacz współczynniki wielomianu 4 stopnia aporoksymującego punkty pomiarowe(x,y), a
następnie oblicz wartość y(5).

x -2 0 1

3 6 8

y 1 3 0

-1 4 7

Kod:

x=[-2 0 1 3 6 8];

y=[1 3 0 -1 4 7];

c=polyfit(x,y,4)

%wyznaczanie współczynników

y5=polyval(x,5)

%obliczanie y(5)

xx=-2.5:0.1:9;

yy=polyval(c,xx);
plot(x,y,’o’,xx,yy,5,y5,’*r’),grid on

Wynik:
c=

-0.0187 0.2.447 -0.5070 -1.4437 2.4623

y5= 1.4386

4.0 Przybliżone rozwiązywanie równań nieliniowych.

Twierdzenie Bolzano-Cauchy`ego

Jeżeli f(x) jest ciągła w przedziale [a,b] i f(a)•f(b)<0 to w przedziale (a,b) znajduje się co
najmniej jeden pierwiastek równania f(x)=0

Twierdzenie

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

14

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Jeżeli w przedziale [a,b] spełnione są założenia tw. Bolzano-Cauchy`ego i dodatkowo

sgn(f’(x))=const

dla

x

∈

[a,b],

to przedział ten jest przedziałem izolacji pierwiastka

równania

f(x)=0

4.1

Metoda bisekcji:

0

)

(

)

(

3

2

3

2

1

1

2

1

<

∗

>

−

=

∨

−

=

+

=

−

−

x

f

x

f

k

k

j

k

j

x

x

x

k

k

k

Błąd metody:

k

k

x

x

x

x

2

1

2

0

−

<

−

=

δ

4.2

Metoda cięciw

)

(

)

(

)

(

)

(

)

(

)

(

1

1

1

1

1

1

1

1

x

f

x

f

x

x

x

f

x

x

x

f

x

f

x

x

x

f

x

x

k

k

k

k

k

k

k

−

−

−

=

−

−

−

=

+

+

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

15

Laboratorium z metod numerycznych – Informatyka Stosowana I .

4.3

Metoda stycznych (Newtona)

b

x

x

f

x

f

a

x

x

f

x

f

k

k

k

k

=

→

<

∗

=

→

<

−

1

1

0

)

(

"

)

(

'

0

)

(

"

)

(

'

k

x

f

x

f

x

x

k

k

k

k

≥

−

=

+1

1

)

(

'

)

(

4.4

Metoda iteracji dla równania x=F(x)

Twierdzenie

Jeżeli w przedziale [a,b] izolacji pierwiastka równania x=F(x) pochodna funkcji F’(x) spelnia

warunek

1

)

(

'

<

≤ M

x

F

to proces iteracji jest zbieżny, przy czym x

0

∈

[a,b].

4.5 Wyprowadzanie metody

)

(

)

(

)

(

)

(

)

(

0

)

(

0

)

(

1

k

k

x

F

x

x

F

x

x

x

f

a

x

F

x

x

x

f

a

x

f

a

x

f

=

=

+

∗

=

=

+

∗

=

∗

=

+

Warunki zapewniające zbieżność

a

x

f

a

x

f

a

x

F

x

F

→

<

+

∗

+

∗

=

<

1

1

)

(

1

)

(

)

(

1

)

(

'

4.6 Metoda iteracji

x

x

x

y

−

+

+

−

=

3

4

2

Kod:
Plik funkcyjny f7.m: [f(x)]

function y=f7(x)

y=-x.^2+4*x+3-sqrt(abs(x));

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

16

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Plik fukncyjny f8.m: [F(x)]

function y=f8(x)

global a
Y=a*f7(x)+x;

Plik skryptowy:

clear all,clc,clf

global a
a=0.3;

fplot('[f7(x),x,f8(x)]',[0 5 -0.5 6])
grid on

hold on
x(1)=1;
b=1;

k=1;
while b>0.001

x(k+1)=f8(x(k));

b=abs(f7(x(k+1)));

k=k+1;

if k>10, break, end

end

n=length(x);
zo=x(n)

plot(zo,0,'ok',x,f8(x),'ok')

Wynik:

4.7

Polecenia MATLABA

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

17

Laboratorium z metod numerycznych – Informatyka Stosowana I .

xo = fzero(‘

plik_funkcyjny’

,

x

0

)

% oblicza miejsce zerowe funkcji podanej w m-

pliku_funkcjnym

%

x

0

punkt startowy

r=roots(p)

% oblicza pierwiastki wielomianu o współczynnikach zadanych wektorem p

fplot(‘

plik_funkcyjny’

,limits)

% wykreśla funkcję podana w m-

pliku_funkcyjnym

w przedziale określonym przez

limits

limits

jest wektorem określającym oś x

– [ x

min

x

max

],

lub osie x i y –

[x

min

x

max

y

min

y

max

].

4.8

Zadania

Zadanie 1

Wyznacz miejsca zerowe funkcji f(x)=0 w przedziale [0,9]

10

2

7

2

.

0

2

3

cos

2

)

(

x

x

e

e

x

x

f

∗

+

∗

⎟

⎠

⎞

⎜

⎝

⎛

−

∗

=

π

Kod:

Plik funkcyjny:

function t=fun2(x)

y=2*cos(pi*x/3-2).*exp(-x/2)-0.2/sqrt(7).*exp(-x/10);

Plik skryptowy:

fplot(‘fun2’,[-1 10 -2 2]), hold on

p1=fzero(‘fun2’,2);

p2=fzero(‘fun2’,6);
p3=fzero(‘fun2’,8);

p=[p1 p2 p3]’

plot(p,zeros(size(p)),’ok’), grid on

Wynik:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

18

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Zadanie 2

Znaleźć pierwiastki wielomianów:

a) p1(x)=2x

3

+6x

2

-20x-30

b) p2(x)=-1.5x

5

+2x

4

+20x

3

-100x

2

+80x+300

Kod:

p1=[2 6 -20 -30];

p2=[-1.5 2 20 -100 80 300];
r1=roots(p1)

r2=roots(p2)
x=-5:0.1:4;

y1=polyval(p1,x);
y2=polyval(p2,x);

subplot(2,1,1),plot(x,y1,r1,0,'o'), grid
subplot(2,1,2),plot(x,y2,r2,0,'o'), grid

Wynik:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

19

Laboratorium z metod numerycznych – Informatyka Stosowana I .

5.0 Całkowanie i różniczkowanie numeryczne.


Dy=diff(y) %zwraca wektor różnic sąsiednich elementów wektora y

1

1

]

[

],

[

1

1

2

1

3

2

1

=

=

−

=

∆

∆

∆

∆

=

=

+

−

n

k

dla

y

y

y

y

y

y

Dy

y

y

y

y

y

k

k

k

n

n

K

K

K

Kod

x=0:0.1:7;
y=sin(x);
n=length(x)

yprim=diff(y)./fidd(x);

%pochodna

pot(x,y,x(1:n-1),yprim), grid

Kod f1.m

fumction y=f1(x)

y=2*x+5*cos(3*x).*exp(-x/2);

Kod

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

20

Laboratorium z metod numerycznych – Informatyka Stosowana I .

x=0:0.1:7;
y=f1(x);
n=length(x)

yprim=diff(y)./diff(x);

%pochodna

plot(x,y,x(1:n-1),yprim),grid

legend(‘funkcja’,’pochodna’);

5.1 Polecenia Matlaba

Dla wielomianów

P=[1 4 7 0 3 1 6]

P1=polyder(P)

Dla funkcji wymiernej

L=[-2 0 4 0 20 0]
M=[ 1 0 4 …]

[L1,M1]=polyder(L,M)

5.2 Przykłady

Wykreślić funkcję i jej pochodną dla x z przedziały -4 do 4

8

2

1

)

(

2

5

3

+

+

−

=

x

x

x

x

W

Rozwiązanie

h=0.05;
x=-4:h:4;

L=[1 0 -1];M=[1 0 2 0 8];
[L1,M1]=polder(L,M)

y=polyvel(L,x)./polyval(M,x);
yprim=polyval(L1,x)./polyval(M1,x);

plot(x,y,x,yprim,’-‘),grid
title(‘Funkcja wymierna i jej pochodna’)
legend(‘funkcja’,’pochodna’)

5.3

Całkowanie numeryczne

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

21

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Wzory do numerycznego całkowania funkcji jednej zmiennej nazywane są

kwadraturami.

Funkcji wielu zmiennych –

kubaturami

Kwadratury interpolacyjne i aproksymacyjne:
Kwadratury Newtona-Cotesa
Kwadratury Gaussa
Kwadratury Lobatto

Metoda prostokątów

∑

∑

−

=

−

=

=

≈

1

1

1

1

)

(

n

k

n

k

k

k

y

h

x

f

h

F

Metoda trapezów

∑

∑

−

=

−

=

+

⎟

⎠

⎞

⎜

⎝

⎛

+

+

=

+

≈

1

1

1

2

1

1

2

2

2

)

(

)

(

n

k

n

k

n

k

k

k

y

y

y

h

h

x

f

x

f

F

Metoda parabol

)

)

(

)

(

4

)

(

(

3

2

/

1

1

2

2

1

2

∑

=

+

−

+

+

≈

n

k

k

k

k

x

f

x

f

x

f

h

F

Metoda Prostokątów

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

22

Laboratorium z metod numerycznych – Informatyka Stosowana I .

5.4

Polecenia Matlaba

∫

=

h

a

dx

x

f

F

)

(

F=quad(‘funk

% Simpson`s

cja’,a,b)

F=quad8(funkcja’a,b)

%Newton-Cotes(dawna)

felp quad

Quadv, quadl, dblquad, triplequad

5.5

Zadania

Oblicz całkę oznaczoną z funkcji f(x) w granicach a=1 b=10 metodą prostokątów,
trapezów I wbudowaną Matlaka

2

*

)

3

cos(

5

2

)

(

x

e

x

x

x

f

+

=

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

23

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Rozwiązanie

h=0.01;

x=1:h:10;
n=length(x)

fplot(‘f1’,[1 10]), grid
y=f1(x);
a=x(1);

b=x(n);
F(1)=h*sum(y(1:n-1));

%prostokątów

F(2)=h*(y(1)+2*sum(y(2:n-1))+y(n))/2;

%trapezów

F(3)=quad(‘f1’,a,b);

%Simpson`s

Wynik=F’

Oblicz całkę metodą prostokątów, trapezów i wbudowaną w Matlaba oraz wyznacz błąd
obliczeń numerycznych

∫

=

π

0

)

sin( dx

x

F

Rozwiązanie

h=0.1;
x=0:h:pi;

n=length(x)
fplot(‘sin’,[0 3*pi]), grid

y=sin(x);
a=x(1);

b=x(n);
format short
F1=h*sum(y:n-1));

%prostokątów

F2=h*(y(1)+2*sum(y(2:n-1))+y(n))/2;

%trapezów

F3=quad(‘sin’,a,b);

%Simpson`s

Wynik=[F1,F2,F3]’
Blad=2-Wynik

Wynik obliczeń

Wyniki=

∫

=

=

π

0

2

)

sin( dx

x

F

1,9954
1,9975
1,9991

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

24

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Błąd=
0,0046

metoda prostokątów

0,0025

metoda trapezów

0,0009

metoda Matlab

Oblicz pole figury leżącej w II-cwiartce układu współrzędnych ograniczonej wielomianem
W(x)=-x

4

+3x+4 oraz osiami OX i OY. Podaj pierwiastki wielomianu W(x), naszkicuj

wykres W(x) dla x z przedziału [-1.5 2]

Rozwiązanie

f2-inline(‘-x.^4+3*x+4’,’x’)

fplot(f2,[-1.5 2]), grid on
p2=[-1 0 0 3 4]

r2=roots(p2)
F2=quad(f2,r2(4),0)

z1=fzero(f2,-1);
z2=fzero(f2,2);
z=[z1;z2];

6.0 Równania różniczkowe zwyczajne (ordinary Differentia Equation ODE)

Równanie n-tego rzędu

(

)

)

1

(

)

2

(

)

1

(

)

(

,...

,

,

,

)

(

−

=

n

n

x

x

x

x

t

f

t

x

Równanie pierwszego rzędu

)

,

(

'

)

(

)

1

(

x

t

f

x

t

x

=

=

• Zagadnienia brzegowe (Boundary Valne Problem BVP)

• Równania różniczkowe cząstkowe (Partial Differentia Equation PDE)

6.1 Metody numeryczne

Numeryczne rozwiązywanie układów równań różniczkowych polega na poszukiwaniu
rozwiązania, startując z punktu, w którym jest ono znane (warunek początkowy). Jeśli znamy
rozwiązanie w punktach t

1

do t

k

to dla kolejnego t

k+1

wyznaczamy y(t

k+1

) stosując odpowiednią

aproksymację.

Prosta metoda Eulera

)

,

(

)

(

)

,

(

1

1

y

x

f

dx

x

dy

y

x

f

x

y

k

k

k

k

=

=

∆

∆

+

+

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

25

Laboratorium z metod numerycznych – Informatyka Stosowana I .

x

y

x

F

y

y

y

y

const

h

x

x

x

x

y

x

y

y

y

y

k

k

k

k

k

k

k

k

k

k

∆

+

=

=

=

=

∆

=

−

=

∆

=

−

=

∆

+

+

+

+

+

*

)

,

(

)

(

1

0

1

1

1

0

0

1

1

6.2 Rozwiązać numerycznie równanie w przedziale 0,3

b

y

a

dx

dy

+

−

=

*

Rozwiązanie analityczne

ax

e

a

b

y

a

b

x

y

−

−

+

=

*

)

(

)

(

0

Dane:

a=2, b=6, y(0)=y

0

=-1

Kod:

clear all, clf
f=inline(‘-2*y+6’,’x’,’y’);

x(1)=0;

% xo

y(1)=-1;

% yo

xk=3;
dx=0.01;

%krok

n=floor((xk-x(1))/dx);

%zaokraglenie

for k=1:n
y(k+1)=y(k)+f(x(k),y(k))*dx;

x(k+1)=x(k)+dx;
end

x1=min(x);x(2)=max(x);y1=min(y);y2=max(y);m=0.2;
yd=3+exp(-2*x)*(-3+y(1));

% rozw. dokładne

b=(y-yd)/(y2-y1)*100;

% blad w procentach

b1=min(b);b3=max(b);
subplot(2,1,1)

plot(x,y,x,yd),grid on, title(‘rozwiazania’)
axis([x1 x2 y1-m y2+m])

subplot(2,1,2)
plot(x,b),grid on,title(‘blad wzgledny w %’)

axis([x1 x2 b1-m b2+m])

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

26

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Wynik:

6.3

Funkcje ODE w Matlabie

[T,Y] = odexx(‘F’, przedzial, y0)

xx – liczba okreslajaca metode

powyższa instrukcja całkuje układ równań różniczkowych postacy y

n

=F(t,y)

w przedzale od

t_pocz

do

t-koniec

przedzial

= [

t_pocz t_koniec

]

z warunkami początkowymi podanymi w wektorze y0

‘F’ jest nazwą ODE pliku funkcyjnego.

Function

F(t,y) musi zawracać kolumnowy wektor.

[T Y] Każdy wiersz w macierzy rozwiązania Y odpowiada argumentowi t zwracanemu w
kolumnowym wektorze T.

Układy równań dobrze uwarunkowanych

ode45 - Runie-Kutta rzędu (4,5)

ode23 - Runie-Kutta rzędu (2,3)
ode113 – Adams-Bashforth-Moulton

Układy równań źle uwarunkowanych (sztywnych)

ode15s - Gear`s

ode23s - Rosenbrock
ode23t - trapez
ode23tb

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

27

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Przykład:
Rozwiązać numerycznie układ równań różniczkowych

0

)

0

(

5

.

0

1

)

0

(

2

2

2

1

1

1

2

=

−

−

=

=

=

y

y

y

dt

dy

y

y

dt

dy

w przedziale t [0,15], przy zadanych warunkach początkowych. Wykreślić przebiegi czasowe i
trajektorię. Podać wartość y

1

i y

2

dla t=1 i t=5

⎥

⎦

⎤

⎢

⎣

⎡

=

⎥

⎦

⎤

⎢

⎣

⎡

=

⎥

⎦

⎤

⎢

⎣

⎡

−

−

=

0

1

)

0

(

*

5

,

0

1

2

0

2

1

y

y

y

y

y

dt

dy

Pliki funkcyjny równanie różniczkowego row1.m

Kod:

function dy=row1(t,y)
dy=zeros(2,1);
dy(1)=2*y(2);

dy(2)=-y(1)-0.5*y(2);

Plik skryptowy rr1.m

Kod:

[t,y]=ode45(‘row1’,[0 14],[1;0]);
subplot(2,1,1)

plot(t,y(:,1),t,y(:,2)),grid on
title(‘przebiegi czasowe’)

subplot(2,1,2)
plot(y(:,1),y(:,2),grid on
title(‘trajektoria’)

y_p=interp1(t,y,[1,5],’linear’)

Wynik:

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

28

Laboratorium z metod numerycznych – Informatyka Stosowana I .

y_p =

0.2550 -0.1284

6.4 Równanie różniczkowe n-tego rzędu

0.2760 -0.5477

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

−

−

1

1

,

,

,

,

n

n

n

n

dt

y

d

dt

dy

y

t

f

dt

y

d

L

podstawiamy

1

1

2

1

−

−

=

=

=

n

n

n

dt

y

d

y

dt

dy

y

y

y

L

otrzymujemy uk

ład równań pierwszego rzędu

⎪

⎪

⎪

⎪

⎩

⎪⎪

⎪

⎪

⎨

⎧

=

2

1

y

dy

=

=

)

,

,

,

,

(

2

1

3

2

n

n

n

y

y

y

t

f

dt

y

d

y

dt

dy

dt

K

M

Przykład:

Rozwiąż równanie w przedziale 0 … 10

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

29

Laboratorium z metod numerycznych – Informatyka Stosowana I .

)

cos(

*

*

'

0

)

0

(

0

)

0

(

'

,

1

)

0

(

)

cos(

*

*

4

1

3

1

3

3

2

2

1

0

1

1

2

1

4

)

1

(

30

t

e

y

y

y

y

dt

d

y

y

dt

d

dt

y

y

y

y

y

y

y

y

y

t

e

y

y

n

n

+

+

−

=

=

=

=

=

=

=

−

=

=

−

−

−

6

y

y

+

y

y

d

=

.5

Równanie Van der Pola

al m

(1)^2)*y(2)];

0

'

*

)

1

(

2

=

+

−

+

y

y

y

y

n

µ

Kod:

function dy=row_vdp(t,y)

glob
dy=[y(2);-y(1)+m*(1-y

d:

r all

;

Ko

clea

global m
m=2;
[t,y]=ode45(‘row_vdp’,[0 20],[2:0]

plot(t,y(:,1)),grid on

e zadanie:

równanie różniczkowe

Kod:

function dx=de(t,x)

dx=[x(2);((2/t^2-1)*x(1))-x(2)/t];

Przykładow
Rozwiąż

Kod:

[t,x]=ode45(‘de’,[0.01 15],[1 0]);
plot(t,x(:,1));grid on
x_5=interp1(t,x(:,1),5,’linear’)

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

Laboratorium z metod numerycznych – Informatyka Stosowana I .

Przepisał: Blick © Materiały dostarczył: qWest © Pomagał: L3Vy ©

PWSZ Tarnów 2007

31

Wynik: