lv_blatt2.m

LU－Zerlegung und Vektorisierung

%Blatt2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,U]=my_LU_opt(A)
% Berechnet eine LU-Zerlegung ohne Pivotizierung (Gauss-Elimination)
% 
% Eingabe: A nxn regulaer Matrix
% Ausgabe: L nxn-untere Dreiecksmatrix
%          U nxn-obere Dreiecksmatrix
% aus dem Buch: Gramlich, Werner Numerische Mathematik mit Matlab
%               b360/ Gram

format long
[m,n]=size(A);

if (m ~= n )
    error('Matrix A ist nicht quadratisch')
end;

for k=1:n-1
    if (A(k,k) == 0)
        error('Pech! Es geht nicht mehr mit dieser Matrix!')
    end;
    A(k+1:n,k)=A(k+1:n,k)/A(k,k);
    A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
end;
L=eye(n,n)+tril(A,-1);
U=triu(A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 function [X]=my_SystemSolver(A,B,Which) 

 % Loest ein System Ax=b 
 %               mit Hilfe LUZerlegung_opt.m und lu.m
 % Eingabe: A      Regulaere nxn Matrix 
 %          b      rechte Seite (nxm Matrix)
 %          Which  =0, falls my_LU_opt.m benutzt wird
 %                 =1, falls lu.m benutzt wird
 % Asgabe:  X Loesung von AX=B (nxm Matrix)
 
 if (cond(A)==Inf) 
     error('Matrix A ist nicht regulaer')
 end;

 [n,m]=size(B);
 X=zeros(n,m);
 Y=zeros(n,m);
  
 % LU-Zerlegung
 if (Which==0)
   [L,U]=my_LU_opt(A);
 else
     [L,U,P]=lu(A);
     B=P*B;
 end;
 
 % Loesung des unteren Dreieckssystems LY=L(UX)B
 if (cond(L)==Inf) 
     error('Matrix L ist nicht regulaer')
 end;
 
 Y=ones(n,1);
 Y(1)=B(1);
 for j=2:n
       Y(j)=(B(j) - L(j,1:j-1)*Y(1:j-1))/L(j,j);
 end;
       
 % Loest das obere Dreiecksystem UX=Y
 if (cond(U)==Inf) 
     error('Matrix U ist nicht regulaer')
 end
  
 X=ones(n,1);
 X(n)=Y(n)/U(n,n);
 for j=n-1:-1:1
       X(j)=(Y(j)-U(j,j+1:n)*X(j+1:n))/U(j,j);
 end; 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
   function [X]=my_TridiagSolver(A,B)
% Eingabe: tridiagonale Matrix A, ein Vektor B mit Ax = B
% Ausgabe: Loesung von Ax=B


[L,U,F] = tridiag_LU(A);
n=length(A(1,:));
 
Y=ones(n,1);
Y(1)=B(1);
for j=2:n
    Y(j)=B(j)-L(j-1)*Y(j-1);
end
 

X=ones(n,1);
X(n)=Y(n)/U(n);
for j=n-1:-1:1
    X(j)=(Y(j)-F(j)*X(j+1))/U(j);
end
     


function [el,u,F] = tridiag_LU(A)
% Eingabe: tridiagonlae Matrix A
% Ausgabe: Vektoren el, u, und F in der Form des Blatts

n=length(A(1,:));
 
 D = diag(A);
 F = diag(A,1);
 E = diag(A,-1);
 
 
 u(1) = A(1,1);
 for j=2:n
     el(j-1) = E(j-1)/u(j-1);
     u(j) = D(j) - el(j-1)*F(j-1);
 end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%