------------------------------------------ APROKSYMACJA function [yi,xi,y_est,x]=aproksymacja_trygonometryczna(yi,p) ni=length(yi); xi=2*[0:ni-1]*pi/ni; xi=xi(:); yi=yi(:); r=(max(xi)-min(xi))/10; x=linspace(min(xi)-r,max(xi)+r,1000); x=x(:); if mod(ni,2)==0 pmax=ni/2; else pmax=(ni-1)/2; end; p=min(pmax,p); wsp_cos=NaN(1,p+1); wsp_sin=NaN(1,p+1); for m=0:p wsp_cos(m+1)=2*sum(cos(m*xi).*yi)/ni; wsp_sin(m+1)=2*sum(sin(m*xi).*yi)/ni; end; wsp_cos(1)=wsp_cos(1)/2; if p==pmax && mod(ni,2)==0 wsp_cos(end)=wsp_cos(end)/2; end; yi_est=zeros(size(xi)); y_est=zeros(size(x)); for m=0:p yi_est=yi_est+wsp_sin(m+1)*sin(m*xi)+wsp_cos(m+1)*cos(m*xi); y_est=y_est+wsp_sin(m+1)*sin(m*x)+wsp_cos(m+1)*cos(m*x); end; ------------------------------------------ BISEKCJA function x=bisekcja(wzor,a,b,fx) x=(a+b)/2; ya=wylicz(wzor,a); yx=wylicz(wzor,x); while abs(yx-fx)>0.00001 if (yx-fx)*(ya-fx)<0 b=x; else a=x; end; x=(a+b)/2; yx=wylicz(wzor,x); end; function y=wylicz(wzor,x) eval(['y=',wzor,';']); ------------------------------------------ CAŁKOWANIE function s=calkowanie_trap(wzor,a,b,n) h=(b-a)/n; s=0; x=a; for i=1:n pole=h*(wylicz(wzor,x)+wylicz(wzor,x+h))/2; x=x+h; s=s+pole; end; function y=wylicz(wzor,x) zbudowany_napis=['y=',wzor,';']; eval(zbudowany_napis); ------------------------- CAŁKOWANIE 2 function s_new=calkowanie_trap2(wzor,a,b) % s=calkowanie_trap2('x^2',3,7) n=2; s_old=calkowanie_trap(wzor,a,b,n); s_new=calkowanie_trap(wzor,a,b,2*n); while abs(s_old-s_new)>0.00001 s_old=s_new; n=2*n; s_new=calkowanie_trap(wzor,a,b,2*n); end; ------------------------------------------ INT LAGRANGE function [y,a]=interp(xi,yi,x) xi=xi(:); yi=yi(:); n=length(xi); X=[]; for i=0:n-1 X=[X,xi.^i]; end; a=X\yi; y=x*0; for i=1:n y=y+a(i)*x.^(i-1); end; ------------------------- PRZYKŁAD function przyklad_interp xi=[-2,1,4,5]; yi=-3+4*xi-2*xi.^2+0.4*xi.^4; x=linspace(min(xi)-2,max(xi)+2,1000); [y,a]=interp(xi,yi,x); plot(xi,yi,'ro'); hold on; %trzymaj wykres plot(x,y); hold off; %zwolnij wykres ------------------------------------------ NEWTON function x=met_Newtona(wzor,x) % x=met_Newtona('(x+3)^4-10000',5) fx=wylicz(wzor,x); while abs(fx)>0.00001 fp=pochodna(wzor,x); delta=fx/fp; x=x-delta; fx=wylicz(wzor,x); end; function y=pochodna(wzor,x) %oszacowanie wart poch w punkcie x eps=0.0001; y=(wylicz(wzor,x+eps)-wylicz(wzor,x-eps))/(2*eps); function y=wylicz(wzor,x) eval(['y=',wzor,';']); ------------------------------------------ ROZKŁAD LU function [L,U,det]=rozkladLU(A) n=size(A,1); U=A; L=eye(n); %Liczenie L i U for j=1:n-1 for i=j+1:n m=U(i,j)/U(j,j); U(i,:)=U(i,:)-m*U(j,:); L(i,j)=m; end; %Liczenie wyznacznika diagU = diag(U); s = size(diagU); det = 1; for i = 1:s det = det*diagU(i); end; end; ------------------------------------------ WARTOŚCI WŁASNE function wlasnosc(A) [V,D] = eig(A); disp('Najwieksza wartosc wlasna: '), disp(max(eig(A))); disp('Wektor: '), disp(V(:,1)); end; ------------------------------------------ MÓJ GAUSS-SEIDEL function x=gauss_seidel(A, b, n) s = size(A,1); x = zeros(s,1); # Wartosc poczatkowa x M = A-diag(diag(A)); # Macierz do obliczen - nie korzystamy z wartosci po przekatnej M = L + U invA = diag(A).^-1; # Macierz przekątnej - potrzebna do obliczeń; Już po odwróceniu! D-1 z przykładu for i=1:n xold = x; # Przechowywanie poprzedniej wartości x for j=1:length(x) # Równanie Gaussa - Seidela x(j) = invA(j)*(b(j) - M(j,:)*x); end end disp(M); end ------------------------------------------ STACHURSKIEGO function odpowiedzi = inter_jacoby(a, b) D = diag(diag(a)); L = a - D; B = inv(D)*L*-1; C = inv(D)*b'; odpowiedzi = b*0; blad = 1; while abs(blad) > 0.0001 odpowiedzi = odpowiedzi * B' + C'; wynik = (odpowiedzi * a') - b; blad = sqrt(sum(wynik.^2)); end end function odpowiedzi = inter_sor(a, b, w) L = triu(a, 1); U = tril(a, -1); D = diag(diag(a)); Bw = inv(D+w*L) * ((1-w)*D-w*U); c = w*(inv(D+w*L)*b'); odpowiedzi = b*0; blad = 1; while abs(blad) > 0.0001 odpowiedzi = odpowiedzi * Bw' + c'; wynik = odpowiedzi * a' - b; blad = sqrt(sum(wynik.^2)); end end