Slad 1
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Całkowanie numeryczne
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 2
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Plan zajęć
1. Całkowanie przybliżone funkcji jednej zmiennej z użyciem
wbudowanych funkcji Scilaba integrate(), inttrap(),
intsplin()
2. Całkowanie funkcji 2 zmiennych
"
wykorzystanie wbudowanej funkcji int2d()
"
użycie metody Monte Carlo
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 3
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1
Utwórz skrypt (zapisując go jako ~/calki.sce) zawierający instrukcje wykonujące
następujące czynności:

f (x) = 1- x2
zdefiniowanie funkcji użytkownika

narysowanie wykresu funkcji na przedziale [-1,1]

przy użyciu funkcji Scilaba integrate() obliczenie wartość całki
1
Ä„
(= )
+"f (x)dx
4
0
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 4
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1
Utwórz skrypt (zapisując go jako ~/calki.sce) zawierający instrukcje wykonujące
następujące czynności:

f (x) = 1- x2
zdefiniowanie funkcji użytkownika

narysowanie wykresu funkcji na przedziale [-1,1]

przy użyciu funkcji Scilaba integrate() obliczenie wartość całki
1
Ä„
(= )
+"f (x)dx
4
1.0
0
0.9
0.8
0.7
// zdefiniowanie funkcji f(x), narysowanie wykresu
0.6
deff('[y]=f(x)','y=sqrt(1-x^2)')
0.5
xp = -1:0.05:1
0.4
plot2d (xp, f(xp), [2]) 0.3
0.2
0.1
// obliczenie całki
0.0
w = integrate('f(x)','x',0,1) -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
w =
0.7853982
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 5
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1 - cd
obliczenie wartości całki przy użyciu funkcji inttrap(), intsplin():

otrzymując wyniki przybliżonych wartości całki dla każdej z funkcji dokonując
podziału przedziału całkowania na 2k części (k=1,2,...,10). Wyniki zapisz w
dwóch wektorach wt[1:10], ws[1:10]

otrzymane wyniki przedstaw na wykresie w zależności od zmiennej
k=1,2,...,10.
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 6
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1 - cd
obliczenie wartości całki przy użyciu funkcji inttrap(), intsplin():

otrzymując wyniki przybliżonych wartości całki dla każdej z funkcji dokonując
podziału przedziału całkowania na 2k części (k=1,2,...,10). Wyniki zapisz w
dwóch wektorach wt[1:10], ws[1:10]

otrzymane wyniki przedstaw na wykresie w zależności od zmiennej
k=1,2,...,10.
// obliczenie całki przy użyciu funkcji inttrap(), intsplin()
for i=1:10
x=linspace(0,1,2^i)
0 .8 0
wt(i)=inttrap(x,f(x))
0 .7 5
ws(i)=intsplin(x,f(x))
end
0 .7 0
// narysowanie wykresu
0 .6 5
k=1:10, plot2d(k, [wt(k) ws(k)])
0 .6 0
0 .5 5
0 .5 0
1 2 3 4 5 6 7 8 9 1 0
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 7
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1 - cd

otrzymane wyniki zapisz do sformatowanego pliku tekstowego ~/w-calki.txt
zapisując w 3 kolumnach (podając wartości całek z dokładnością do 6 miejsc po przecinku):

k,

wynik otrzymany za pomocÄ… inttrap(),

wynik otrzymany za pomocÄ… intsplin().
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 8
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 1 - cd

otrzymane wyniki zapisz do sformatowanego pliku tekstowego ~/w-calki.txt
zapisując w 3 kolumnach (podając wartości całek z dokładnością do 6 miejsc po przecinku):

k,

wynik otrzymany za pomocÄ… inttrap(),

wynik otrzymany za pomocÄ… intsplin().
// zapisanie wyników do pliku tekstowego
u1 = file('open','~/w-calki.txt','new')
for k=1:10
write(u1,[k, wt(k), ws(k)],'(f3.0,2X,f9.6,2X,f9.6)')
end
file('close',u1);
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 9
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2
W otwartym skrypcie dopisz instrukcje służące do obliczenia,
korzystając z funkcji int2d() wartości przybliżone całek dla :
+"fdS
+"zdS
+"fdH
S
S
H
gdzie

z(x,y)=1, S=[0,2]x[0,4] (wartość dokładna = 8)

f(x,y)=xy2, S=[0,2]x[0,4] (wartość dokładna = 42,666)

f(x,y)=xy2 , H= jest obszarem pomiędzy prostą y=2x a parabolą y=x2
(wartość dokładna = 32/5)
funkcja I = int2d(X,Y,f) oblicza całkę po obszarze podzielonym na trójkąty

X macierz prostokątna wymiaru 3 x N - każda kolumna przedstawia
współrzędne X-owe jednego trójkąta

Y macierz prostokątna wymiaru 3 x N - każda kolumna przedstawia
współrzędne Y-owe jednego trójkąta

f zdefiniowana funkcja podcałkowa
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 10
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla z(x,y)=1, S=[0,2]x[0,4] - dokonując podziału obszaru całkowania na
2 trójkąty
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 11
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla z(x,y)=1, S=[0,2]x[0,4] - dokonując podziału obszaru całkowania na
2 trójkąty
// podział obszaru całkowania [0,2]x[0,4] na 2
trójkąty
X=[0,0;2,0;2,2];
Y=[0,0;0,4;4,4];
// zdefiniowanie funkcji podcałkowej
deff('zp=z(x,y)','zp=1')
// obliczenie całki
I=int2d(X,Y,z)
disp(I)
> 8.
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 12
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla f(x,y)=xy2, S=[0,2]x[0,4], wykorzystując już utworzony podział

