martes, 27 de septiembre de 2011

programa LU

function lu()
n=input('introduzca la cantidad de ecu ');
for i=1:n
    for j=1:n
        A(i,j)=input('los elementos de la matriz A ');
    end
end

for i=1:n
    b(i)=input('los elementos del vector B ');
end
U=zeros(n);
L=zeros(n);
for j=1:n
    L(j,j)=1;
end
for j=1:n
    U(1,j)=A(1,j);
end
for i=2:n
    for j=1:n
        for k=1:i-1
            s1=0;
            if k==1
                s1=0;
            else
                for p=1:k-1
                s1=s1+L(i,p)*U(p,k);
                end
            end
            L(i,k)=(A(i,k)-s1)/U(k,k);
        end
        for k=i:n
             s2=0;
           for p=1:i-1
               s2=s2+L(i,p)*U(p,k);
           end
           U(i,k)=A(i,k)-s2;
        end
    end
end
disp(A);
disp(L);
disp(U);
c=inv(L)*b';
x=inv(U)*c;
disp(x);


y se debe hacer correr es de la siguiente forma:

lu()
introduzca la cantidad de ecu 3
los elementos de la matriz A 4
los elementos de la matriz A -2
los elementos de la matriz A -1
los elementos de la matriz A 5
los elementos de la matriz A 1
los elementos de la matriz A -1
los elementos de la matriz A 1
los elementos de la matriz A 2
los elementos de la matriz A -1
los elementos del vector B 9
los elementos del vector B 7
los elementos del vector B 12


     4    -2    -1
     5     1    -1
     1     2    -1


    1.0000         0         0
    1.2500    1.0000         0
    0.2500    0.7143    1.0000


    4.0000   -2.0000   -1.0000
         0    3.5000    0.2500
         0         0   -0.9286


   -1.3077
   -0.2308
  -13.7692

programa seidel

le codigo de seidel es de la siguiente forma:

function z=seidel(a,b,Es)
D=diag(diag(a));
L=(tril(a)-D)*(-1);
U=(triu(a)-D)*(-1);
m=inv(D)*(L+U);
c=inv(D)*b';
n=length(m);
x0=zeros(n,1);
Em=99999999999;
cont=0;
s=0;
while Em>=Es
    for i=1:n
        for j=1:n
            s=s+(m(i,j)*x0(j));
        end
        x(i)=s+c(i);
        s=0;
        E(i)=abs(x(i)-x0(i));
        x0(i)=x(i);
    end
    cont=cont+1;
    for i=1:n
        Em=max(E(i));
    end
end
disp('numero de iteraciones')
disp(cont)
disp('valor del error')
disp(Em)
disp('vector x')
disp(x0)

lo mismo que jacobi A es la matriz, B es el vector, Es es el error deseado

programa jacobi

el codigo del programa es:

function jacobi(a,b,x0,Es)
D=diag(diag(a));
L=(tril(a)-D)*(-1);
U=(triu(a)-D)*(-1);
m=inv(D)*(L+U);
c=inv(D)*b;
E=99999999999;
cont=0;
n=length(m);
while E>=Es
x=m*x0+c;
cont=cont+1;
E=max(abs(x-x0));
x0=x;
end
disp('numero de iteraciones')
disp(cont)
disp('valor del error')
disp(E)
disp('vector x')
disp(x0)


para ejecuytar se tiene los siguiente
A es la matriz, B es el vector, xo es el vector nulo, Es es el error deseado