Transcript MATwyk6
Slide 1
Wykład 6
https://play.google.com
Dr Aneta Polewko-Klim
Slide 2
Obliczenia numeryczne.
Wykonując obliczenia numeryczne zakładamy:
operujemy na liczbach, dopuszczając
możliwość występowania niedokładności w
ich reprezentacji,
podczas obliczeń mogą występować pewne
niedokładności.
Slide 3
Miejsca zerowe i minima funkcji.
Slide 4
Opcje
Wszystkie te funkcje można wywołać z opcjonalnym
parametrem opcje,
np.: x1 = fzero(‘cos(x) - 0.2’, 3, opcje)
Opcje określają parametry wywołania tych funkcji.
Ustawienia tych parametrów można zmienić poprzez
funkcję optimset.
opcje = optimset(parametr, wartość, ... )
Argumenty funkcji optimset to pary
parametr – wartość,
Parametr jest łańcuchem znakowym.
Slide 5
Funkcja optimset - lista parametrów i
ich możliwych wartości (kilka):
Slide 6
Wywołanie z nazwą funkcji
(wyświetla parametry dla konkretnej funkcji)
Slide 7
Przykłady parametrów funkcji
optimset:
Slide 8
Miejsca zerowe
Przykład:
>> fzero('cos(2*x)-.2',8)
ans =
8.7401
>> cos(2*ans)-.2
ans =
2.8866e-015
Slide 9
>> op=optimset('Diagnostics','on',...
'Display','Iter');
Looking for a zero in the interval [7.0949,
8.9051]
>> fzero('cos(2*x)-.2',8,op)
Func-count
x
f(x)
1
8
-1.15766
2
7.77373
-1.18715
3
8.22627 -0.935369
4
7.68
-1.14007
5
8.32
-0.7962
6
7.54745
-1.01789
7
8.45255 -0.565028
8
7.36 -0.750391
9
8.64
-0.19876
10
7.0949 -0.252615
11
8.9051
0.30677
Procedure
initial
search
search
search
search
search
search
search
search
search
search
12
7.91238
-1.19319
interpolation
13
8.70207 -0.0749546
interpolation
14
8.74193 0.00367162
interpolation
15
8.74007 2.49516e-005
interpolation
16
8.74006 -1.30267e-009
interpolation
17
8.74006 2.88658e-015
interpolation
18
8.74006 -4.05231e-015
interpolation
Zero found in the interval: [7.0949, 8.9051].
ans =
8.7401
Slide 10
Minimum – funkcja 1 zmiennej
Przykład
>>
funkcja=inline('sin(2*x)+cos(x)/2');
>>fplot(funkcja,[-2*pi, 2*pi])
>> x1=fminbnd(funkcja,-2,4)
x1 = 2.4375
>> funkcja(x1)
ans =
-1.3679
Slide 11
Minimum – funkcja 2 zmiennej
Do znalezienia minimum musimy
Podać postać tej funkcji lub
funkcję inline oraz punkt startowy.
Jak widać na rysunku funkcja
ta ma dwa minima. W zależności
od tego jaki punkt startowy
podamy, znajdziemy albo
pierwsze albo drugie.
Przykład
>> fun=inline('-exp(-(x(1)-2).^2-x(2).^2)-exp(-(x(1)+2).^2-x(2).^2)');
>> x0=fminsearch(fun,[2 2])
x0 = 2.0000 -0.0000
>> fx0=fun(x0)
fx0 = -1.0000
Slide 12
Pierwiastki wielomianów.
Ogólna postać wielomianu:
W(x) = a1xn + a2xn-1 + ... + anx + an+1
współczynniki wielomianu a1, a2, ... są uszeregowane według
malejących potęg zmiennej x.
Slide 13
Przykład:
W(x) = x2 + 3x – 4
>> roots([1 3 -4])
ans =
-4
1
>> poly(ans)
ans =
1 3 -4
W(x) = 3x5 – 2x4 + 4x3 – 6x + 1
>> roots([3 -2 4 0 -6 1])
ans =
0.2084 + 1.4463i
0.2084 - 1.4463i
-0.9198
1.0000
0.1697
Slide 14
Układy równań liniowych
Zapis macierzowy, każdy układ równań możemy przedstawić:
A*X = B
lub
X*A = B
Matlab dysponuje dwoma operatorami dzielenia, które odpowiednio
rozwiązują przypadki:
A*X = B mamy X = A\B – musi być równa liczba wierszy A i B
X*A = B mamy X = B/A – musi być równa liczba kolumn A i B
W technice częściej spotyka się pierwszy z tych przypadków.
Macierz A w ogólnym przypadku może mieć w wierszy i k kolumn.
Istnieją więc trzy możliwości:
w = k – macierz kwadratowa, jeśli jest nieosobliwa X można
wyznaczyć dokładnie;
w > k – układ jest nadokreślony;
w < k – układ jest nieoznaczony (rozwiązanie zawiera najwyżej w
niezerowych elementów;
Slide 15
Układ kwadratowy.
Mamy układ trzech równań
z trzema niewiadomymi.
Musimy zdefiniować
macierze A i B, a następnie
obliczyć wynik.
2 x1 3 x 2 4 x 3 5
x1 2 x 2 x 3 1
3x 2x x 5
1
2
3
Przykład
>> A=[2 -3 -4; -1 2 1 ; 3 -2 -1]
A=
2 -3 -4
-1 2 1
3 -2 -1
>> B=[-5 -1 5]'
B=
-5
-1
5
>> X=A\B
X=
2.0000
-1.0000
3.0000
Slide 16
Układ nieoznaczony cz.1
Liczba równań jest mniejsza
od liczby niewiadomych
6 x1 8 x 2 7 x 3 3 x 4 1
3 x1 5 x 2 4 x 3 x 4 2
Przykład
>> A=[6 8 7 3;3 5 4 1]
A=
6
8
7
3
5
4
>> B=[1 2]'
B=
1
2
>> X=A\B
X=
0
0.7143
0
-1.5714
3
1
Slide 17
Układ nieoznaczony cz.2
Ogólne rozwiązanie takiego
układu otrzymamy dodając
do otrzymanego rozwiązania
iloczyn macierzy Z będącej
jądrem macierzy A i
dowolnego wektora q.
Jądro: Z = null(A) A*Z = 0
Przykład cd.
>> Z=null(A)
Z=
-0.4253 -0.6519
-0.3981 0.3942
0.8126 -0.1607
0.0163 0.6276
>> q=[2 17]'
q= 2
17
>> X1=X+Z*q
X1 = -11.9324
6.6186
-1.1066
9.1306
Slide 18
Układ nadokreślony.
Gdy chcemy otrzymać n
współczynników funkcji
opisującej przebieg
jakiegoś zjawiska.
Dokonujemy wtedy np.
dużo więcej niż n
pomiarów i na ich
podstawie szukamy
wyniku.
x1 x 2 2 x 3 1
x1 2 x 2 x 3 2
3 x1 x 2 5 x 3 3
2 x1 2 x 2 3 x 3 4
Przykład
>> A=[1 -1 2; 1 -2 -1; 3 -1 5; -2 2 3]
A = 1 -1 2
1 -2 -1
3 -1 5
-2 2 3
>> B=[1 2 3 -4]'
B=1
2
3
-4
>> X=A\B
X = 1.4286
-0.1429
-0.2857
Slide 19
Interpolacja.
Zagadnienie dotyczy przechodzenia funkcji przez zadane punkty,
tzw. węzły interpolacji. Standardowe procedury Matlaba realizują
interpolację za pomocą wielomianów pierwszego i trzeciego
stopnia, oraz funkcjami sklejanymi trzeciego stopnia. Funkcje
interpolacyjne:
Slide 20
Metody interpolacji
‘nearest’— interpolacja najbliższy sąsiad
‘linear’ – interpolacja liniowa;
‘cubic’ – interpolacja wielomianem
trzeciego stopnia;
‘spline’ – interpolacja funkcjami sklejanymi
trzeciego stopnia;
Slide 21
Przykład:
x=0:.4:2;
y=exp(-.5*x).*cos(5*x);
xi=0:.05:2;
y1=interp1(x,y,xi,'nearest');
y2=interp1(x,y,xi,'linear');
y3=interp1(x,y,xi,'cubic');
y4=interp1(x,y,xi,'spline');
plot(x,y,'*',xi,exp(-.5*xi).*cos(5*xi),...
xi,y1,':',xi,y2,':',xi,y3,':',xi,y4,':')
legend('wezly','funkcja','nearest',…
'linear','cubic','spline')
Slide 22
Aproksymacja.
Zagadnienie aproksymacji polega na znalezieniu pewnej
funkcji, której wykres z pewnym przybliżeniem
przechodzi w pobliżu pewnych punktów.
Sama funkcja minimalizuje wartość pewnego kryterium
dopasowania, np. aproksymacja średniokwadratowa,
czyli minimalizujemy błąd średniokwadratowy.
Standardową metodą aproksymacji w Matlabie jest
aproksymacja wielomianem:
polyfit(x,y,n)
polecenie to znajduje wektor współczynników
wielomianu aproksymującego dane zawarte w wektorach
x i y. Wartość wielomianu możemy obliczyć funkcją
polyval.
Slide 23
Aproksymacja.
Przykład 1:
x=1:8;
y=[1.45 1.95 2.43 3.02 3.53 4 4.4 5.1];
a=polyfit(x,y,1)
xx=1:.01:8;
yy=polyval(a,xx);
plot(x,y,'o,xx,yy)
Otrzymaliśmy funkcję y=0.5121x+0.9304
Slide 24
Przykład 2:
x=0:.4:2*pi;
y1=sin(x);
blad=rand(1,length(x))/10;
y=y1+blad;
a1=polyfit(x,y,2)
a2=polyfit(x,y,3)
xx=0:.01:2*pi;
yy1=polyval(a1,xx);
yy2=polyval(a2,xx);
plot(x,y,'o',xx,sin(xx),':',xx,yy1,xx,yy2)
Dla wielomianu drugiego stopnia mamy:
y = – 0.0323x2 – 0.0982x + 0.7252
Dla wielomianu trzeciego stopnia mamy:
y = 0.0921x3 – 0.8611 x2 + 1.8274x – 0.0793
Slide 25
Pochodna funkcji
• Do znajdowania pochodnej możemy posłużyć się
funkcją:
diff
Zwraca ona wektor różnic sąsiednich elementów wektora
będącego jej argumentem.
• Czyli pochodna funkcji danej w postaci tablicy (np.
pomiary), obliczymy jako:
diff(y)./diff(x).
• Należy przy tym pamiętać, że liczba elementów
wektora zwracanego przez funkcję diff, jest o jeden
mniejsza, niż liczba elementów argumentu.
Slide 26
Pochodna funkcji
Przykład
x=72:.5:90;
y=[7237 7224 7203 7183 7164 7149 7129 7107 7060 7011 6912
6746 6517 ...
6225 5875 5449 4950 4376 3789 3202 2614 2082 1639 1290
1000 ...
758 581 470 394 328 284 241 191 146 103 52 0];
yprim=diff(y)./diff(x);
subplot(2,1,1),plot(x,y,'+'),grid on
ylabel('Dane')
subplot(2,1,2),plot(x, [yprim 0],'o'),grid on
ylabel('Pochodna')
Slide 27
Pochodna wielomianu
• Jeśli mamy dany wielomian, do obliczenia wartości
pochodnej możemy użyć funkcji:
polyder
• Zwraca ona współczynniki wielomianu pochodnego.
Wywołanie:
a_prim = polyder(a)
[L_prim, M_prim] = polyder(L,M)
Slide 28
Pochodna wielomianu
a=[-1 3 5 -2];
x=0:.05:3;
y=polyval(a,x);
a1=polyder(a);
y1=polyval(a1,x);
subplot(2,1,1);
plot(x,y)
grid on
ylabel('Dane')
subplot(2,1,2)
plot(x,y1)
grid on
ylabel('Pochodna')
Slide 29
Pochodna wielomianu - funkcja wymierna
Dodatkowo należy musimy zdefiniować licznik i mianownik:
Przykład:
L=[1 -3];
M=[1 5 -2 1];
x=-3:.05:3;
y=polyval(L,x)./polyval(M,x);
[L1,M1]=polyder(L,M);
y1=polyval(L1,x)./polyval(M1,x);
subplot(2,1,1),plot(x,y),grid on
ylabel('Dane')
subplot(2,1,2),plot(x,y1),grid on
ylabel('Pochodna')
y
x3
x 5x 2x 1
3
2
Slide 30
Całkowanie – całki
oznaczone.
Funkcja
b
f ( x ) dx
a
Opis
Q=quad(f,a,b)
Q=quad(f,a,b,tol)
Q=quad(f,a,b,tol,tr)
oblicza wyrażenie za pomocą prostej
kwadratury Simpsona
Q=quad8(f,a,b,tol,tr)
oblicza wyrażenie za pomocą złożonej
kwadratury Newtona rzędu 8; daje
lepsze wyniki niż poprzednia
Parametry:
f – łańcuch znaków określający nazwę funkcji;
a, b – przedział całkowania;
tol – tolerancja błędu (domyślnie 10-3);
tr – jeśli jest różny od 0, rysowany jest wykres obrazujący postęp w
obliczeniach lub wypisywane kolejne iteracje.
Slide 31
Całkowanie
Przykład:
>> quad('sin(x)',1,10)
ans =
1.3794
>> quad8('sin(x)',1,10)
ans =
1.3794
Gdy dodamy kolejne parametry:
>> quad('sin(x)',1,10,1e-3,1)
5
7
9
11
13
15
17
19
21
23
25
ans =
1.0000000000
1.0000000000
1.0000000000
1.0000000000
2.1250000000
3.2500000000
3.2500000000
4.3750000000
5.5000000000
5.5000000000
7.7500000000
1.3793
9.00000000e+000 2.1980581837
4.50000000e+000 -0.1666251940
2.25000000e+000 1.5343093586
1.12500000e+000 1.0665674753
1.12500000e+000 0.4678628302
2.25000000e+000 -1.7026633720
1.12500000e+000 -0.6631045444
1.12500000e+000 -1.0396930457
4.50000000e+000 1.5317252037
2.25000000e+000 0.6048270787
2.25000000e+000 0.9427905376
10
sin(
1
x ) dx
Slide 32
Całkowanie
Przykład:
W przypadku funkcji quad8 otrzymamy wykres:
>> quad8('sin(x)',1,10,1e-3,1)
ans =
1.3794
10
sin(
1
x ) dx
Slide 33
Równania różniczkowe
Równania różniczkowe możemy podzielić na 3 grupy:
- równania różniczkowe zwyczajne gdzie szukamy rozwiązania
równania różniczkowego dla zadanego warunku początkowego;
- równania różniczkowe zwyczajne gdzie szukamy rozwiązania
równania dla zadanych warunków granicznych;
- równania różniczkowe cząstkowe
Slide 34
Równania różniczkowe I rzędu
W rozwiązaniach wielu problemów fizycznych, ekonomicznych
wymagana jest znajomość funkcji y=y(t) przy znajomości funkcji
y'=f(t, y) oraz warunków początkowych y(a)=y0,
gdzie: a oraz y0 są liczbami rzeczywistymi
f funkcją w postaci jawnej.
W takim przypadku mamy do czynienia z równaniem różniczkowym
zwyczajnym pierwszego rzędu
Slide 35
Slide 36
Przykład:
Szukamy rozwiązania numerycznych y = y(t) dla wartości
t = 0, .25, .5, .75, 1 dla równania różniczkowego y' = -2ty2,
przy warunku początkowym y(0) = 1.
dy = @(t,y)(-2*t.*y(1).^2);
tspan = [0 .25 .5 .75 1];
y0 = 1;
[t1 y1] = ode23('eq1', tspan, y0);
[t2 y2] = ode45('eq1', tspan, y0);
[t1 y1 y2]
% porównanie wyników
Wykład 6
https://play.google.com
Dr Aneta Polewko-Klim
Slide 2
Obliczenia numeryczne.
Wykonując obliczenia numeryczne zakładamy:
operujemy na liczbach, dopuszczając
możliwość występowania niedokładności w
ich reprezentacji,
podczas obliczeń mogą występować pewne
niedokładności.
Slide 3
Miejsca zerowe i minima funkcji.
Slide 4
Opcje
Wszystkie te funkcje można wywołać z opcjonalnym
parametrem opcje,
np.: x1 = fzero(‘cos(x) - 0.2’, 3, opcje)
Opcje określają parametry wywołania tych funkcji.
Ustawienia tych parametrów można zmienić poprzez
funkcję optimset.
opcje = optimset(parametr, wartość, ... )
Argumenty funkcji optimset to pary
parametr – wartość,
Parametr jest łańcuchem znakowym.
Slide 5
Funkcja optimset - lista parametrów i
ich możliwych wartości (kilka):
Slide 6
Wywołanie z nazwą funkcji
(wyświetla parametry dla konkretnej funkcji)
Slide 7
Przykłady parametrów funkcji
optimset:
Slide 8
Miejsca zerowe
Przykład:
>> fzero('cos(2*x)-.2',8)
ans =
8.7401
>> cos(2*ans)-.2
ans =
2.8866e-015
Slide 9
>> op=optimset('Diagnostics','on',...
'Display','Iter');
Looking for a zero in the interval [7.0949,
8.9051]
>> fzero('cos(2*x)-.2',8,op)
Func-count
x
f(x)
1
8
-1.15766
2
7.77373
-1.18715
3
8.22627 -0.935369
4
7.68
-1.14007
5
8.32
-0.7962
6
7.54745
-1.01789
7
8.45255 -0.565028
8
7.36 -0.750391
9
8.64
-0.19876
10
7.0949 -0.252615
11
8.9051
0.30677
Procedure
initial
search
search
search
search
search
search
search
search
search
search
12
7.91238
-1.19319
interpolation
13
8.70207 -0.0749546
interpolation
14
8.74193 0.00367162
interpolation
15
8.74007 2.49516e-005
interpolation
16
8.74006 -1.30267e-009
interpolation
17
8.74006 2.88658e-015
interpolation
18
8.74006 -4.05231e-015
interpolation
Zero found in the interval: [7.0949, 8.9051].
ans =
8.7401
Slide 10
Minimum – funkcja 1 zmiennej
Przykład
>>
funkcja=inline('sin(2*x)+cos(x)/2');
>>fplot(funkcja,[-2*pi, 2*pi])
>> x1=fminbnd(funkcja,-2,4)
x1 = 2.4375
>> funkcja(x1)
ans =
-1.3679
Slide 11
Minimum – funkcja 2 zmiennej
Do znalezienia minimum musimy
Podać postać tej funkcji lub
funkcję inline oraz punkt startowy.
Jak widać na rysunku funkcja
ta ma dwa minima. W zależności
od tego jaki punkt startowy
podamy, znajdziemy albo
pierwsze albo drugie.
Przykład
>> fun=inline('-exp(-(x(1)-2).^2-x(2).^2)-exp(-(x(1)+2).^2-x(2).^2)');
>> x0=fminsearch(fun,[2 2])
x0 = 2.0000 -0.0000
>> fx0=fun(x0)
fx0 = -1.0000
Slide 12
Pierwiastki wielomianów.
Ogólna postać wielomianu:
W(x) = a1xn + a2xn-1 + ... + anx + an+1
współczynniki wielomianu a1, a2, ... są uszeregowane według
malejących potęg zmiennej x.
Slide 13
Przykład:
W(x) = x2 + 3x – 4
>> roots([1 3 -4])
ans =
-4
1
>> poly(ans)
ans =
1 3 -4
W(x) = 3x5 – 2x4 + 4x3 – 6x + 1
>> roots([3 -2 4 0 -6 1])
ans =
0.2084 + 1.4463i
0.2084 - 1.4463i
-0.9198
1.0000
0.1697
Slide 14
Układy równań liniowych
Zapis macierzowy, każdy układ równań możemy przedstawić:
A*X = B
lub
X*A = B
Matlab dysponuje dwoma operatorami dzielenia, które odpowiednio
rozwiązują przypadki:
A*X = B mamy X = A\B – musi być równa liczba wierszy A i B
X*A = B mamy X = B/A – musi być równa liczba kolumn A i B
W technice częściej spotyka się pierwszy z tych przypadków.
Macierz A w ogólnym przypadku może mieć w wierszy i k kolumn.
Istnieją więc trzy możliwości:
w = k – macierz kwadratowa, jeśli jest nieosobliwa X można
wyznaczyć dokładnie;
w > k – układ jest nadokreślony;
w < k – układ jest nieoznaczony (rozwiązanie zawiera najwyżej w
niezerowych elementów;
Slide 15
Układ kwadratowy.
Mamy układ trzech równań
z trzema niewiadomymi.
Musimy zdefiniować
macierze A i B, a następnie
obliczyć wynik.
2 x1 3 x 2 4 x 3 5
x1 2 x 2 x 3 1
3x 2x x 5
1
2
3
Przykład
>> A=[2 -3 -4; -1 2 1 ; 3 -2 -1]
A=
2 -3 -4
-1 2 1
3 -2 -1
>> B=[-5 -1 5]'
B=
-5
-1
5
>> X=A\B
X=
2.0000
-1.0000
3.0000
Slide 16
Układ nieoznaczony cz.1
Liczba równań jest mniejsza
od liczby niewiadomych
6 x1 8 x 2 7 x 3 3 x 4 1
3 x1 5 x 2 4 x 3 x 4 2
Przykład
>> A=[6 8 7 3;3 5 4 1]
A=
6
8
7
3
5
4
>> B=[1 2]'
B=
1
2
>> X=A\B
X=
0
0.7143
0
-1.5714
3
1
Slide 17
Układ nieoznaczony cz.2
Ogólne rozwiązanie takiego
układu otrzymamy dodając
do otrzymanego rozwiązania
iloczyn macierzy Z będącej
jądrem macierzy A i
dowolnego wektora q.
Jądro: Z = null(A) A*Z = 0
Przykład cd.
>> Z=null(A)
Z=
-0.4253 -0.6519
-0.3981 0.3942
0.8126 -0.1607
0.0163 0.6276
>> q=[2 17]'
q= 2
17
>> X1=X+Z*q
X1 = -11.9324
6.6186
-1.1066
9.1306
Slide 18
Układ nadokreślony.
Gdy chcemy otrzymać n
współczynników funkcji
opisującej przebieg
jakiegoś zjawiska.
Dokonujemy wtedy np.
dużo więcej niż n
pomiarów i na ich
podstawie szukamy
wyniku.
x1 x 2 2 x 3 1
x1 2 x 2 x 3 2
3 x1 x 2 5 x 3 3
2 x1 2 x 2 3 x 3 4
Przykład
>> A=[1 -1 2; 1 -2 -1; 3 -1 5; -2 2 3]
A = 1 -1 2
1 -2 -1
3 -1 5
-2 2 3
>> B=[1 2 3 -4]'
B=1
2
3
-4
>> X=A\B
X = 1.4286
-0.1429
-0.2857
Slide 19
Interpolacja.
Zagadnienie dotyczy przechodzenia funkcji przez zadane punkty,
tzw. węzły interpolacji. Standardowe procedury Matlaba realizują
interpolację za pomocą wielomianów pierwszego i trzeciego
stopnia, oraz funkcjami sklejanymi trzeciego stopnia. Funkcje
interpolacyjne:
Slide 20
Metody interpolacji
‘nearest’— interpolacja najbliższy sąsiad
‘linear’ – interpolacja liniowa;
‘cubic’ – interpolacja wielomianem
trzeciego stopnia;
‘spline’ – interpolacja funkcjami sklejanymi
trzeciego stopnia;
Slide 21
Przykład:
x=0:.4:2;
y=exp(-.5*x).*cos(5*x);
xi=0:.05:2;
y1=interp1(x,y,xi,'nearest');
y2=interp1(x,y,xi,'linear');
y3=interp1(x,y,xi,'cubic');
y4=interp1(x,y,xi,'spline');
plot(x,y,'*',xi,exp(-.5*xi).*cos(5*xi),...
xi,y1,':',xi,y2,':',xi,y3,':',xi,y4,':')
legend('wezly','funkcja','nearest',…
'linear','cubic','spline')
Slide 22
Aproksymacja.
Zagadnienie aproksymacji polega na znalezieniu pewnej
funkcji, której wykres z pewnym przybliżeniem
przechodzi w pobliżu pewnych punktów.
Sama funkcja minimalizuje wartość pewnego kryterium
dopasowania, np. aproksymacja średniokwadratowa,
czyli minimalizujemy błąd średniokwadratowy.
Standardową metodą aproksymacji w Matlabie jest
aproksymacja wielomianem:
polyfit(x,y,n)
polecenie to znajduje wektor współczynników
wielomianu aproksymującego dane zawarte w wektorach
x i y. Wartość wielomianu możemy obliczyć funkcją
polyval.
Slide 23
Aproksymacja.
Przykład 1:
x=1:8;
y=[1.45 1.95 2.43 3.02 3.53 4 4.4 5.1];
a=polyfit(x,y,1)
xx=1:.01:8;
yy=polyval(a,xx);
plot(x,y,'o,xx,yy)
Otrzymaliśmy funkcję y=0.5121x+0.9304
Slide 24
Przykład 2:
x=0:.4:2*pi;
y1=sin(x);
blad=rand(1,length(x))/10;
y=y1+blad;
a1=polyfit(x,y,2)
a2=polyfit(x,y,3)
xx=0:.01:2*pi;
yy1=polyval(a1,xx);
yy2=polyval(a2,xx);
plot(x,y,'o',xx,sin(xx),':',xx,yy1,xx,yy2)
Dla wielomianu drugiego stopnia mamy:
y = – 0.0323x2 – 0.0982x + 0.7252
Dla wielomianu trzeciego stopnia mamy:
y = 0.0921x3 – 0.8611 x2 + 1.8274x – 0.0793
Slide 25
Pochodna funkcji
• Do znajdowania pochodnej możemy posłużyć się
funkcją:
diff
Zwraca ona wektor różnic sąsiednich elementów wektora
będącego jej argumentem.
• Czyli pochodna funkcji danej w postaci tablicy (np.
pomiary), obliczymy jako:
diff(y)./diff(x).
• Należy przy tym pamiętać, że liczba elementów
wektora zwracanego przez funkcję diff, jest o jeden
mniejsza, niż liczba elementów argumentu.
Slide 26
Pochodna funkcji
Przykład
x=72:.5:90;
y=[7237 7224 7203 7183 7164 7149 7129 7107 7060 7011 6912
6746 6517 ...
6225 5875 5449 4950 4376 3789 3202 2614 2082 1639 1290
1000 ...
758 581 470 394 328 284 241 191 146 103 52 0];
yprim=diff(y)./diff(x);
subplot(2,1,1),plot(x,y,'+'),grid on
ylabel('Dane')
subplot(2,1,2),plot(x, [yprim 0],'o'),grid on
ylabel('Pochodna')
Slide 27
Pochodna wielomianu
• Jeśli mamy dany wielomian, do obliczenia wartości
pochodnej możemy użyć funkcji:
polyder
• Zwraca ona współczynniki wielomianu pochodnego.
Wywołanie:
a_prim = polyder(a)
[L_prim, M_prim] = polyder(L,M)
Slide 28
Pochodna wielomianu
a=[-1 3 5 -2];
x=0:.05:3;
y=polyval(a,x);
a1=polyder(a);
y1=polyval(a1,x);
subplot(2,1,1);
plot(x,y)
grid on
ylabel('Dane')
subplot(2,1,2)
plot(x,y1)
grid on
ylabel('Pochodna')
Slide 29
Pochodna wielomianu - funkcja wymierna
Dodatkowo należy musimy zdefiniować licznik i mianownik:
Przykład:
L=[1 -3];
M=[1 5 -2 1];
x=-3:.05:3;
y=polyval(L,x)./polyval(M,x);
[L1,M1]=polyder(L,M);
y1=polyval(L1,x)./polyval(M1,x);
subplot(2,1,1),plot(x,y),grid on
ylabel('Dane')
subplot(2,1,2),plot(x,y1),grid on
ylabel('Pochodna')
y
x3
x 5x 2x 1
3
2
Slide 30
Całkowanie – całki
oznaczone.
Funkcja
b
f ( x ) dx
a
Opis
Q=quad(f,a,b)
Q=quad(f,a,b,tol)
Q=quad(f,a,b,tol,tr)
oblicza wyrażenie za pomocą prostej
kwadratury Simpsona
Q=quad8(f,a,b,tol,tr)
oblicza wyrażenie za pomocą złożonej
kwadratury Newtona rzędu 8; daje
lepsze wyniki niż poprzednia
Parametry:
f – łańcuch znaków określający nazwę funkcji;
a, b – przedział całkowania;
tol – tolerancja błędu (domyślnie 10-3);
tr – jeśli jest różny od 0, rysowany jest wykres obrazujący postęp w
obliczeniach lub wypisywane kolejne iteracje.
Slide 31
Całkowanie
Przykład:
>> quad('sin(x)',1,10)
ans =
1.3794
>> quad8('sin(x)',1,10)
ans =
1.3794
Gdy dodamy kolejne parametry:
>> quad('sin(x)',1,10,1e-3,1)
5
7
9
11
13
15
17
19
21
23
25
ans =
1.0000000000
1.0000000000
1.0000000000
1.0000000000
2.1250000000
3.2500000000
3.2500000000
4.3750000000
5.5000000000
5.5000000000
7.7500000000
1.3793
9.00000000e+000 2.1980581837
4.50000000e+000 -0.1666251940
2.25000000e+000 1.5343093586
1.12500000e+000 1.0665674753
1.12500000e+000 0.4678628302
2.25000000e+000 -1.7026633720
1.12500000e+000 -0.6631045444
1.12500000e+000 -1.0396930457
4.50000000e+000 1.5317252037
2.25000000e+000 0.6048270787
2.25000000e+000 0.9427905376
10
sin(
1
x ) dx
Slide 32
Całkowanie
Przykład:
W przypadku funkcji quad8 otrzymamy wykres:
>> quad8('sin(x)',1,10,1e-3,1)
ans =
1.3794
10
sin(
1
x ) dx
Slide 33
Równania różniczkowe
Równania różniczkowe możemy podzielić na 3 grupy:
- równania różniczkowe zwyczajne gdzie szukamy rozwiązania
równania różniczkowego dla zadanego warunku początkowego;
- równania różniczkowe zwyczajne gdzie szukamy rozwiązania
równania dla zadanych warunków granicznych;
- równania różniczkowe cząstkowe
Slide 34
Równania różniczkowe I rzędu
W rozwiązaniach wielu problemów fizycznych, ekonomicznych
wymagana jest znajomość funkcji y=y(t) przy znajomości funkcji
y'=f(t, y) oraz warunków początkowych y(a)=y0,
gdzie: a oraz y0 są liczbami rzeczywistymi
f funkcją w postaci jawnej.
W takim przypadku mamy do czynienia z równaniem różniczkowym
zwyczajnym pierwszego rzędu
Slide 35
Slide 36
Przykład:
Szukamy rozwiązania numerycznych y = y(t) dla wartości
t = 0, .25, .5, .75, 1 dla równania różniczkowego y' = -2ty2,
przy warunku początkowym y(0) = 1.
dy = @(t,y)(-2*t.*y(1).^2);
tspan = [0 .25 .5 .75 1];
y0 = 1;
[t1 y1] = ode23('eq1', tspan, y0);
[t2 y2] = ode45('eq1', tspan, y0);
[t1 y1 y2]
% porównanie wyników