Dla obszaru całkowania pomiędzy prostą y=2x a parabolą y=x2 wpisując obszar w
prostokąt [0,2]x[0,4], definiując nową funkcję podcałkową i wykorzystując już
utworzony podział
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 13
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla z(x,y)=xy2, S=[0,2]x[0,4], wykorzystując już utworzony podział

Dla obszaru całkowania pomiędzy prostą y=2x a parabolą y=x2 wpisując obszar w
prostokąt [0,2]x[0,4], definiując nową funkcję podcałkową i wykorzystując już
utworzony podział
// zdefiniowanie funkcji podcałkowej
deff('z=f(x,y)','z=x*y^2')
I=int2d(X,Y,f) // obliczenie całki
disp(I)
> 42.666667
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 14
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla z(x,y)=xy2, S=[0,2]x[0,4], wykorzystując już utworzony podział

Dla obszaru całkowania pomiędzy prostą y=2x a parabolą y=x2 wpisując obszar w
prostokąt [0,2]x[0,4], definiując nową funkcję podcałkową i wykorzystując już
utworzony podział
// zdefiniowanie funkcji podcałkowej
deff('z=f(x,y)','z=x*y^2')
I=int2d(X,Y,f) // obliczenie całki
disp(I)
> 42.666667
function [z]=fn(x,y)
z=0
if (y<=2*x)&(y>=x^2) then
z=f(x,y)
end
endfunction
I=int2d(X,Y,fn) // obliczenie całki
disp(I)
> 4.8634489 (wartość dokładna 32/5 = 6.4)
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 15
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla obszaru całkowania pomiędzy prostą y=2x a parabolą y=x2 wpisując obszar w
prostokąt [0,2]x[0,4], dzieląc prostokąt na trójkąty prostokątne o
przyprostokÄ…tnych = 1
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 16
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Dla obszaru całkowania pomiędzy prostą y=2x a parabolą y=x2 wpisując obszar w
prostokąt [0,2]x[0,4], dzieląc prostokąt na trójkąty prostokątne o
przyprostokÄ…tnych = 1
// podział obszaru [0,2]x[0,4] na trójkąty prostokątne o boku =1
for i=1:2
for j=1:4
k=2*i-1+4*(j-1)
X(1,k)=i-1, X(2,k)=i, X(3,k)=i
Y(1,k)=j-1, Y(2,k)=j-1, Y(3,k)=j
k=2*i + 4*(j-1)
X(1,k)=i-1, X(2,k)=i-1, X(3,k)=i
Y(1,k)=j-1, Y(2,k)=j, Y(3,k)=j
end
end
I=int2d(X,Y,fn), disp(I)
> 6.5729887 (wartość dokładna 32/5 = 6.4)
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 17
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd
Oszacować przy użyciu Metody Monte Carlo, korzystając z generatora
liczb losowych wartość całki
+"zdS dla funkcji :
S
z(x,y)=xy2 , S jest obszarem pomiędzy prostą y=2x a parabolą y=x2 (wartość dokładna = 32/5)

znajdując N0 =10000 punktów losowych w obszarze S0 = obszar prostokątny [0,2] x [0,4]
// wygenerowanie punktów losowych z obszaru [0,2] x [0,4]
N_0=1000
xy =rand(2,N_0)
xy(1,:)=2*xy(1,:), xy(2,:)=4*xy(2,:)
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 18
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Obliczając N - liczbę punktów losowych leżących w obszarze S i korzystając z formuły
przybliżyć wartość całki (suma we wzorze obliczana jest tylko dla punktów leżących w
obszarze całkowania S)
N N
1 N 1 S0 N
z(xi , yi ) H" S0 Å" z(xi , yi ) = z(xi , yi )
" " "
+"zdS H" S N
N0 N N0 i= 1
i= 1 i= 1
S
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 19
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Ćwiczenie 2 - cd

Obliczając N - liczbę punktów losowych leżących w obszarze S i korzystając z formuły
przybliżyć wartość całki (suma we wzorze obliczana jest tylko dla punktów leżących w
obszarze całkowania S)
N N
1 N 1 S0 N
z(xi , yi ) H" S0 Å" z(xi , yi ) = z(xi , yi )
" " "
+"zdS H" S N
N0 N N0 i= 1
i= 1 i= 1
S
N=0, sumka = 0, S0 = 8
for i=1:N_0
x= xy(1,i), y= xy(2,i)
// sprawdzenie czy losowy punkt należy do obszaru S
if (y<=2*x)&(y>=x^2) then
sumka = sumka + f(x,y)
N = N+1
end
end
calka = S0 * sumka/N_0
disp(calka), disp(N),
> 6.2687125 (wartość całki - dokładna 6.4)
> 172
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.
Slad 20
j
27 styczeń 2009
SciLab w obliczeniach numerycznych - część 4
Funkcje SciLaba
Funkcje przybliżonego całkowania:

int2d()  obliczenie całki z funkcji 2 zmiennych po obszarze opisanym siatką
trójkątów

int3d()  obliczenie całki z funkcji 3 zmiennych, obszar całkowania opisany
siatką czworościanów

integrate(), intg()  obliczenie całki z funkcji jednej zmiennej metodą
kwadratur

intsplin()  obliczenie całki z funkcji sklejanej (jednej zmiennej) interpolującej
zbiór punktów

inttrap()  obliczenie całki z funkcji (jednej zmiennej) interpolującej zbiór
punktów  wzór trapezów
Instytut K onstrukcji Budowlanych Metody Obliczeniowe
Zakład K omputerowego Wspomagania Projektowania Semestr II.