I. PREZENTAREA GENERALĂ GENERALĂ A PROGRAMULUI SCILAB
I. 1. Ce este ScilabScilab-ul? Structural, SCILAB-ul este, cu mici modificări, o variantă gratuită a MATLAB-ului având versiuni pentru Unix şi pentru Windows. La adresa www.scilab.org se găseşte documentaţie gratuită pe INTERNET referitoare la Scilab precum şi programul corespunzător. Pentru început, Scilab-ul se poate folosi pentru efectuarea calculelor elementare, operaţiilor cu vectori şi matrice (cu elemente reale şi complexe) şi reprezentărilor grafice de curbe şi suprafeţe. De asemenea, Scilabul poate fi un mediu agreabil de realizare a calculelor numerice pentru că dispune de o serie de metode uzuale din acest domeniu, de exemplu: - rezolvarea sistemelor liniare, - calculul valorilor şi vectorilor proprii, - rezolvarea ecuaţiilor diferenţiale, - rezolvarea sistemelor neliniare, - noţiuni de bază din algebra liniară, - algoritmi de optimizare Programarea în Scilab se efectuează destul de uşor, pentru anumite secvenţe mai complicate se pot folosi subprograme realizate în C sau fortran (77). Scilab-ul include aplicaţii specifice, numite TOOLBOX-uri - colecţii extinse de funcţii SCILAB care dezvoltă mediul de programare de la o versiune la alta, pentru a rezolva probleme din domenii variate.
I. 2. Mediul Scilab Mediul Scilab lucrează cu mai multe tipuri de ferestre: o fereastră de comenzi, o fereastră pentru fişiere sau grafice, o fereastră pentru ‘help’. Când se lansează în execuţie se deschide implicit fereastra de comandă care, conţine meniurile File, Edit, Preferences, Control, Editor, Applications şi ? (help) şi o bară de butoane unde Buton
Semnificaţie Deschide o nouă fereastră de comandă Deschide Scipad-ul (editorul Scilab-ului) Deschide un fişier Copiează Lipeşte Schimbă directorul Consolă Scilab Alege font
În
fereastra
de
comenzi,
instrucţiunile
se
scriu
pe
prompterul --> şi se execută apăsând Enter; instrucţiunile se execută de la stanga la dreapta şi de sus în jos în ordinea în care s-au scris.
Meniul File conţine:
unde, un fişier specificat este executat cu Exec şi încărcat cu Load. Load Activat meniul Editor face să se dechidă editorul Scilab-ului în care se pot scrie programe care se pot executa folosind Execute/Load into Scilab. Scilab
Folosind Scicos din Applications se pot realiza diverse diagrame şi simulări în fereastra care se deschide.
Pentru realizarea diagramelor se pot alege obiecte folosind opţiunea Palettes din Edit (sau clic dreapta cu mouse-ul), de exemplu:
Se pot introduce instrucţiuni cu ajutorul ferestrei de dialog obţinute prin accesarea opţiunii Context (clic dreapta cu mouseul).
Într-o fereastră asemănătoare se pot edita şi personaliza graficele.
II. INTRODUCERE ÎN SCILAB SCILAB
II. 1. Noţiuni Scilab Scilabul este un pachet de programe de înaltă performanţă, dedicat calcului numeric şi reprezentărilor grafice în domeniul ştiinţei şi
ingineriei,
precum
şi
simulărilor
stocastice.
Matricea
este
elementul de bază; cu aceasta se pot rezolva probleme fără a fi necesară scrierea unui program într-un limbaj de programare. Tabel 1. Funcţii SCILAB de interes general: Denumire funcţie exit, quit help who cd dir clc error exist disp clear $ apropos type(a) pause
Acţiune funcţie Comanda pentru ieşirea din SCILAB Furnizează informaţii “on-line” despre SCILAB Listează variabilele curente Schimbă directorul de lucru Afişează lista directorului Şterge fereastra de comenzi Afişează mesaje de eroare Verifică dăcă o variabilă sau un fişier există Afişează o matrice sau un text Şterge variabile şi funcţii Indexul ultimului element al unei matrice sau a unui vector Caută în ‘help’ un cuvânt specificat Tipul variabilei a Intrerupe executia unei functii sau a unui script
Tabel 2. Operatorii aritmetici SCILAB:
Operaţia Adunare Scădere Înmulţire Împărţire la stânga Împărţire la dreapta Ridicare la putere Transpunere
Scalari
Matrice Tablouri
+ * \ / ^
+ * \ / ^ ‘
+ .* .\ ./ .^ .’
Priorita tea 4 4 3 3 3 2 1
Tabel 3. Operatorii relaţionali SCILAB: Operator < > <= >= == ~=
Semnificaţie Mai mic Mai mare Mai mic sau egal Mai mare sau egal Identic Diferit
Observatie: Observatie Operatorii aritmetici au prioritate fata de operatorii relaţionali. Tabel 4. Operatorii Operatorii logici SCILAB: Operator NU ŞI SAU
Simbol SCILAB ~ & |
Prioritatea 1 2 3
Tabel 5. Constante uzuale SCILAB:
Constanta
%i %pi %e %eps %inf %s %t (%f)
Semnificaţiea Numarul imaginar
−1
π = 3.1415926...
Constanta lui Euler e=2.7182818... eps ≈ 2.2 ⋅ 10−16 ∞
Expresie polinomiala s=poly(0,’s’) Variabila booleana true (false)
Tabel 6. Variabile SCILAB:
Nr.crt. 1.
2.
3.
4.
5.
Tipul variabilei
Exemplu
Numeric - scalar si a=2.5; b=2.5e3 matrice A=[1 2 3;4 5 6] (real sau complex) Sir de caractere -->sir="Acesta este un sir de caractere"; -->length(sir) ans = 31. Boolean -->2*%t ans = 2. -->7*%f ans = 0. Listă -->l=list(1,2,3,4,5) (list, mlist, tlist) l = l(1) 1. l(2) 2. l(3) 3. l(4) 4. l(5) 5. Polinomial -->p=poly([1 2 3 4],'x','coeff') p = 2 3 1 + 2x + 3x + 4x
Tabel 7. Funcţii uzuale SCILAB:
Apelarea funcţiei
Semnificaţia matematică
sin(x) sinc(x)
sin x
cos(x) tan(x) cotg(x) sinh(x) cosh(x) tanh(x) asin(x) acos(x) atan(x) exp(x) log(x), log10(x) floor(x) ceil(x) Re(z), Im(z) lenght(v) size(a) max(v) min(v) sqrt(x)
cos x tg x ctg(x) sh x ch x th x arcsin x arccos x arctg x ex ln(x), lg(x) [x] [x]+1 Re z, Im z Numărul elem. vectorului v Dimensiunea matricei a Cel mai mare elem. al vect. v Cel mai mic elem. al vect. v
sin( x) x
x
Exemple 1. Să se calculeze: 5
11 π a) − 2 + ln(5) − lg(21) + 2 4
SoluŃie: -->(%pi/4)^5-sqrt(2)+log(5)-log10(21)+2^11 ans =
2047.1719
b) sin(x)sin(x)-tg2(x)+2arctg(3x)(x)+2arctg(3x)-ch(x) pentru x dat.
2. Fie numerele complexe z1=2+5i si z2=z2=-3+i, să se calculeze:
z1+z2, z1*z2, z1/z2, z17. SoluŃie: -->z1=2+5*%i z1 =
2. + 5.i
-->z2=-3+%i z2 = - 3. + i -->z1+z2 ans =
-1. + 6.i
3. Să se calculeze valoarea expresiilor urmatoare pentru un x dat: x3 + 1 ; 5
a) 2 sin( π ⋅ x ) −
SoluŃie: -->x=2; 2*sin(%pi*x)-sqrt((x^3+1)/5) ans = - 1.3416408 3
1 ; x − 1
b) 2 x 2 +
1 c) 2arctg( x ) + ln x − 1
3
2
SoluŃie: -->x=5; 2*atan(x^2)+log(1/(x-1)^3) ans = - 1.0972478 d)
x 5 + 8 − lg( x − 1) ;
1 e) 2 arcsin( π ⋅ x ) + 2 x − 1 2
2
1 ; x − 1
f) 7 sin x 2 + log 5
(
)
π⋅ x 2 − ch x + 2 ; 5
g) tg
2
h) P = 3 x 2 + 2 x + 1 ;
i) sh(2 x + ln( x ))8 ;
3
Exemplu listă: listă multe obiecte în Scilab se introduc cu funcţiile tlist sau mlist; mlist de exemplu, dacă se defineşte o variabilă cu un anumit număr de câmpuri, aceasta se realizează uşor utilizând operatorul „..”. Variabila x de tip tlist se defineşte uşor astfel: -->x.culoare=4; -->x.valoare=rand(1,3); -->x.nume='un_nume'; -->x x = culoare: 4 valoare: [0.8784126,0.1138360,0.1998338] nume: "un_nume"
II. 2. Manevrarea vectorilor şi matricelor
Definirea matricelor: Introducerea explicită a listei de elemente; Generarea prin instrucţiuni şi funcţii; Crearea de fişiere script; Încărcarea din fişiere de date externe.Introducerea explicită a
unei matrice: Elementele unei linii sunt separate prin virgulă sau spaţiu; Lista elementelor matricei trebuie încadrată între paranteze drepte; Liniile se separă prin punct şi virgulă (‘;’). Exemplu pentru definirea unei matrice: matrice
-->A=[1 1 1;2 4 8;3 9 27] Exemplu pentru definirea unui vector: vector
-->b=[ 2 10 44 190] b = ! 2.
10.
44.
190. !
Tabel 8. Funcţii pentru construirea unui vector sau unei matrice tip Funcţia eye
diag
zeros
Semnificaţia Matricea unitate
Exemplu
I=eye(4,4) I = ! 1. 0. 0. ! ! 0. 1. 0. ! ! 0. 0. 1. ! Matrice B=diag(b) diagonală B = ! 2. 0. 0. ! 0. 10. 0. ! 0. 0. 44. ! 0. 0. 0. Matrice cu -->C=zeros(3,2) toate C = elementele 0 ! 0. 0. !
0. ! 0. ! 0. ! 190. !
ones
Matrice cu toate elementele 1
rand
Matrice cu elemente aleatoare
triu
Partea triunghiulară superioară a unei matrice
tril
Partea triunghiulară inferioară a unei matrice
linspace
Vector cu incrementare constantă între doua valori
poly
Matrice polinomiala
! 0. 0. ! ! 0. 0. ! -->D=ones(2,3) D = ! 1. 1. 1. ! ! 1. 1. 1. ! -->M=rand(3,4) M = ! 0.2113249 0.3303271 0.8497452 0.0683740 ! ! 0.7560439 0.6653811 0.6857310 0.5608486 ! ! 0.0002211 0.6283918 0.8782165 0.6623569 ! -->U=triu(D) U = ! 1. 1. 1. ! ! 0. 1. 1. ! ! 0. 0. 1. ! -->V=tril(D) V = ! 1. 0. 0. ! ! 1. 1. 0. ! ! 1. 1. 1. ! -->v=linspace(0,1,5) v = ! 0. 0.25 0.5 0.75 1. ! Observaţie: acelas lucru se obtine daca -->x=0:0.25:1 x = ! 0. 0.25 0.5 0.75 1. ! Pentru -->y=1:5 se obtine y =! 1. 2. 3. 4. 5. ! Definirea unui polinom prin radacini: -->p=poly([1 3],'s') p = 3 - 4s + s2 Definirea unui polinom prin coeficienti:
-->poly([1 2],'s','c') ans = 1 + 2s Exemple şi exrciţii
1 2 si b = (5 6) , apoi săa se a = 1. Să se definească matricele 3 4 1 2 construiască matricea c = 3 4 folosind matricele 5 6
a si b.
SoluŃie: -->a=[1 2;3 4];b=[5 6]; -->c=[a;b] c =! 1.
2. !
! 3.
4. !
! 5.
6. !
2. Să se construiască matricea d prin extragerea liniilor 2 şi 3 din matricea c păstrând toate coloanele (prima coloană) a matricei c.
SoluŃie: -->d=c(2:3,:) d =! 3.
4. !
! 5.
6. !
-->d=c(2:3,1) d =! 3. ! ! 5. !
3. Să se calculeze suma elementelor matricei a.
SoluŃie: -->sum(a) ans =
10.
4. Să se calculeze a5. 5. Să se calculeze aT. 6. Să se calculeze a-1. 7. Să se calculeze urma matricei a. 8. Să se calculeze rangul matricei a. 9. Fie A o matrice patratică; ce reprezintă diag(diag(A))?
10.
Fie x o matrice (vector), scrieŃi instrucŃiunea care calculează
matricea y a carui element (i,j) este egal cu f(xi,j) unde: a) f(x) = 2x2-3x+1;
SoluŃie: -->x=[2 -5 7 4 -1]; -->2*x^2-3*x+1 ans =! 3.
66.
78.
21.
6. !
b) f ( x ) = 2 x 2 − 3 x + 1 ; c) f(x) = (x(x-1)/(x+4);
SoluŃie: -->(x-1)./(x+4) ans =! 0.1666667
6.
0.5454545
0.375 - 0.6666667 !
Observaţie: Observaţie atenţie la operaţiile element cu element când se lucrează cu un vector sau o matrice; trebuie adăugat un punct în faţa operatorilor de multiplicare, *, impărţire, /, sau ridicare la putere, ^. Exemplul 1: 1 la exerciţiul de mai sus dacă nu se pune punctul înaitea semnului / scilabul dă un raspuns eronat: -->(x-1)/(x+4) ans =
0.4155844
Exemplul 2: 2 Dacă se definesc două matrice pătratice, A şi B, observaţi diferenţa dintre A*B (se înmulţeşte matricea A cu matricea B) şi A .*B (se înmulţeşte termen cu termen A(i,j)*B(i,j)) -->A=[1 2 3;2 3 4;3 4 5]; -->B=[5 6 7;6 7 8;7 8 9]; -->A*B ans = 38.
44.
50.
56.
65.
74.
74.
86.
98.
-->A .*B
ans = 5.
12.
21.
12.
21.
32.
21.
32.
45.
d) f ( x ) = 11.
1 1 + x2
Să se definească o matrice pătratică de ordinul 5, a; să se indice
matricele triunghiulare corespunzatoare acesteia respectiv inferioară).
(superioară,
II. 3. Programare în Scilab
SCILAB-ul poate lucra în urmatoarele moduri:
- În modul linie de comandă; fiecare linie este prelucrată imediat şi rezultatele sunt afişate; - Cu programe conţinute în fişiere. Un fişier constă dintr-o succesiune de instrucţiuni SCILAB, cu posibilitatea apelării altor fişiere precum şi a apelării recursive. Un program SCILAB poate fi scris sub forma fişierelor:
- Script: Script conţine o secventă de comenzi SCILAB; se execută în fereastra
de
comandă
prin
apelarea
numelui
fişierului:
exec(‘nume_fisier’) sau direct, în editorul Scilab-ului folosind din meniul ‘Execute’ opţiunea ‘Load into Scilab’.Funcţie Funcţie: Funcţie prima linie conţine cuvântul function, iar ultima linie conţine cuvântul
endfunction; poate lucra cu argumente. Se pot defini mai multe funcţii în acelaşi fişier.
Tabel 9. Instrucţiuni şi funcţii de control Instrucţiune if else elseif case for while break return error end
Semnificaţia Instrucţiune pentru execuţia condiţională Instrucţiune asociată cu if Instrucţiune asociată cu if Instrucţiune folosită pentru a face o alegere Instrucţiune pentru crearea ciclurilor cu număr specificat de paşi Instrucţiune pentru crearea ciclurilor cu condiţie logică Instrucţiune pentru terminarea forţată întrun ciclu Returnează execuţia la funcţia precizată Instrucţiune pentru afişarea unui mesaj de eroare Instrucţiune pentru încheierea ciclurilor for,
while şi if Tabel 10. Sintaxa instrucţiunilor SCILAB Instrucţiune
Sintaxă if expresie_logică grup instrucţiuni end if expresie_logică grup instrucţiuni A else grup instrucţiuni B end if expresie_logică_1 grup instrucţiuni A elseif expresie_logică_2 grup instrucţiuni B end select expr_0 case expr_1 Instr_1 … case expr_n Instr_n else Instr end for index = expresie \\ expresie = iniţial: pas: final grup_instrucţiuni end while expresie grup_instrucţiuni end
if else elseif
case
for
while
II. 3.1. Fişiere de tip script
Definirea şi reprezentarea grafică a funcţiilor cu acoladă: cu ajutorul Scipad-ului (editorul Scilab-ului) se editează un fişier script care se salvează, de exemplu, Functie.sce şi conţine următoarele instrucţiuni:
a=-4;b=4;d=0.01;c=0; x=a:d:b; k=max(size(x)); for i=1:k if (x(i)>=a)&(x(i)=c)&(x(i)<=%pi)
y(i)=sin(x(i)); elseif (x(i)>%pi)&(x(i)<=b) y(i)=(x(i))^1/3-%pi^1/3; end end plot(x,y) În fereastra de comandă se execută fişierul folosind instrucţiunea:
exec('functie.sce') şi se va obţine într-o nouă fereastră graficul funcţiei: x + sin( x ) daca − 4 ≤ x ≤ 0 f ( x ) = sin( x ) daca 0 ≤ x ≤ π 3 3 x − π daca π ≤ x ≤ 4
(
)
1
0
-1
-2
-3
-4 -4
-3
-2
-1
0
1
2
3
4
Observaţie: Fişierul se poate executa direct, în editorul Scilab-ului folosind din meniul ‘Execute’ opţiunea ‘Load into Scilab’.
Calculul limitei şirurilor definite recurent: să se calculeze limita şirului lim c n , unde c1 = 10 si c n+1 =
n→∞
2 + c n2 ,∀n ≥ 1 . 2c n
Se scrie programul în fişierul numit SirRecurent.sce:
x=10; c=(2+x^2)/(2*x); while abs(c-x)>0.00001 x=c; c=(2+x^2)/(2*x);
end c Se lansează în execuţie în fereastra de comandă scriind:
-->exec("SirRecurent.sce") şi se obtine rezultatul:
c = 1.4142136 Exemple
1. Să se genereze o matrice A, cu n linii şi n+1 coloane, ale cărei elemente sunt:
2, daca i = j . A = − 1, daca i − j = 1 0, in rest
Soluţie: n=10; for i=1:n for j=1:n+1 if i==j A(i,j)=2; elseif abs(i-j)==1 A(i,j)=-1; else A(i,j)=0; end end end A In fereastra de comandă o sa apară matricea definită mai sus
A = ! 2. - 1. ! - 1.
0.
0.
0.
0.
0.
0.
0.
0.
0. !
2. - 1.
0.
0.
0.
0.
0.
0.
0.
0. !
! 0. - 1.
2. - 1.
0.
0.
0.
0.
0.
0.
0. !
2. - 1.
0.
0.
0.
0.
0.
0. !
2. - 1.
0.
0.
0.
0.
0. !
2. - 1.
0.
0.
0.
0. !
2. - 1.
0.
0.
0. !
2. - 1.
0.
0. !
2. - 1.
0. !
! 0.
0. - 1.
! 0.
0.
0. - 1.
! 0.
0.
0.
0. - 1.
! 0.
0.
0.
0.
0. - 1.
! 0.
0.
0.
0.
0.
0. - 1.
! 0.
0.
0.
0.
0.
0.
0. - 1.
! 0.
0.
0.
0.
0.
0.
0.
0. - 1.
2. - 1. !
2. Să se scrie un program utilizând while care calculează suma elementelor vectorului x=(5 2 –9 10 –1 9 –1) până când întâlneşte un număr mai mare ca 8.
Soluţie: x=[5 2 -9 10 -1 9 -1]; s=0; k=1; while (x(k)<=8)&(k<=length(x)) s=s+x(k); k=k+1; end s In fereastra de comandă o să apară rezultatul:
--> --> s = - 2. 3. Să se calculeze limita următorului şir definit recurent: lim x n , unde x 0 = 1 si x n +1 =
n→∞
Soluţie: x=1; eps=0.00001; c=x/sqrt(x^2+1); while abs(c-x)>eps
xn x n2
+1
,∀n ≥ 0
x=c; c=x/sqrt(x^2+1); end c In fereastra de comandă o să apară rezultatul:
c =
0.0271363 10
10
k =1
k =1
4. Să se calculeze şi să se afişeze A = ∑ (2k − 1)3 + ∏ 3k 2
Solutie: s=0;p=1; for i=1:10 s=s+(2*i-1)^3; p=p*3*i^2; end a=s+p In fereastra de comandă o să apară rezultatul:
a =
7.776D+17
Exerciţii
1. Să se calculeze suma elementelor unui vector care sunt cuprinse între valorile a şi b cu a
Soluţie: v=[2 5 7 -2 5 8 0 12 32]; s=0; a=-1;b=7; for i=1:length(v) if (v(i)>a)&(v(i)
Soluţia va aparea în fereastra de comandă:
--> s =
12.
2. Să se calculeze produsul elementelor diferite de zero ale unui vector. 3. Se dă şirul de numere x1, x2, …, xn. Să se determine numărul de elemente pozitive şi să se calculeze produsul lor. 4. Să se calculeze media aritmetică a elementelor unui vector care sunt mai mari decat o valoare dată. 5. Se dă şirul de numere x1, x2, …, xn. Să se calculeze: - media aritmetică a numerelor pozitive; - suma pătratelor numerelor negative 6. Să se ordoneze crescător elementele unui vector. 7. Să se calculeze produsul scalar al vectorilor x = ( x1, x2, …, xn) şi y = (y1, y2, …,yn). 8. Se dau două siruri de câte n numere: a1, a2, …, an şi b1, b2, …, bn. Să se formeze şirul c1, c2, …, cn în care termenii se obţin după următorul procedeu: pentru orice i, i = 1, n , ai + bi , daca ai ⋅ bi < 0 ci = 2 max{ai , bi }, daca ai ⋅ bi ≥ 0
II. 3.2. Fişiere de tip funcţie
Fişierele funcţie sunt utilizate pentru crearea unor noi funcţii SCILAB. Forma generală a unui fişier funcţie este:
(parametri_intrare) trare) function [parametri_ieşire] = nume_funcţie (parametri_in …..instrucţiuni de definire a funcţiei… endfunction unde:
function: function cuvânt cheie care declară fişierul ca fişier funcţie (obligatoriu);
nume_functie: nume_functie numele funcţiei, adică numele sub care se salvează fişierul, fără extensie; nu poate fi identic cu cel al unui fişier existent;
parametri_iesire: parametrii de ieşire trebuie separaţi cu virgulă şi cuprinsi între paranteze drepte. Dacă funcţia nu are parametri de ieşire parantezele drepte şi semnul egal (=) nu mai au sens;
parametri_intrare: parametrii de intrare trebuie separaţi cu virgulă şi cuprinsi între paranteze rotunde. Daca funcţia nu are parametri de intrare parantezele rotunde nu mai au sens.
endfunction: instrucţiunea de terminare a unei funcţii.
Observaţii: Observaţii Un fişier de tip funcţie poate să contină mai multe funcţii. Pentru utilizarea unei funcţii trebuie folosită, mai intâi, comanda:
getf(“NumeFunctie.sci”) Spre deosebire de fişierele script care au extensia ‘sce’, fişierele de tip funcţie au extensia ‘sci’. Exemplul 1:
Calculul factorialului: se vor defini în acelaşi fişier (functii.sci) două funcţii care calculează factorialul unui număr dat; prima nu verifică dacă argumentul este întreg pozitiv, a doua funcţie face acest lucru înainte de a face calculul
function [f]=factorial(n) f=prod(1:n) endfunction function [f]=factorial1(n) //dacă n nu este întreg şi pozitiv se transmite un mesaj de eroare, //apoi se transformă numărul într-un număr întreg pozitiv if (n-floor(n)~=0)|n<0 then n=floor(abs(n)) warning('argumentul nu este întreg pozitiv; se va calcula '+sprintf("%d",n)+"!") end f=prod(1:n) endfunction În fereastra de comandă se vor folosi aceste funcţii astfel:
-->getf('functii.sci') -->factorial(5) ans =
120.
-->factorial(0.5) ans =
1.
-->factorial(-0.5) ans =
1.
Dacă argumentul nu este întreg şi pozitiv funcţia factorial întoarce un rezultat oarecare, 1. Dacă vom folosi funcţia factorial1 se va obţine:
-->factorial1(5) ans =
120.
-->factorial1(0.5)
WARNING:argumentul nu este întreg pozitiv; se va calcula 0! ans =
1.
-->factorial1(-5) WARNING:argumentul nu este întreg pozitiv; se va calcula 5! ans =
120.
Recursivitate: o funcţie se poate apela pe ea însăşi. Exemplul 2: 2 calculul factorialului recursiv
function [f]=fact(n) if n<=1 then f=1 else f=n*fact(n-1) end endfunction Exemplul 3: 3 exemplul următor calculează polinoamele Chebyshev care verifică relaţia de recurenţă: Tn +1 ( x ) = 2 xTn ( x ) − Tn −1 ( x ) .
function ch=chebyshev(x,n) if n==0 ch=1; elseif n==1 ch=x; else ch=2*x*chebyshev(x,n-1)-chebyshev(x,n-2); end endfunction Observaţie: Observaţie O funcţie se poate defini direct prin intermediul comenzii deff
deff(‘[y1, y2,…, ]=nume_functie(x1, x2, …)’,text)
Exemplu: deff(‘[y]=f(x)’,’y=sin(x).*cos(x)’)
Exerciţii
1. Să se calculeze şi să se afişeze suma pătratelor elementelor mai mari decat 1 ale unei matrice pătratice.
Soluţie: function [s]=sumaP(a) [m,n]=size(a) s=0; for i=1:m for j=1:n if a(i,j)>=1 s=s+a(i,j)^2 end end end endfunction 2. Să se calculeze media aritmetică a elementelor de pe diagonala principală a unei matrice pătratice de ordinul n, cu n<10
Soluţie: function [ma]=media(a) n=size(a,'r') ma=0; for i=1:n ma=ma+a(i,i) end ma=ma/n endfunction
3. Să se determine valorile maxime şi minime ale unui vector şi ale unei matrice, apoi să se determine şi indicele (poziţia acestora).
Soluţie: -->a=ceil(10*rand(3,3)) a = ! 3.
4.
9. !
! 8.
7.
7. !
! 1.
7.
9. !
-->[m,n]=max(a) n = ! 1.
3. !
m = 9. -->[mi,p]=min(a) p = ! 3. mi =
1. ! 1.
4. Să se determine valorile şi vectorii proprii ale matricei 0 − 2 1 1 − 2 1 . 0 1 − 2
Soluţie: In determinarea soluţiilor nebanale ale ecuaţiei Ax = λx unde:
A – este o matrice patratică de ordinul n, x – este vector coloană de ordinul n, iar λ - este un scalar,
valorile x, respectiv λ , care satisfac ecuaţia de mai sus se numesc vectori proprii respectiv valori proprii. -->a=[-2 1 0;1 -2 1;0 1 -2]; -->[v,d]=spec(a) d =! - 3.4142136
0.
0.
!
!
0.
!
0.
- 2.
0.
!
0. - 0.5857864 !
v =! 0.5
- 0.7071068
- 0.5
! - 0.7071068
2.082D-17 - 0.7071068 !
! 0.5
0.7071068
- 0.5
! !
5. Să se calculeze suma elementelor diferite de o valoare dată, v, ale unei matrice.
Soluţie: function [s]=sdv(a,v) [m,n]=size(a) s=0; for i=1:m for j=1:n if a(i,j)~=v s=s+a(i,j) end end end endfunction 6. Să se determine elementul maxim al unui matrice de numere reale. 7. Să se calculeze produsul elementelor diferite de zero de pe diagonala principală (secundară) ale unei matrice.
Soluţie: function [p]=prodDS(a) n=size(a,'c') p=1; for i=1:n if a(i,n-i+1)~=0 p=p*a(i,n-i+1) end end
endfunction 8. Să se calculeze suma elementelor unei matrice situate deasupra (dedesubtul) diagonalei principale (secundare).
Soluţie: function [s]=sdeasupra(a) n=size(a,'r') s=0; for i=1:n-1 for j=i+1:n s=s+a(i,j) end end endfunction 9. Să se realizeze un program care calculează produsul a două matrice.
Solutie: function [c]=f(a,b) [m,n]=size(a) [n,p]=size(b) for i=1:m for j=1:p c(i,j)=0 for k=1:n c(i,j)=c(i,j)+a(i,k)*b(k,j) end end end endfunction 10.
Să se determine numărul de coloane cu toate elementele
negative dintr-o matrice de tip m × n .
Solutie: function [cn]=colneg(a)
[m,n]=size(a); cn=0; j=1; while j0 j=j+1; i=1; else i=i+1; if i>m cn=cn+1; j=j+1; end end end end endfunction 11.
Să se determine o linie a unei matrice de tip m × n care
conţine elementul a cărui valoare absolută este mai mare decât toate celelalte elemente ale matricei.
Solutie: function [L]=linie(a) [m,n]=size(a) v=a(1,1) for i=1:m for j=1:n if abs(a(i,j))>v L=i end end end
endfunction 12.
Se dau numerele aij , i = 1,2,...,100; j = 1,2,...,50 si b j j = 1,2,...,50 . Să
se realizeze un program care calculează:
50
ci = ∑ aij b j , i = 1,2,...,100 . j =1
Soluţie: function [c]=f(a,b) [m,n]=size(a) for i=1:m c(i)=0 for j=1:n c(i)=c(i)+a(i,j)*b(j) end end endfunction
III. 1. Rezolvarea ecuaţiilor
III. III. 1.1. Rezolvarea ecuaţiilor algebrice
Forma canonică a unei ecuaţii algebrice este: an ⋅ x n + an −1 ⋅ x n −1 + ... + a1 ⋅ x + a0 = 0
unde n este gradul polinomului (pentru a rezolva ecuaţia în Scilab n ≤ 100) şi ak sunt coeficienţii polinomului care pot fi numere reale
sau complexe. Pentru a rezolva o ecuaţie algebrică mai întâi, se defineşte expresia polinomului cu ajutorul funcţiei poly apoi, cu ajutorul funcţiei roots roots se determină rădăcinile ecuaţiei. Exemplu: Exemplu Să se rezolve ecuaţia x 4 − 12 x 3 + 22 x 2 − 20 x = 0 . -->p=poly([0,-20,22,-12,1],'x','coef') p = 2
3
4
- 20x + 22x - 12x + x -->sol=roots(p) sol = ! 0
!
! 1. + i
!
! 1. - i
!
! 10.
!
III. 1.2. Rezolvarea ecuaţiilor transcendente Orice ecuaţie se poate scrie sub forma f(x) = g(x) echivalentă cu f(x) – g(x) = 0. După ce se reprezintă grafic funcţia f(x) – g(x) pe
intervalul de interes putem să evaluăm valorile aproximative în care funcţia se anulează. Se alege drept valoare iniţială pentru x valoarea pentru care, conform graficului, funcţia se anulează. Pentru a determina soluţia ecuaţiei se apelează funcţia fsolve a Scilab-ului. Exemplu: Exemplu Să se resolve ecuaţia sin( x ) = 1 − 2 x --> x=0:0.1:2; --> y1=sin(x); --> y2=1-2.*x; --> plot(x,y1-y2) Se va obţine următorul grafic:
1.0
0.5
0.0 0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
-0.5
-1.0
-1.5
-2.0
-2.5
-3.0
Se observă că soluţia ecuaţiei se află aproape de x = 0.3; pentru aflarea soluţiei cu o eroare dată se va folosi următoarea secvenţă: function y=f(x) y=sin(x)+2.*x-1 endfunction [s]=fsolve(0.3,f,10^(-10))
obţinându-se pentru soluţie valoarea: s =
0.3354180
Pentru verificare se calculează f(s) care se obţine în acest caz 0. Exerciţii: Exerciţii Sa se rezolve urmatoarele ecuatii: 1. x4-2x3-x+2=0; 2. x3+3x+1=0; 3. ( x − 1)2 − 2 ⋅ e x ; 4.
x x − 2 ⋅ cos π ⋅ . 2
III. 2. Rezolvarea sistemelor de ecuaţii liniare liniare Să se rezolve următoarele sisteme liniare de ecuaţii: a)
9 x1 + x2 − x3 = 18 x1 + 9 x2 − x4 = 19 − x1 + 9 x3 + x4 = 8 − x2 + x3 + 9 x4 = 9
b)
x1 + 2 x 2 + 3 x 3 + 4 x 4 2 x 1 + 3 x 2 + 4 x 3 + x 4 3 x 1 + 4 x 2 + 2 x 3 + x 4 4 x 1 + x 2 + 2 x 3 + 3 x 4
= 11 = 12 = 13 = 14
Solutie; rezolvarea sistemului de la punctual a) in Scilab -->a=[9 1 -1 0;1 9 0 -1;-1 0 9 1;0 -1 1 9]; b=[18;19;8;9]; -->x=a^(-1)*b x =! 1.8831169 ! ! 2.025974 ! ! 0.9740260 ! ! 1.1168831 ! 1.
Să se calculeze factorul Cholesky al matricei: 1 −1 0 9 9 0 − 1 ; 1 A= −1 0 9 1 0 −1 1 9
Factorizarea Cholesky constă în descompunera matricei într-un produs de forma A=RtR; (R=chol(A)) apoi să se rezolve sistemul de la ex.1 punctul a.
Solutie: -->a=[9 1 -1 0;1 9 0 -1;-1 0 9 1;0 -1 1 9]; b=[18;19;8;9];
-->R=chol(a) R =! 3.
0.3333333 - 0.3333333
! 0.
2.981424
! 0.
0.
2.981191
! 0.
0.
0.
0.
!
0.0372678 - 0.3354102 ! 0.3396294 ! 2.9617819 !
-->y=(R')^(-1)*b y =! 6.
!
! 5.7019733 ! ! 3.2830838 ! ! 3.3079642 ! -->x=R^(-1)*y x =! 1.8831169 ! ! 2.025974 ! ! 0.9740260 ! ! 1.1168831 ! 2.
Să se factorizeze LU (lower – upper) matricea: 1 −1 0 9 9 0 − 1 ; 1 A= −1 0 9 1 0 −1 1 9
Factorizarea LU: o matrice pătrată este descompusă sub forma produsului a două matrice triunghiulare, una inferior triunghiulară, cu elemente 1 pe diagonala principală şi cealaltă superior triunghiulară ([L,U]=lu(A)) apoi să se rezolve sistemul de la ex.1 punctul a.
Solutie: -->a=[9 1 -1 0;1 9 0 -1;-1 0 9 1;0 -1 1 9]; b=[18;19;8;9]; -->[L U]=lu(a) U =! 9.
1.
- 1.
0.
! 0.
8.8888889
! 0.
0.
8.8875
! 0.
0.
0.
L =! 1.
0.
!
0.1111111 - 1. 1.0125
! !
8.7721519 ! 0.
0. !
! 0.1111111
1.
! - 0.1111111
0.0125
! 0.
- 0.1125
0.
0. ! 1.
0. !
0.1139241
1. !
-->y=L^(-1)*b y =! 18.
!
! 17.
!
! 9.7875
!
! 9.7974684 ! -->x=U^(-1)*y x =! 1.8831169 ! ! 2.025974 ! ! 0.9740260 ! ! 1.1168831 ! III. 3. 3. Rezolvarea sistemelor de ecuaţii neliniar neliniare liniare x 2 + axy + y 2 + x + sin( x ) − 4a = 0 Fie sistemul 2 x 2 + ay 3 − 3,5 = 0
Să se rezolve sistemul pentru a=1 folosind ca punct de start (1,1).
Rezolvarea cu funcţia fsolve a SCILAB-ului: function y=f(x) y=zeros(2,1); y(1)=x(1)^2+x(1)*x(2)+x(2)^2+x(1)+sin(x(1))-4; y(2)=2*x(1)^2+x(2)^3-3.5; endfunction [s]=fsolve([1;1],f,0.0001) Se obtine rezultatul: s =! 0.4905521 ! ! 1.4452428 !
Utilizarea Calculatoarelor: Laboratoarele 6 si 7 Reprezentări grafice în Scilab
Tabel: Reprezentări grafice în SCILAB: SCILAB: Funcţia plot
Semnificaţia Reprezintă grafic o curba în 2D unind punctele de coordonate (x,y) plot2d Reprezintă grafic o funcţie în 2D fplot2d Reprezintă grafic o funcţie în 2D definite anterior paramfplot2d Reprezintă grafic cu animaţie o funcţie în 2D plot3d Reprezintă grafic o suprafaţă în 3D contour2d Reprezintă grafic liniile de contur în2D contour Reprezintă grafic liniile de contur ale unei suprafeţe Observaţie: Observaţie Instrucţiunea plot va desena într-o altă fereastră graficul; pentru mai multe informaţii şi exemple tastaţi în fereastra de comandă: apropos Graphics. Exemplu plot x=linspace(0,2*%pi,101); y=exp(x).*sin(4*x); plot(x,y,'x','y','y=exp(x)*sin(4x)') y
y=exp(x)*sin(4x)
200
100
0
-100
-200
-300
x
-400 0
1
2
3
4
5
6
7
Exemplu plot2d
t=(0:0.1:6*%pi); plot2d(t',sin(t)'); xtitle('plot2d and xgrid ','t','sin(t)'); xgrid(); sin(t)
plot2d and xgrid
1.0 0.8
0.6 0.4
0.2
0.0 -0.2
-0.4 -0.6
-0.8 t
-1.0 0
2
4
6
8
10
12
14
16
18
20
Variante ale instrucţiunii plot2d: plot2d2, plot2d3 Aceste instrucţiuni se utilizează exact ca plot2d: se foloseşte aceeaşi sintaxă şi aceleaşi argumente opţionale. • plot2d2: permite desenarea unei funcţii în scară; în loc să traseze un segment de dreaptă între punctele ( xi , yi )
şi
( xi +1 , yi +1 ) , plot2d2 trasează un segment orinzontal între ( xi , yi ) şi ( xi +1 , yi ) , apoi un segment vertical între ( xi +1 , yi ) şi ( xi +1 , yi +1 ) . Exemplu: Exemplu n=10; x=(0:n)'; y=x; xbasc() plot2d2(x,y,style=2,frameflag=5,rect=[0,-1,n+1,n+1]) xtitle("Exemplu de grafic cu plot2d2")
Exemplu de grafic cu plot2d2 11
9
7
5
3
1
-1 0
2
4
6
8
10
12
• plot2d3 desenează diagrame cu bare verticale: pentru fiecare punct ( xi , yi ) plot2d3 plot2d3 trasează un segment vertical între ( xi ,0) şi ( xi , yi ) . Exemplu: Exemplu n=6; x=(0:n)'; y=binomial(0.5,n)'; xbasc() plot2d3(x,y,style=2,frameflag=5,rect=[-1,0,n+1,0.32]) xtitle("Exemplu de grafic cu plot2d3") Exemplu de grafic cu plot2d3 0.32
0.28
0.24
0.20
0.16
0.12
0.08
0.04
0.00 -1
0
1
2
3
4
5
6
7
Exemplu fplot2d deff("[y]=f(x)","y=sin(x)+cos(x)") x=[0:0.1:10]*%pi/10; fplot2d(x,f) clf(); fplot2d(1:10,'parab') 100 90
80 70
60
50 40
30 20
10 0 1
2
3
4
5
6
7
8
9
10
Exemplu paramfplot2d deff('y=f(x,t)','y=t*sin(x)') x=linspace(0,2*%pi,50);theta=0:0.05:1; paramfplot2d(f,x,theta); 1.0 0.8
0.6 0.4
0.2
0.0 -0.2
-0.4 -0.6
-0.8 -1.0 0
1
2
3
4
5
6
7
Trasarea pe acelaşi grafic a mai multor multor curbe care nu au acelaşi număr de puncte x1=linspace(0,1,61)'; x2=linspace(0,1,31)'; x3=linspace(0.1,0.9,12)'; y1=x1.*(1-x1).*cos(2*%pi*x1); y2=x2.*(1-x2); y3=x3.*(1-x3)+0.1*(rand(x3)-0.5);//la fel cu y2 dar cu o mica perturbatie ymin=min([y1;y2;y3]);ymax=max([y1;y2;y3]); dy=(ymax-ymin)*0.05; dreptunghi=[0,ymin-dy,1,ymax+dy]; xbasc() plot2d(x1,y1,style=1,frameflag=5,rect=dreptunghi) plot2d(x2,y2,style=2,frameflag=0) plot2d(x3,y3,style=-1,frameflag=0) xtitle("Mai multe curbe pe acelasi grafic")
Mai multe curbe pe acelasi grafic 0.4
0.3 +
+ +
+ +
0.2
+
+
+
+
+ +
0.1 +
0.0
-0.1
-0.2
-0.3 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Exemplu Exemplu plot3d t=[0:0.3:2*%pi]'; z=sin(t)*cos(t'); plot3d(t,t,z)
1.0 Z
0 -1.0 0
0 1
1 2
2 3
3 4
4
Y
5
X
5 6
6
Exemplu plot3d1: colorează numai după valorile lui z x=linspace(0,2*%pi,31); z=cos(x)'*cos(x); C=hotcolormap(32); xset("colormap",C) xset("hidden3d",30) xbasc() plot3d1(x,x,z,flag=[1 4 4])
1.0 Z
0 -1.0 0
0 1
1 2
2 3
3 4
Y
4 5
5 6
6 7
7
X
Exemplu: Exemplu Reprezentarea color, în plan a liniilor de contur ale unei suprafeţe folosind Sgrayplot: x=-10:10; y=-10:10;m =rand(21,21); Sgrayplot(x,y,m,"111",[-20,-20,20,20]) t=-%pi:0.1:%pi; m=sin(t)'*cos(t);clf() Sgrayplot(t,t,m) 3
2
1
0
-1
-2
-3 -3
-2
-1
0
1
2
3
+
Exemplu: Exemplu Reprezentarea liniilor de contur în plan a unei suprafeţe contour2d(1:10,1:10,rand(10,10),5,rect=[0,0,11,11]) xset("fpf","%.2f") clf() contour2d(1:10,1:10,rand(10,10),5,rect=[0,0,11,11])
11.0 0.17
0.66 0.66
8.8
0.50 0.33 0.83 0.66 0.50 0.33 0.17
0.33 0.66 0.83
6.6
0.50
0.66
0.33 0.33 0.17
0.17 0.33 0.50 0.66 0.83 0.17
0.66
0.50
0.17 0.83
0.170.33 0.66 4.4
0.33 0.66 0.500.66
0.33
0.33
0.17
0.33
0.50 0.66
2.2
0.33 0.50 0.66
0.33 0.17
0.83 0.66
0.50 0.66 0.83 0.50 0.33
0.50 0.33 0.17 0.66 0.83
0.0 0.0
2.2
4.4
6.6
8.8
11.0
Exemplu contour t=%pi*[-10:10]/10; deff("[z]=surf(x,y)","z=sin(x)*cos(y)"); z=feval(t,t,surf); rect=[-%pi,%pi,-%pi,%pi,-1,1]; contour(t,t,z,10,35,45," ",[0,1,0],rect) xset("fpf","%.2f") ;xbasc() contour(t,t,z,10,flag=[0,1,0],ebox=rect)
Reprezentări grafice în plan 1. Să se definească urmatoarele funcţii cu acoladă: − x , − 4 ≤ x < 0 f ( x) = x, 0 ≤ x ≤ 4
si
sin( x ), {− 10 ≤ x < −4}∪ {4 ≤ x ≤ 10} f ( x) = cos(5 x ), − 4 ≤ x ≤ 4
şi să se reprezinte grafic.
Soluţie function [y]=f(x) y=abs(x) endfunction x=linspace(-10,10,100); plot2d(x,f(x)) xtitle('Reprezentarea grafica a functiei modul')
Reprezentarea grafica a functiei modul 10 9
8 7
6
5 4
3 2
1 0 -10
-8
-6
-4
-2
0
2
4
6
8
10
2. Să se reprezinte grafic conturul liniilor de nivel ale funcţiei 2 2 z = x ⋅ e − x − y in domeniul − 2 ≤ x ≤ 2, − 2 ≤ y ≤ 3 .
Soluţie x=linspace(-2,2,10); y=linspace(-2,3,10); for i=1:10 for j=1:10 z(i,j)=x(i)*exp(-x(i)^2-y(j)^2); end end contour2d(x,y,z,10) 1.9
1.5 -0.037 0.037 1.1 -0.259
0.7
-0.111
0.185
-0.333
0.111
0.259 -0.185
0.3
-0.1
-0.5
-0.9
-1.3
-1.7 -2.0
-1.6
-1.2
-0.8
-0.4
0.0
0.4
0.8
1.2
1.6
2.0
3. Folosind ecuaţiile parametrice să se reprezinte grafic în plan urmatoarele curbe:
x = a + r cos( t ) , t ∈ [0,2π] y = b + r sin( t )
a) cercul: ( x − a )2 + ( y − b)2 = r 2 ;
Soluţie coordonate polare: a=1;b=2;r=2; t=linspace(0,2*%pi,100); x=r*cos(t)+a; y=r*sin(t)+b; plot2d(x,y,frameflag=4,axesflag=5) Reprezentarea grafică a cercului de centru (1, 2) şi rază 2: coordinate polare 4.0 3.6
3.2 2.8
2.4
2.0 1.6
1.2 0.8
0.4 0.0 -1.60
b) elipsa:
x2 a
2
+
-0.73
0.13
1.00
1.87
y2
x = a cos( t ) = 1; , t ∈ [0,2π] b y = b sin( t ) 2
Soluţie a=2;b=1; t=linspace(0,2*%pi,100); x=a*cos(t); y=b*sin(t); plot2d(x,y,frameflag=4,axesflag=5)
2.73
3.60
Reprezentarea grafică a elipsei 1.500
1.071
0.643
0.214
-2.0
-1.6
-1.2
-0.8
-0.4
0.0
0.4
0.8
1.2
1.6
2.0
-0.214
-0.643
-1.071
-1.500
c) hiperbola:
x2 a
2
−
y2
x = ach( t ) π π = 1; , t ∈ [− , ] (ramura din dreapta) 4 4 b y = bsh( t ) 2
Soluţie a=2;b=1; t=linspace(-%pi/4,%pi/4,100); x=a*cosh(t); y=b*sinh(t); plot2d(x,y,frameflag=4,axesflag=5) Reprezentarea grafică a hiperbolei, ramura din dreapta
0.9
0.7
0.5
0.3
0.1
-0.1
1.200
1.660
2.120
2.580
3.040
3.500
-0.3
-0.5
-0.7
-0.9
d) astroïda:
2 2 2 x = a cos 3 ( t ) 3 3 3 , t ∈ [ 0, 2 π ] x +y =a ; 3
y = a sin ( t )
Soluţie a=2; t=linspace(0,2*%pi,100); x=a*(cos(t))^3; y=a*(sin(t))^3; plot2d(x,y,frameflag=4,axesflag=5) 2.0 1.6
1.2 0.8
0.4
-2.60
-1.73
0.0 0.00
-0.87
0.87
1.73
-0.4
-0.8 -1.2
-1.6 -2.0
x = a ( t − sin( t )) , t ∈ [0,4π] y = a (1 − cos( t )
e) cicloïda:
Soluţie a=2; t=linspace(0,4*%pi,100); x=a*(t-sin(t)); y=a*(1-cos(t)); plot2d(x,y,frameflag=4,axesflag=5) Reprezentarea grafică a cicloïdei 12.50
9.00
5.50
2.00
0 -1.50
-5.00
-8.50
4
8
12
16
20
24
28
2.60
x = a ( 2 cos( t ) − cos( 2t )) , t ∈ [0,2π] y = a ( 2 sin( t ) − sin( 2t ))
f) cardioïda:
Soluţie a=2; t=linspace(0,2*%pi,100); x=a*(2*cos(t)-cos(2*t)); y=a*(2*sin(t)-sin(2*t)); plot2d(x,y,frameflag=4,axesflag=5) Reprezentarea grafică a cardioïdei
6
4
2
-9.00
-6.86
-4.71
-2.57
0 -0.43
1.71
3.86
6.00
-2
-4
-6
x = a cos(2t ) cos(t )
g) lemniscata lui Bernoulli:
y = a cos(2t ) sin(t )
, t ∈ [0, 2π ]
Soluţie a=2; t=linspace(0,2*%pi,100); x=a*sqrt(cos(2*t)).*cos(t); y=a*sqrt(cos(2*t)).*sin(t); plot2d(x,y,frameflag=4,axesflag=5) Reprezentarea grafică a lemniscatei lui Bernoulli
1.5
1.0
0.5
0.0 -2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-0.5
-1.0
-1.5
Reprezentări grafice în spaţiu 1. Să se reprezinte grafic în 3D liniile de contur ale funcţiilor a. z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2)
Soluţie t=linspace(-3,3,100); deff("[z]=surf(x,y)","z =3*(1-x).^2.*exp(-(x.^2)-(y+1).^2)"); z=feval(t,t,surf); rect=[-3,3,-3,3,-2,2]; xset("fpf","%.2f") xbasc() contour(t,t,z,10,flag=[0,1,0],ebox=rect)
b. z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) - 1/3*exp(-(x+1).^2 - y.^2) 2. Să se reprezinte grafic în 3D funcţia z = x 2 + y 2 , x , y ∈ [−3,3] .
Soluţie x=linspace(-3,3,100); y=linspace(-3,3,100); for i=1:length(x) for j=1:length(y) z(i,j)=x(i)^2+y(j)^2; end end plot3d(x,y,z)
3. Să se reprezinte grafic în 3D urmatoarele cuadrice: a. Sfera: x 2 + y 2 + z 2 = 4
Soluţie u=linspace(0,2*%pi,100); v=linspace(0,%pi,100); r=4; for i=1:length(u) for j=1:length(v)
x(i,j)=r*sin(v(j))*cos(u(i)); y(i,j)=r*sin(u(i))*sin(v(j)); z(i,j)=r*cos(v(j)); end end plot3d(x,y,z)
b. Elipsoidul:
x2 y2 z2 =1 + + 4 9 16
Soluţie u=linspace(0,2*%pi,100); v=linspace(0,%pi,100); a=2;b=3;c=4; for i=1:length(u) for j=1:length(v) x(i,j)=a*sin(v(j))*cos(u(i)); y(i,j)=b*sin(u(i))*sin(v(j)); z(i,j)=c*cos(v(j)); end
end plot3d(x,y,z)
x 2 y2 z2 c. hiperboloidul cu o panză: =1 + − 4 9 16
d. hiperboloidul cu două panze: e. paraboloidul eliptic: x 2 +
x2 y2 z2 +1 = 0 + − 4 9 16
y2 = 2z 4
Soluţie h=10; u=linspace(0,2*%pi,100); v=linspace(0,h,100); a=2; for i=1:length(u) for j=1:length(v) x(i,j)=v(j)*cos(u(i)); y(i,j)=v(j)*sin(u(i)); z(i,j)=1/2*(x(i,j)^2+y(i,j)^2/4);
end end plot3d(x,y,z)
f. paraboloidul hyperbolic:
x 2 y2 = 2z − 4 9
Soluţie h=10; u=linspace(0,2*%pi,100); v=linspace(0,h,100); a=2; for i=1:length(u) for j=1:length(v) x(i,j)=v(j)*cos(u(i)); y(i,j)=v(j)*sin(u(i)); z(i,j)=1/2*(x(i,j)^2/4-y(i,j)^2/9); end end plot3d(x,y,z)
4. Să se reprezinte grafic cilindrul: x 2 + y 2 = a 2 , z ∈ [0, h]
Soluţie h=10; u=linspace(0,2*%pi,100); v=linspace(0,h,100); a=2; for i=1:length(u) for j=1:length(v) x(i,j)=a*cos(u(i)); y(i,j)=a*sin(u(i)); z(i,j)=v(j); end end plot3d(x,y,z)
5. Să se reprezinte grafic conul: x 2 + y 2 = z 2 , z ∈ [0, h]
Solutie h=10; u=linspace(0,2*%pi,100); v=linspace(0,h,100); a=2; for i=1:length(u) for j=1:length(v) x(i,j)=v(j)*cos(u(i)); y(i,j)=v(j)*sin(u(i)); z(i,j)=sqrt(x(i,j)^2+y(i,j)^2); end end plot3d(x,y,z)
6. Să se reprezinte grafic paraboloidul: x 2 + y 2 = z , z ∈ [0, h]
Soluţie h=10; u=linspace(0,2*%pi,100); v=linspace(0,h,100); a=2; for i=1:length(u) for j=1:length(v) x(i,j)=v(j)*cos(u(i)); y(i,j)=v(j)*sin(u(i)); z(i,j)=x(i,j)^2+y(i,j)^2; end end plot3d(x,y,z)
Utilizarea Calculatoarelor: Laboratorul Laboratorul al 88-lea I. Integrare numerică
I.1. Formule de cuadratură
Fie f : [a , b] → ℜ . Pentru calculul aproximativ al integralei b
I = ∫ f ( x )dx se pot folosi următoarele formule: a
Formula trapezelor : dacă h =
b−a şi n
a=x0,x1,…,xn=b este o diviziune a intervalului [a,b]
atunci: I≈
n −1 h n h ( f ( x ) + f ( x ) ) = f ( a ) + f ( b ) + 2 f ( x i ) ∑ ∑ i −1 i 2 i =1 2 i =1
( b − a )3 max| f ( 2 ) Eroarea e majorată de er ≤ . 12n 2
Formula lui Simpson : dacă h =
b−a şi 2n
a=x0, x1,…, x2n=b este o diviziune a intervalului
[a,b] atunci I≈
n n −1 h n ( f ( x2i − 2 ) + 4 f ( x2i −1 ) + f ( x2i )) = h f (a ) + f (b) + 4 ∑ f ( x2i −1 ) + 2 ∑ f ( x2i ) ∑ 3 i =1 3 i =1 i =1
Eroarea este majorată de: er ≤
(b − a)5 2880
max f ( 4 )
1 . n4
Exemplu: Exemplu
Fie
functia
f(x)=(2x+5)/(x2+1). Să
se
calculeze
5
integrala I = ∫ f ( x )dx folosind cele două formule pentru n=2, n=10, 2
n=100 şi să se compare rezultatele. Rezolvare cu funcţia Scilab-ului În Scilab există functia integrate – integrare prin cuadraturi, care are formula
[x]=integrate(expr,v,x0,x1 [,ea [,er]]) unde expr: expresie Scilab (de obicei expresia funcţiei ca şir de caractere!!),
v: şir de caractere, variabila de integrare, x0,x1: numere reale, limitele de integrare, ea,er: numere reale, valoare erorii absolute. Valoarea implicită e 0; er: număr real, valoarea erorii relative. Valoarea implicită 1.d-8 x1
Această funcţie calculează I =
∫ f (v )dv cu o acurateţe ce satisface x0
relaţia:
| I − x |<= max(ea, er* | I |) Exemple de folosire: integrate('sin(x)','x',0,%pi)
integrate(['if x==0 then 1,';'else sin(x)/x,end'],'x',0,%pi)
Rezolvare exemplu
1.
Integrare cu funcţia Scilab
-->integrate('(2*x+5)/(x^2+1)','x',2,5) ans = 2.9799189
2.
Aproximarea integralei cu formula trapezelor
function f=f(x) f=(2*x+5)/(x^2+1); endfunction function I=IntegrareTrapeze(a,b,n) h=(b-a)/n; for i=1:n+1 x(i)=a+(i-1)*h; end s=0; for i=2:n s=s+f(x(i)); end I=h/2*(f(a)+f(b)+2*s); endfunction IntegrareTrapeze(2,5,10) ans = 2.9866206 3.
Aproximarea integralei cu formula lui Simpson
Pentru f definit deasupra avem: function I=IntegrareSimpson(a,b,n)
h=(b-a)/(2*n); for i=1:2*n+1 x(i)=a+(i-1)*h; end s=0; for i=1:n s=s+f(x(2*i)); end s1=0;
for i=2:n s1=s1+f(x(2*i-1)); end I=h/3*(f(a)+f(b)+4*s+2*s1); endfunction IntegrareSimpson(2,5,5) ans =
2.9799626
Exerciţii: Să se aproximeze integralele: 1
∫
1.
0 5
dx
;
1+ x2 dx
∫ x ln 2 x ;
2.
2
3
∫
3.
1 x x
dx 2
+ 5x + 1
I. 2. Formule de Cubatura Cubatura bd
Pentru calculul aproximativ al integralei I = ∫ ∫ f ( x , y )dxdy se pot ac
folosi următoarele formule: Formula trapezelor b−a d −c , r= , x i = a + i ⋅ h, i = 0, n şi y j = c + j ⋅ r , j = 0, m atunci: n m h⋅ r n m I≈ ∑ ∑ f xi , y j + f x i −1 , y j + f x i , y j −1 + f x i −1 , y j −1 4 i =1 j =1
dacă h =
( (
) (
) (
) (
))
Formula lui Simpson dacă h =
b−a d −c , r= şi x i = a + i ⋅ h, i = 0,2n şi y j = c + j ⋅ r , j = 0,2m 2m 2n
atunci I≈
[(
h⋅ r n m ∑ ∑ f x 2 i − 2 , y 2 j − 2 + f x 2 i − 2 , y 2 j + f x 2 i , y 2 j − 2 + f x 2i , y 2 j 9 i =1 j =1
) ( ) ( ) ( ) + 4( f (x 2i , y 2 j −1 ) + f (x 2i −1 , y 2 j ) + f (x 2i − 2 , y 2 j −1 ) + f (x 2i −1 , y 2 j − 2 )) + 16 f (x 2i − 1 , y 2 j − 1 ) ] .
1 Exemplu: Fie funcţia f ( x , y ) = . Să se calculeze Exemplu x⋅ y
4.4 2.6
∫ ∫ f ( x, y )dxdy , 4 2
apoi să se calculeze integrala folosind cele doua formule pentru diferite valori ale lui n şi m şi să se compare rezultatele.
Rezolvare în Scilab function f=f(x,y) f=1/(x*y); endfunction function I=IntegrareCubTrapeze(a,b,c,d,n,m) h=(b-a)/n; r=(d-c)/m; for i=1:n+1 x(i)=a+(i-1)*h; end for j=1:m+1 y(j)=c+(j-1)*r; end for i=1:n+1 for j=1:m+1 M(i,j)=f(x(i),y(j)) end end s1=0; for i=2:n+1 for j=2:m+1 s1=s1+(M(i,j)+M(i-1,j)+M(i,j-1)+M(i-1,j-1));
end end I=h*r*s1/4 endfunction function II=IntegrareCubSimson(a,b,c,d,n,m) h=(b-a)/(2*n); r=(d-c)/(2*m); for i=1:2*n+1 x(i)=a+(i-1)*h; end for j=1:2*m+1 y(j)=c+(j-1)*r; end for i=1:2*n+1 for j=1:2*m+1 M(i,j)=f(x(i),y(j)) end end s1=0; for i=2:n+1 for j=2:m+1 s2=M(2*i-3,2*j-3)+M(2*i-3,2*j-1)+M(2*i-1,2*j3)+M(2*i-1,2*j-1); s3=M(2*i-1,2*j-2)+M(2*i-2,2*j-1)+M(2*i-3,2*j2)+M(2*i-2,2*j-3); s1=s1+s2+4*s3+16*M(2*i-2,2*j-2); end end II=h*r*s1/9 endfunction IntegrareCubTrapeze(4,4.4,2,2.6,2,2) IntegrareCubSimson(4,4.4,2,2.6,2,2) ans = ans =
0.0250060 0.0250060
MODELARE ŞI SIMULARE ÎN SCILAB
Una din problemele fundamentale din multe domenii ştiinţifice şi inginereşti este problema modelării şi simulării. Scilab-ul furnizează o serie de mijloace pentru dezvoltarea şi simularea diferitelor tipuri de modele. 1. Problema Cauchy pentru ecuaţii diferenţiale ordinare Problema
Cauchy
pentru
ecuaŃii
diferenŃiale
constă
în
y ′( x ) = f ( x , y( x )) . y( x 0 ) = y 0
determinarea unei funcŃii y(x) astfel ca
Se imparte intervalul [x0, x0+L] pe care se caută soluŃia în subintervale egale de lungime h: x00, în funcŃie de valorile (xk,yk), k
y(x care se face să fie cât mai mică. Adică, dacă y( xk)=yk pentru k
diferenŃa |y(xi)-yi| este
mică pentru h mic. Valorile soluŃiei în punctele xi se aproximează prin numerele yi care se determină cu formulele de mai jos.
1.1. 1.1. Metode uniuni-pas
Se numesc astfel metodele în care pentru calculul lui yi se folosesc doar valorile pentru x şi y la pasul i-1. Avem y ( x + h ) = y ( x ) + ∫ y' ( x + t )dt = h
0
h ∫0 f ( x , y( x + t ))dt .
(1)
In funcŃie de cum se aproximează integrala se obŃin diverse metode pentru estimarea lui y(x+h).
1.1.1. 1.1. Metoda Euler: 1. h
y ( x + h) = y( x ) + ∫ y' ( x + t )dt =
(2)
0
( ) = y(x ) + hf (x, y(x )) + O(h )
= y ( x ) + y' ( x )h + O h
2
2
adică
( )
y ( x i −1 + h) = y ( x i −1 ) + f ( xi −1 , y ( xi −1 )) h + O h 2 .
Putem calcula aproximativ atunci pe yi astfel: y i = y i −1 + h ⋅ f ( x i −1 , y i −1 ), i = 1, n
(3)
In acest caz s=2, iar metoda se numeşte metoda lui Euler.
1.1.2. 1.2. Metoda predictor1. predictor-corector Să calculăm integrala (1) prin metoda trapezelor. Găsim : y( x + h) = y ( x ) +
( )
h ( y' ( x ) + y' ( x + h)) + O h3 = 2
( )
h = y ( x ) + ( f ( x , y ( x )) + f ( x + h, y ( x + h))) + O h3 2
(4)
Prin urmare această metoda ne conduce la: y i = y i −1 +
h ( f ( xi −1 , yi −1 ) + f ( xi , yi )) 2
(5)
Formula (5) conŃine pe yi implicit, dar pentru h mic, membrul II din (5) este o contracŃie ca funcŃie de yi, deci se poate determina yi din (5) prin aproximaŃii succesive: y p = y i −1 + hf ( xi −1 , yi −1 ) predictor i c h p yi = yi −1 + f ( xi −1 , yi −1 ) + f xi , yi corector 2 y p = y c noul predictor i i
(
(
))
(6)
Prin repetarea ultimelor doua ecuaŃii din (5) până | y ip − y ic | devine suficient de mic se ajunge la yi din (5) cu precizia dorită. Metoda se
numeşte
predictor-corector
şi
aproximează
soluŃia
ecuaŃiei
diferenŃiale mai bine ca metoda Euler: s=3.
1.1.3. 1.3. Metoda Euler perfectionată: 1. Daca în (4) locuim y(x+h) cu y* astfel ca |y(x+h) |y(x+h)-y*|=O(h2) atunci (4) se mai scrie: y ( x + h) = y ( x ) +
(
(
)) ( )
h f ( x , y ( x )) + f x + h, y* + O h3 2
(7)
Acest lucru poate fi realizat dacă luăm y*=y(x)+hf(x,y(x)), (conform cu (2)), caz în care (7) devine: y ( x + h) = y ( x ) +
( )
h ( f ( x , y( x )) + f ( x + h, ( y ( x ) + hf ( x , y ( x )))) + O h3 2
(8)
De aici deducem urmatoarea formulă de calcul aproximativ: y i = y i −1 +
h ( f ( xi −1 , yi −1 ) + f ( xi , yi −1 + h ⋅ f ( xi −1 , yi −1 )), i = 1, n 2
(9)
cunoscuta sub numele de metoda Euler perfecŃionată. Eroarea la un pas este O(h3).
1.1.4 1.4 Metode de tip Runge1. Runge-Kutta Aceste metode constau îin a determina o aproximare pentru
y(x+h) astfel: a) Se determină k1, k2,..kq prin: k1 (h) = hf ( x , y ( x )) k 2 (h) = hf ( x + α 2 h, y ( x ) + β 21k1 (h)) ....... kq = hf x + α q h, y ( x ) + βq ,1k1 (h) + βq ,2 k 2 (h) + .. + βq ,q −1kq −1 (h)
(
(10)
)
b) Se determină z(h) ca aproximare pentru y(x+h) prin: z (h) = y ( x ) + p1k1 (h) + p2 k 2 (h) + ... + pq kq (h)
(11)
CoeficienŃii α,β, p se determină din condiŃia ca r(h)=y(x+h)r(h)=y(x+h)-
z(h)=O(hs) pentru un s cât mai mare. Se gaseşte:
a) pentru q=2, s=3 rezultă α 2 =
1 1 β2,1 = p1 = t p2 = 1 − t . Pentru 2t 2t
t=3/4, y(x) inlocuit cu yi-1 si z(h) cu yi+1 se obŃin formulele: k1 = hf ( x , yi −1 ) 2 2 k 2 = hf x i −1 + h, yi −1 + k1 3 3 1 3 yi = yi −1 + k1 + k 2 4 4
(12)
numite formule Runge-Kutta de ordin 2 1 2
1 2
ii) pentru q=4, s=5 se gaseşte de exemplu α 2 = ,β21 = , α3 =
1 1 ,β31 = 0,β32 = 2 2
α 4 = 1,β41 = 0,β42 = 0,β43 = 1 ,
ceea ce conduce la
formulele: k1 = hf ( xi −1 , yi −1 ) k = hf x + h , y + k1 i −1 i −1 2 2 2 h k2 k3 = hf xi −1 + , yi −1 + 2 2 k4 = hf ( xi −1 + h, yi −1 + k3 ) y = y + 1 (k + 2k + 2k + k ), i = 1, n i −1 1 2 3 4 i 6
(13)
cunoscute ca formulele Runge-Kutta de ordinul 4. In ce priveşte convergenŃa lui yi la y(xi) avem teorema: Teorema. Teorema Fie f(x,y) de clasa C1 si |fy(x,y)|
[x0,x0+L] divizat în părŃi egale. Fie |y(xi+1)-yi|
(
max | y ( x i ) − yi |≤ e LC C1h s −1 + iδ + r0
L 0≤ i ≤ h
)
(14)
Rezolvare cu funcţia ScilabScilab-ului: ului exista funcŃii de rezolvare a ecuaŃiilor diferenşiale, de exemplu ode – rezolvarea de ecuaşii diferenşiale ordinare
y=ode(y0,t0,t,f) [y,w,iw]=ode([type],y0,t0,t [,rtol [,atol]],f [,jac] [,w,iw]) [y,rd,w,iw]=ode("root",y0,t0,t [,rtol [,atol]],f [,jac],ng,g [,w,iw])
y=ode("discrete",y0,k0,kvect,f) Parametri: y0: y0 vector real sau matrice (condiŃiile iniŃiale). t0: t0 scalar real (timpul iniŃial). t: vector real (timpuri la care soluŃia e calculată). f: extern (funcŃie sau şir de caractere sau listă). type: type unul din şirurile de caractere: "adams" "stiff" "rk" "rkf" "fix" "discrete" “roots" rtol,atol: rtol,atol constante reale sau vectori reali de aceeaşi dimensiune ca şi y. jac: jac extern (funcŃie sau şir de caractere sau listă). w,iw: w,iw vectori reali. ng: ng întreg. g: extern (funcŃie sau şir de caractere sau listă). k0: k0 întreg (timp iniŃial). kvect : vector de întregi. Pentru mai multe detalii să se apeleze help(“ode”), şi să se vadă funcŃiile date la See Also. Exemplu y′( x ) = x * y Se cere valoarea lui y pentru x=1. y (0) = 1
Fie
SoluŃie: function f=f(x,y) f=x*y; endfunction function y=Euler(n,x,xx,y)
h=(xx-x)/n; for i=1:n do y=y+h*f(x,y); x=x+h; end endfunction Euler(10,0,1,1) function y1=EulerPerfectionata(n,x,xx,y1) h=(xx-x)/n; for i=1:n do y1=y1+h/2*(f(x,y1)+f(x+h,y1+h*f(x,y1))); x=x+h; end endfunction EulerPerfectionata(10,0,1,1) function y2=RungeKutta(n,x,xx,y2) h=(xx-x)/n; for i=1:n do k1=h*f(x,y2); k2=h*f(x+h/2,y2+k1/2); k3=h*f(x+h/2,y2+k2/2); k4=h*f(x+h,y2+k3); y2=y2+1/6*(k1+2*k2+2*k3+k4); x=x+h; end endfunction RungeKutta(10,0,1,1) ode(1,0,1,f) ans =
1.5471104
ans =
1.6478813
ans =
1.648721
ans =
1.6487224 1.6487224
Laboratorul al 88-lea Funcţii de intrare/ieşire a datelor Implicit, în Scilab rezultatele apar în fereastra de comandă, iar graficele într-o fereastră grafică. Sunt situaţii în care rezultatul se doreşte salvat într-un fişier sau să fie scos la imprimantă, de asemenea datele de intrare se pot prelua dintr-un fişier sau se pot introduce de la tastatură. In acest capitol se vor descrie modalităţi de intrare ieşire a diferitelor tipuri de date. Afişarea variabilelor Se defineşte variabila urmatoare: -->a=[%pi,4,%inf] Tabel 1. Afişarea variabilelor Modalităţi de afisare Exemplu Afişare in fereastra de comanda a = (implicit) 3.1415927 4. Inf Afişare utilizând disp -->disp(a) 3.1415927 4. Inf Afişare utilizând print -->print(%io(2),a) a = 3.1415927 4. Inf Scrierea valorii variabilei a intr-un -->print('a.txt',a) fişier, de exemplu, a.txt Print poate fi folosit pentru a forţa -->function y=g(x); a=2*%pi, afişarea unei variabile in interiorul unei print(%io(2),a), y=x+a, endfunction funcţii -->g(2) a = 6.2831853 Utilizarea interfeţelor grafice In Scilab există o serie de instrucţiuni pentru ferestre de dialog predefinite, dintre acestea putem aminti • x_message; x_message deschide o ferestră de mesaj • x_choose; x_choose deschide o fereastră ce conţine o listă de elemente • x_dialog; x_dialog deschide o fereastră de dialog pentru a permite utilizatorului să dea un răspuns • x_choices; permite alegerea unei opţiuni dintr-o listă Exemplul 1 -->x_message(['Ce -->x_message(['Ce ati vrea';'sa invatati?'])
optiune=x_choose(['Windows','Linux','Alt optiune=x_choose(['Windows','Linux','Alt sistem de operare'],'Sisteme de operare') ObservaŃie Observa ie: ie Dacă se alege o opŃiune din listă apare, ca răspuns, numărul de ordine lin listă a opŃiunii alese optiune = 2.
-->linux=['Red -->linux=['Red Hat Linux','Mandriva Linux','Mepis Linux','Mepis Linux']; -->l1=list('Sistemul -->l1=list('Sistemul de operare Linux',1,linux); -->l2=list('Programe',2,['MuPAD','Scilab','Maxima']); -->l2=list('Programe',2,['MuPAD','Scilab','Maxima']); -->x_choices('Linux -->x_choices('Linux si programe sub Linux',list(l1,l2)) ans = 2. 2. Se aleg opŃiunile
După ce se aleg opŃiunile se apasă Finish. Finish Exemplul 2 In exemplul următor se realizează un script în care se foloseşte o fereastră de dialog pentru introducerea a patru valori de intrare şi o ferestra de afişare a rezultatului.
text=['no pentru 532';'ne pentru 532';'no pentru 1064';'ne pentru 1064']; m=x_mdialog('Introduceti valorile indicilor de refractie',text,[' ';' ';' ';' ']); no1=evstr(m(1)); ne1=evstr(m(2)); no2=evstr(m(3)); ne2=evstr(m(4)); teta=asin(sqrt(abs((no1^(-2)-no2^(-2))/(ne2^(-2)-no2^(-2)))))*180/%pi; x_message('Unghiul de adaptare de faza este: '+string(teta)+' grade')
Citirea datelor de la tastatura si afisarea rezultatelor folosind ferestre de dialog: Citirea unei valori: valoare=evstr(x_mdialog('Introduceti o valoare: ',['v='],[' '])); Citirea unui vector: a=evstr(x_matrix('introduceti elementele unei matrice de dimensiune 1x5: ',zeros(1,5))) Citirea unei matrice: a=evstr(x_matrix('introduceti elementele unei matrice de dimensiune 3x3: ',zeros(3,3))) Afisarea unei matrice: x_matrix(‘Elementele matricei s sunt’,s)
Exemplu Să se calculeze suma elementelor unui vector, care sunt cuprinse intre valorile a şi b cu aa)&(v(i)
Observatie: Sa se modifice programul de mai sus astfel incat pentru introducerea valorilor a si b sa se foloseasca o singura fereastra de dialog Exercitii: Sa se rezolve exercitiile de mai jos folosind pentru introducerea si afisarea datelor ferestre de dialog. 1. Să se calculeze produsul elementelor diferite de zero ale unui vector. 2. Se dă şirul de numere x1, x2, …, xn. Să se determine numărul de elemente pozitive şi să se calculeze produsul lor. 3. Să se calculeze media aritmetică a elementelor unui vector care sunt mai mari decat o valoare dată. 4. Sa se determine elementul maxim al unui matrice de numere reale (se va defini o functie in editorul scilab-ului). 5. Să se calculeze produsul elementelor diferite de zero de pe diagonala principală (secundară) ale unei matrice. 6. Să se calculeze suma elementelor unei matrice situate deasupra (dedesubtul) diagonalei principale (secundare). 7. Sa se realizeze un program care calculeaza produsul a doua matrice. 8. Se dau numerele aij , i = 1,2,...,100; j = 1,2,...,50 si b j j = 1,2,...,50 . Sa se realizeze un 50
program care calculeaza:
ci = ∑ aij b j , i = 1,2,...,100 . j =1
9. Să se calculeze suma elementelor unui vector care sunt cuprinse intre valorile a şi b cu a
a i + bi , daca a i ⋅ bi < 0 ci = 2 max{a i , bi }, daca a i ⋅ bi ≥ 0
Laboratorul al 99-lea
Fişiere de date în Scilab Tabel 1. Funcţii de citire/scriere ale datelor în fişiere Scilab Procedura
Semnificaţie
[x]=read(file[x]=read(file-desc,m,n,[format]) [x]=read(file[x]=read(file-desc,m,n,k,format)
Citeşte din fişier; filefile-desc: desc şir de caractere specificând numele fişierului sau valoare întreagă specificând unitatea logică m, n: numere întregi n reprezentând dimensiunile matricei. Dacă nu se dă numărul de linii, se citeşte tot fişierul format: format şir de caractere ce specifică un format ”Fortran” k: întreg sau vector de întregi Scriere în fişier.
write(fisier,text_de_scris)
Alte instrucţiuni pentru fişiere: [unit,[err]]=file(‘open’,nume_fisier, [unit,[err]] file(‘open’,nume_fisier, Deschidere fişier: [status]) status poate fi “new” (pentru fişier nou), “old” (pentru fişier deja existent), “unknown” (când nu se ştie exact) unit este un ‘identificator’ pentru fişier, err conţine posibilele erori. file(‘close’, unit) Inchide fisierul descris de unit file(‘last’,unit) Pozitionare la sfârşit de fişier. file(‘rewind’,unit) Pozitionare la început de fişier. Pentru mai multe detalii să se apeleze apropos comanda. comanda
Exerciţii 1. Să se creeze un fişier în care să se scrie elementele a patru vectori creaţi aleator, apoi să se citească în variabila x primii doi vectori. u=file('open',TMPDIR+'/foo','unknown') for k=1:4 a=rand(1,4) write(u,a) end file('rewind',u) x=read(u,2,4) file('close',u) 2. Fie a=5, b=3. Să se scrie în fişierul Date aceste date. Apoi să se citească
conţinutul
fişierului.
Se
prelucrează
datele:
se
calculează c=a+b. Să se salveze valoarea lui c în acelaşi fişier.
Solutie: -->a=3;b=5; -->u=file('open','Date','unknown'); -->write(u,a); -->write(u,b); -->file('rewind',u) -->x=read(u,2,1); -->c=sum(x) c =
8.
3. Să se definească matricele a ∈ Μ 3×6 , b ∈ Μ 4×6 în SCILAB, să se salveze în fişierul Matrice, apoi să se citească conţinutul fişierului şi să se srie în matricea c. Să se vizualizeze cele trei matrice.
Soluţie:
a=[1 2 3 4 5 6;2 3 4 5 6 7;3 4 5 6 7 8]; b=[4 5 6 7 8 9;4 5 6 7 8 9;4 5 6 7 8 9;4 5 6 7 8 9]; u=file('open','Matrice','unknown'); write(u,a); write(u,b); file('rewind',u) c=read(u,7,6) file('close',u) In fereastra de comanda va aparea: c =! 1.
2.
3.
4.
5.
6. !
! 2.
3.
4.
5.
6.
7. !
! 3.
4.
5.
6.
7.
8. !
! 4.
5.
6.
7.
8.
9. !
! 4.
5.
6.
7.
8.
9. !
! 4.
5.
6.
7.
8.
9. !
! 4.
5.
6.
7.
8.
9. !
4. Să se calculeze: a=
10
1
∑ ( i − 7 )2
i =1 i≠7
b=
20
11
−
∏ ( )2 k−2
k =1 k≠2
Să se salveze în fişierul Valori variabilele a şi b, apoi să se citească conţinutul fişierului.
Solutie: s=0; for i=1:10 if i~=7 s=s+1/(i-7); end end p=1; for i=1:20
if i~=2 p=p+11/((i-2)^2); end end u=file('open','Valori','unknown'); write(u,s); write(u,p); file('rewind',u) r=read(u,1,2) file('close',u) In fereastra de comandă va apărea:
r =! - 0.6166667
29.499825 !
5. Să se realizeze un program în SCILAB care: încarcă în matricea Matrice apoi recalculează fiecare element a conţinutul fişierului Matrice, al matricei (mai puţin cele de pe frontieră) ca media aritmetică a elementelor vecine. Procedeul se va repeta până când diferenţa dintre valoarea nou calculată şi cea anterioară este mai mică decât un epsilon dat sau după maxim 1000 de iteraţii.
Soluţie: u=file('open','C:\Program Files\scilab-3.0\matrice.txt', 'unknown') ;a=read(u,4,5);file('close',u); [m,n]=size(a); an=zeros(m,n);k=0;ermax=1; while (ermax>10^(-10))&(k<1000) for i=2:m-1 for j=2:n-1 an(i,j)=(a(i-1,j)+a(i+1,j)+a(i,j-1)+a(i,j+1))/4; er=abs(an(i,j)-a(i,j)); if er>ermax
ermax=er; end a(i,j)=an(i,j); end end k=k+1; end a;x=1:4;y=1:5; contour2d(x,y,a,10); xtitle('Liniile de contur!') In urma lansării în execuţie a scriptului de mai sus, în fereastra de comandă va aparea matricea a, iar în ferestra grafică vor apărea liniile de contur Linii le de contur! 5.0
0.818
1.636 2.455
4.6
3.273 4.091
4.2 5.727 3.8
7.364 4.909 6.545 4.909
3.4
5.727
3.0
2.6
2.2 3.273 1.8
4.091
2.455
1.4
1.636
1.0 1.0
1.4
1.8
2.2
2.6
3.0
3.4
3.8
4.2
Laboratorul al 1010-lea
Modelare si simulare in SCICOS 1. Introducere Simularea în SCILAB se realizează cu opţiunea Scicos din meniul Applications care este o aplicaţie a programului cu facilităţi importante, permiţând obţinerea unor modele precise ale unor sisteme complexe. Scicos este un utilitar Scilab gratuit (opensource) inclus în pachetul Scilab, ce oferă multe din functionalităţile puse la dispoziţie de Simulink şi MATRIXx. Bazat pe un formalism deschis motivat de limbaje sincrone extins la dinamica timpului continuu, Scicos poate fi folosit pentru a modela şi simula sisteme dinamice hibride. El utilizează infrastructura de calcul a SCILAB-ului, respectiv organizarea matriceala a variabilelor. Principalul avantaj este interfaţa comodă cu utilizatorul, acesta având la dispozitie
blocuri ce realizează diferite funcţii: matematice, de conectare, de vizualizare etc. Prin interconectarea acestora, pe baza modelelor matematice
ale
sistemelor
simulate,
se
construiesc
modele
complexe. La randul lor, acestea pot fi grupate, creandu-se noi blocuri, ce pot fi în continuare interconectate. Blocurile sunt organizate în biblioteci. Un bloc Scicos poate avea 2 tipuri de intrări şi ieşiri: intrare regulară (de obicei aşezată pe margini) ieşire regulară (de asemenea pe margini) intrare de activare (de obicei aşezată sus) ieşire de activare (de obicei aşezată jos) Intrările şi ieşirile regulare sunt folosite pentru a comunica date de la un bloc la altul prin legături, intrările şi ieşirile de activare sunt conectate de legături de activare care transmit informaţia de control (activare).
Descrierea
modului
de
lansare
a
opţiunii
Scicos
şi
a
componentei bibliotecilor se va realiza considerandu-se varianta SCILAB 3.0. 2. Construcţia diagramelor Folosind Scicos din Applications Applications se pot realiza diverse diagrame şi simulări în fereastra care se deschide.
Deschiderea bibliotecilor se realizează făcând click pe opţiunea Palettes din Edit (sau clic dreapta cu mouse-ul),
Pentru realizarea diagramelor se pot alege obiecte folosind: • Sources: Sources blocuri ce reprezintă surse de semnale (sinusoidal, constant, generator de semnal, etc.).
• Sinks: Sinks blocuri de vizualizare a semnalelor, pot fi modificate rezoluţiile pe verticală şi orizontală, în funcţie de domeniile semnalelor vizualizate.
• Liniar: Liniar blocuri de calcule matematice lineare.
• NonNon-linear: linear blocuri de calcule matematice nelineare, funcţii trigonometrice, funcţii SCILAB.
• Events: Events blocuri referitoare la evenimente
• Threshold: Threshold butoane de toleranţă
• Others: Others butoane relaţionale, logice, etc.
• Branching: Branching butoane de transformare.
• Electrical Electrical: cal blocuri pentru circuite electrice
• Hydraulics: Hydraulics blocuri hidraulice
2.1. Realizarea unui model Deschiderea unei noi ferestre de modelare se poate face alegând opţiunea New al meniului Diagram. Diagram Plasarea blocurilor în noua schemă se realizează prin drag-area = tragerea acestora (apasarea butonului din stânga al mouse-ului pe blocul necesar) şi poziţionarea blocului în noua schemă. Unele blocuri au posibilitatea actualizării parametrilor, aceştia având valori implicite pentru blocurile luate din biblioteci. Facând dublu click pe fiecare bloc, se va deschide o fereastra de dialog în care se introduc noile valori ale parametrilor blocului respectiv. Interconectarea blocurilor se realizează prin unirea (cu butonul din stânga) apăsată al unei borne de ieşire a unui bloc cu o bornă de intrare a altui bloc ( urmariţi modificarea tipului
de cursor pentru a vedea cand poate fi eliberat butonul mouse-ului. Un punct de conexiune (conectare) a unei ieşiri la intrările mai multor blocuri se realizează făcând click dreapta pe prima legatură şi tragerea acesteia spre celelalte intrări. Se pot introduce instrucţiuni cu ajutorul ferestrei de dialog obţinute prin accesarea optiunii Context (clic dreapta cu mouse-ul).
După realizarea schemei bloc corespunzatoare modelului matematic se plasează blocurile de vizualizare; (cel mai frecvent Scope„Osciloscop” din biblioteca Sinks). Acestea trebuie activate (dublu
click), deschizându-se fereastra ce conţine ecranul osciloscopului, putându-se în acest moment modifica configurarea osciloscopului. Pentru a vizualiza mai multe semnale în acelaşi sistem de axe, semnalele vor fi multiplexate, conectate la intrarile unui bloc Mux, ieşirea acestuia conectându-se la un osciloscop având un singur sistem de axe. 2.2 Stabilirea parametrilor simulării După realizarea modelului se selectează parametrii simulării din fereastra Set Block properties (meniul Simulate/Setup) momentul începerii simulării, durata simularii, metoda de integrare, pasul maxim de integrare, eroarea.
In ceea ce priveste metoda de integrare, Scicos prezintă iniţial în fereastra de modificare a parametrilor simulării metoda implicit aleasă în funcţie de structura modelului. Aceasta poate fi schimbată, alegandu-se între o metoda cu pas variabil de integrare şi una cu pas fix. Metoda de integrare cu pas variabil implicit aleasă este ode45, ceea ce constituie metoda de integrare Runge-Kutta de ordinul 5, ce oferă rezultate bune pentru majoritatea modelelor de tip continuu. Metodele de integrare cu pas fix sunt variante ale celor cu pas variabil. Lansarea în execuţie se face din meniul Simulate/Run. Salvarea unui model se poate realiza cu comanda din meniul
Diagram-Save As..., specificandu-se directorul şi numele sub care va fi salvat (Nu faceţi salvări decât în directorul propriu de lucru!). Exemple Scicos din Scilab este asemănător cu simulink din mathlab, un simulator şi constructor de sistem dinamic, este putin dificil de folosit deoarece trebuie asociat semnalul de activare de timp unor anumite blocuri. O dată ce conceptul de activare de timp este manevrat restul este uşor. Scicos are interfacţa grafică uzuală pentru a alege şi a modifica dinamic sisteme de blocuri. Blocuri noi pot fi definite de utilizator in C, Fortran sau Scilab. De fapt Scicos nu este doar o versiune de diagrame de blocuri a unui sistem dinamic de
simulare, este construit în jurul unui formalism de sisteme hibride, precum activarea timpului, etc. Folosind exemple simple, se va arăta cum modelele hibride (timp continuu sau discret) pot fi construite în acest formalism. Se va arăta, de asemenea, că acest formalism este strans legat de limbajele sincrone. Scicos include un copilator, acesta foloseşte descrierea modelului, de obicei compilat de editorul Scicos, pentru a construi tabele de program (orar) care apoi pot fi folosite de simulator şi de funcţii de generare de cod. Include, de asemenea,
un simulator, simulatorul
Scicos foloseste tabele de orar şi alte informaţii puse la dispoziţie de compilator pentru a rula simularea. Simulatorul este de natură hibridă, adică are de-a face cu sisteme de timp sau eveniment discrete sau continui. Pentru partea de timp continuu foloseste o parte de rezolvare ODE LSODAR sau DAE DASKR, depinzând de natura sistemului de timp continuu considerat, şi include şi un generator de cod: Scicos poate genera cod C pentru realizarea comportamentului unor subsiteme. Aceste subsisteme nu trebuie să includă componente de timp continuu.
Exemplul 1
Un script simplu în Scicos pentru simularea răspunsului sinusoidalei într-un circuit, capturează cateva lucruri esenţiale în folosirea la inceput a Scicosului. Diagrama de blocuri şi rezultatul în modul de reprezentare Scicos este prezentat în aplicaţia urmatoare.
După ce se contruieşte modelul fixând blocurile şi legăturile corespunzatoare se definesc proprietăţile fiecarui bloc în parte, astfel pentru blocurile: • sinusoid generator
• num(s)/den(s)
• scope
• ceas
Se lansează simularea
şi se observă efectele: apare graficul din fereastra de mai jos şi în fereastra de comandă apare valoarea corespunzătoare blocului de calcul: 0.890625.
Tema: Tema să se modifice schema modelului conform imaginii de mai jos apoi să se observe efectele.
