Trebuie piratat
https://sites.google.com/site/pascanraisa/calcul-numeric---seriile-33-34
Script: Home > New Script
Trebuie sa dam run la acel script si cand ne intreaba daca sa facem din el workspace sau ceva de genul, dam da.
Command History:
Home > Layout > History > Docked
clear all
Comentariu: %
Sectiune: %% enter
Rulam sectiune:
Editor > Run Section
fprintf(...)
warning( fprintf(...) )
disp(MyMatrix)
floor(x)
ceil(x)
sign(x) % 1 sau -1 sau 0
func2str(f) % function to string
diff( f(x), x ) % Derivata lui f in x, daca f e expresie in functie de x! (cu syms, vezi mai jos)
simplify(f) % Reduce termeni dintr-o expresie simbolica f
a = 2
b = 3
c = a + b
v = [1 2 3]
max(v) / min(x)
dot(v, w) Produs scalar
v = 1:2:10 De la 1 la 10 cu pasul 3
v = linspace(low, high, n) n elemente
M = [1 2 3; 4 5 6; 7 8 9]
M(1, 1) Elementul 1, 1 din variabila a
M(1, :) Prima linie, toate coloanele
M(1, 2:3) De la 2, pana la 3, inclusiv
M([1 3], [1 3])
A * B
A ^ 2 (A * A)
det(A)
inv(A)
A .* B Element cu element, ori elementul respectiv din b
A .^ 3 Ridicarea la putere a fiecarui element din A
A([1 2],:) = A([2 1],:)
fplot(f, [-5, 5])
plot(x, y, 'o r -') x = vector pentru x-axis, y la fel, o = cerculete, r = red
plot(x, y, 'o r', 'linewidth', 3, 'proprietate', 'valoare')
Vezi help plot pentru mai multe optiuni
2 grafice in acelasi grafic:
plot(...)
hold on
plot(...)
legend('grafic 1', 'grafic 2')
grid on % Face display si la grid in grafic
xlabel('axa ox')
ylabel('axa oy')
title('graficul meu')
plot(x1, y1, x2, y2, ...)
fill(xuri, yuri, color)
if ....
...
elseif ...
...
end
for i = 1:length(v)
...
end
f = @(x) x + 2
f = inline('x+2', 'x')
New > Function
function [y] = myfunc(x)
...
end
syms x y
f = x^3 - 6 * x^2
g = 8 * x^3
f + g
Astea merg cu orice, cum ar fi determinatn, inversa, etc
Substitutia se face asa:
subs(f, {x}, {2})
Si calculeaza rezultatul
Ca sa convertim o expresie cu necunoscute intr-o functie, facem:
F = matlabFunction(f, 'vars', {x})
F = matlabFunction(x^3 + x^2, 'vars', {x})
Daca scriem 'help bisectie' in command window, ne apare asta
function [x_aprox] = bisectie(f, domeniul, epsilon, n_max_iteratii, message)
% Functia de bisectie
%
% Synopsis: r = bisectie(f, domeniul)
% r = bisectie(f, domeniul, epsilon)
% r = bisectie(f, domeniul, epsilon, n_max_iteratii)
% r = bisectie(f, domeniul, epsilon, n_max_iteratii, message)
1: Functia 'myfunc' trebuie sa fie intr-in file numit 'myfunc,m'
2: Si obviously trebuie sa stea in acelasi folder ca scriptul
3: Trebuie ca folderul sa fie deschis in matlab! (run la script > change folder)
MetodaBisectiei(f, left, right)
if not hasRoot(f, left, right): exit
middle = right - (right-left)/2
if hasRoot(f, left, middle):
MetodaBisectiei(f, left, middle)
else:
MetodaBisectiei(f, middle, right)
Newton(f, fDerivat, x)
newX = x - f(x) / fDerivat(x)
Newton(f, fDerivat, newX)
Matlab:
function [xaprox] = Newton(f, df, x0, epsilon, nmaxiter)
x(1) = x0
k = 1
while true
k = k + 1;
x(k) = x(k-1) - f(x(k-1)) / df(x(k-1))
if k > nmaxiter
fprintf('Am ajuns la numarul maxim de iteratii\n')
break
end
if abs(x(k) - x(k-1)) / abs(x(k-1)) < epsilon
break
end
end
xaprox = x(k)
end
function [xaprox] = Secanta(f, a, b, epsilon)
% Pornim cu x1 si x2 aleatoare pe intervalul [a, b]
x(1) = a + (b-a) * rand(1)
x(2) = a + (b-a) * rand(1)
k = 2
while abs(x(k) - x(k-1)) / abs(x(k-1) > epsilon
k = k + 1
x(k) = x(k-2) * f(x(k-1)) - x(k-1) * f(x(k-2)) / (f(x(k-1)) - f(x(k-2)));
if x(k) < a | x(k) > b
clear x
x(1) = a + (b-a) * rand(1)
x(2) = a + (b-a) * rand(1)
k = 2
end
end
xaprox = x(k)
end
SubsDesc(m, v):
x[n] = v[n]/m[n][n]
suma = x[j]*m[i][j] + .. // Faci substitutia gen
x[i] = (b[i] - suma) / m[i][i]
Matlab:
function [x] = SubsDesc(A, b)
n = siz(A, 1)
x = []
for i = 1:n
if A(i,i) == 0
error("Error")
return
end
x(n) = 1/A(n,n) * b(n)
for k = n-1:-1:1
sum = 0
for j = k+1:n
sum = sum + A(k,j)*x(j)
end
x(k) = 1/A(k,k) * (b(k) - sum)
end
end
function [x]= Gauss(A, b, metoda, toleranta, verbose)
verbose = 1;
toleranta = 10^2 * eps;
metoda = 'FP'; % Sau 'PP' sau 'PT'
fprintf('%d',size(A,1))
fprintf('%d',size(b))
A = [A, b]; % A devine A extinsa (cu bara deasupra) cu b in dreapta
n = size(A, 1);
index = linspace(1,n,n)
for k = 1:n-1
switch metoda
case 'FP'
indice_pivot = -1
for i = k:n
if abs(A(i,k)) > toleranta
indice_pivot = i;
break
end
end
if indice_pivot == -1
fprintf('Sistem incompatibil au compatibil nedeterminat');
return
end
case 'PP'
max = -1;
indice_pivot = -1;
for i = k:n
if abs( A(i,k) ) > max
max = abs( A(i,k) );
indice_pivot = i;
end
end
if max <= toleranta % Daca max e aproape 0
fprintf('Sistem incompatibil au compatibil nedeterminat');
return
end
case 'PT'
max = abs(A(k,k));
indice_pivot = k;
indice_pivot_2 = k;
for i = k:n
for j = k:n
if abs(A(i,j)) > max
max = abs(A(i,j));
indice_pivot = i;
indice_pivot_2 = j;
end
end
end
if max <= toleranta
fprintf('Sistem incompatibil au compatibil nedeterminat');
return
end
end %endswitch
if indice_pivot ~= k
A([indice_pivot,k], :) = A([k,indice_pivot], :) % Swap(Linia1, Linia2) basically
end
if metoda == 'PT'
A(:, [indice_pivot_2, k]) = A(:, [k, indice_pivot_2])
index([indice_pivot_2, k]) = index([k, indice_pivot_2])
end
for l = k+1:n
mlk = A(l,k) / A(k,k)
A(l,:) = A(l,:) - mlk * A(k,:)
end
fprintf('Matricea extinsa curenta la iteratia %d', k);
format short
disp(A)
end %endfor
y = SubsDesc(A(:,1:n), A(:,n+1));
for i = 1:n
x(index(i)) = y(i)
end
end %endfunction
function [r, A_transformat] = Rang(A, tol)
% Rang - Calculeaza rangul unei matrici
% Synopsis - r = Rang(A)
% [r, A_transformat] = R(A)
% [r, A_transformat] = R(A, tol)
% Input - tol = toleranta (implicit 10^(-16)
% Output - r = rangul
% - A_transformat = matricea obtinuta in urma transformarilor
% elementare (matricea equelon?)
if nargin < 2
tol = 10^(-16);
end
[m, n] = size(A);
h = 1; % Linia, de la 1 la m
k = 1; % Coloana, de la 1 la an
r = 0; % Incepem de la 0
while h<=m && k<=n
max = 0;
pivot = -1;
for j = h:m % Gasim maximul pe coloana curenta k
if abs(A(j,k)) >= max
max = abs(A(j,k));
pivot = j; % Retinem indicele unde s-a gasit
end
end
if max <= tol
k = k + 1;
continue
end
if pivot ~= h
A([pivot,h],:) = A([h,pivot],:); % Interschimbam liniile p si h
disp("My A")
disp(A)
end
for l = h+1:m
mlk = A(l,k)/A(h,k);
A(l,:) = A(l,:) - mlk * A(h,:);
end
h = h + 1;
k = k + 1; % End loop
r = r + 1;
end
if nargout > 1
A_transformat = A;
end
end
function [L, U, W] = FactorizareLU(A, tol)
n = size(A, 1)
L = eye(n) % In
w = 1:n
for k = 1:n-1
max = 0;
pivot = -1;
for j = k:n % Gasim maximul pe coloana curenta k
if abs(A(j,k)) >= max
max = abs(A(j,k));
pivot = j; % Retinem indicele unde s-a gasit
end
end
if abs(A(p,k)) <= tol
error("ERROR: A nu admite factorizare LU!");
return;
end
if p ~= k
A([p,k],:) = A([k,p],:)
w([p,k]) = w([k,p])
end
for l = k+1:n
L(l,k) = A(l,k) / A(k,k)
A(l,:) = A(l,K) - L(l,k) * A(k,:);
end
if k > 1 && p ~= k
% Interschimbam sublinii din matricea L
L([p,k], 1:k-1) = L([k,p], 1:k-1);
end
end
if abs(A(n,n)) <= tol
error("A nu admite factorizre LU");
end
end
function [L, x] = Cholesky(A, b)
% Cholesky - Rezolva sistemele Ax = B cu matricea A simetrica si
% pozitiv definita
% Synopsis
%
% L = FactCholesky(A)
% [L, x] = FactCholesky(A, b)
%
% Input A = Matrice simetrica si pozitiv definita
% b = Vectorul termenilor liberi
%
% Output L = Matrice inferior triunghiulara cu elementele de pe
% diagonala principala strict pozitive (l(k,k) > 0)
% A = L * L'
% x = Solutia sistemului
alpha = A(1,1); % Prima coloana din matricea L, expresia de sub radical
if alpha <= 0
warning("Matricea A nu e pozitiv definita");
return;
end
n = size(A,1);
L = eye(n);
L(1,1) = sqrt(A(1,1));
L(2:n, 1) = A(2:n, 1) / L(1,1);
for k = 2:n
alpha = A(k,k) - sum( L(k, 1:(k-1)) .^ 2);
if alpha <= 0
warning("Matricea A nu e pozitiv definita");
return;
end
L(k,k) = sqrt(alpha);
for i = (k+1):n
suma = sum( L(i, 1:(k-1)) );
suma = suma .* L(k, 1:(k-1));
L(i,k) = (A(i,k) - suma) / L(k,k);
end
end
if nargin < 2
return;
end
y = L \ b; % magie
x = L' \ y;
end
Test it
A = [
1 2 3;
2 5 8;
3 8 14
];
b = [ -5 -14 -25 ]';
L1 = Cholesky(A);
[L2,x] = Cholesky(A, b);
disp(L1);
disp(" ");
disp(L2);
function [lambda] = Jacobi1(A, Eps)
% Jacobi1 calculeaza valorile proprii ale unei
% matrici simetrice
% Synopsis lambda = Jacobi(A, Eps)
% Input: A = matrice patratica simetrica
% (valorile proprii ale unei
% matrici simetrice sunt reale)
% Eps = 10^(-16)
% Output: lambda = valorile proprii
n = size(A, 1);
modul = sqrt(sum(sum(A(1:n, 1:n).^2)) - sum(diag(A).^2));
modul
while modul >= Eps
maxim = abs(A(1, 2));
p = 1;
q = 2;
for i=1:n-1
for j=i+1:n
if abs(A(i, j)) > maxim
maxim = abs(A(i, j));
p = i;
q = j;
end
end
end
if A(p, p) == A(q, q)
teta = pi/4;
else
teta = 1/2 * (atan(2*A(p, q) / (A(q, q) - A(p, p))));
end
c = cos(teta);
s = sin(teta);
for j=1:n
if j ~= p && j ~= q
u = A(p,j) * c - A(q, j) *s;
v = A(p, j) *s + A(q, j) * c;
A(p, j) = u;
A(q, j) = v;
A(j, p) = u;
A(j, q) = v;
end
end
u = c^2*A(p, p)-2*c*s*A(p, q)+s^2*A(q, q);
v = s^2*A(p, p)+2*c*s*A(p,q)+c^2*A(q,q);
A(p, p) = u;
A(q, q) = v;
A(p, q) = 0;
A(q, p) = 0;
modul = sqrt(sum(sum(A(1:n, 1:n).^2)) - sum(diag(A).^2));
end
lambda = diag(A);
end
function [normap] = normap(A, p) cf?
%normap calculeaza norma p a matricei A
% Synopsis:
% normap = normap(A, p);
% Input:
% A = matrice patratica
% p = tipul normei (1, 2, inf)
% Output:
% normap = norma p a matricei A
n = size(A, 1);
if p == 1
for j = 1:n
v(j) = sum(abs(A(:,j)));
end
normap = max(v);
elseif p == 2
lambda = Jacobi1(transpose(A) * A, 10^(-10));
normap = max(sqrt(lambda));
elseif p == Inf
for j = 1:n
v(j) = sum(abs(A(j,:)));
end
normap = max(v);
end
function [xaprox, N] = MetodaJacobi(A, a, p, Eps)
% MetodaJacobi - metoda iterativa de rezolvare a sistemelor te tipul:
% Ax = a, A matrice patratica (ordin n) si inversabila
% si este convergenta pentru
% matrici cu proprietatea ca (raza spectrala) Rho(I(n) - A) < 1
% sau ||I(n) - A|| < 1 (asta implica si ce e mai sus, deci e mai ok)
% Unde Rho(B) = maximul valorilor proprii
% Rho(B) = max(abs(lambda(i)))
% (metodele celelalte, Gauss, etc, erau directe. Asta e iterativa)
% Input: A - matricea asociata sistemului (Ax = a), patratica, invs
% a - vectorul termenilor liberi
% p - tipul normei
% Eps - optional, eroarea de aproximare, implicit = 10^(-16)
% Output: xaprox - vectorul ce reprezinta solutia aproximativa
% N - numarul de iteratii
if nargin < 3
Eps = 10^(-16);
end
I = eye(size(A))
q = norm(I - A, p);
if q >= 1
fprintf("Metoda Jacobi nu asigura convergenta")
xaprox = [];
return
end
x = zeros(size(A, 1), 1);
x(:, 1) = zeros(size(a));
k = 1;
B = I - A;
b = a;
while true
k = k + 1;
disp(k);
x(:, k) = B * x(:, k-1) + b;
if norm((x(:, k) - x(:, k-1)), p) / norm(x(:, k-1), p) < Eps
break
end
end
xaprox = x(:, k);
N = k;
end
function [ xaprox, N ] = MetodaJacobi( A, a, p, Eps )
% MetodaJacobi este o metoda iterativa de rezolvare a sistemelor
% Ax = a si este convergenta pentru matrici cu proprietatea ca
% (raza spectrala)-> rho(In - A) < 1 (sau exista o norma a.i.
% ||In - A|| < 1)
% Input:
% A - matricea asociata sistemului ( Ax = a )
% - A patratica de ordin n, inversabila
% a - vectorul termenilor liberi
% p - tipul normei
% Eps (optional) - eroarea de aproximare
% - implicit Eps = 10^(-16)
% Output:
% -xaprox - vectorul care reprezinta solutia aproximativa
% -N - numarul de iteratii necesar pentru calculul solutiei aprox
% https://codeshare.io/294gpM
if nargin < 4
Eps = 10^(-16);
end
I = eye(size(A));
q = norm(I - A, p);
if q >= 1
fprintf("Metoda Jacobi nu asigura convergenta.")
xaprox = []
return
end
x = zeros(size(A, 1), 1);
x(:, 1) = zeros(size(a));
k = 1;
B = I - A;
b = a;
while true
k = k + 1;
x(:, k) = B * x(:, k-1) + b;
if norm((x(:, k) - x(:, k-1)), p) / norm(x(:, k-1), p) < Eps
break
end
end
xaprox = x(:, k);
N = k;
end
function [ xaprox, N ] = MetJacobirelaxata( A, a,sig, Eps )
%MetJacobirelaxata este o metoda iterativa de rezolvare a sistemelor
% de forma Ax = a cu A sim si poz definita
%INPUT: A - matricea asociata sistemului ( Ax = a )
% - A patratica de ordin n, simetrica si poz definita
% a - vectorul termenilor liberi
% sig- optional implicit sig = 2/(lam1+lamn)
% lamn cea mai mare val proprie lamn cea mai mica val proprie
% Eps (optional) - eroarea de aproximare
% - implicit Eps = 10^(-16)
%OUTPUT: -xaprox - vectorul care reprezinta solutia aproximativa
% -N - numarul de iteratii necesar pentru calculul solutiei aprox
if nargin < 4
Eps = 10^(-6);
end
if nargin < 3
lam = eigs(A);
sig = 2/ (max(lam) + min(lam));
end
B = eye(size(A)) - sig*A;
b = sig * a;
x = zeros(size(A, 1), 1);
x(:, 1) = zeros(size(a));
k = 1;
while true
k= k + 1;
x(:,k) = B * x(:,k-1) + b;
if sqrt(dot(A * (x(:,k) - x(:,k-1)),x(:,k) - x(:,k-1)))/sqrt(dot(A * x(:,k-1),x(:,k-1))) < Eps
break
end
end
xaprox = x(:, k);
N = k;
end