Capítulo
3: RESOLUCIÓN DE
SISTEMAS DE ECUACIONES LINEALES
Por
defecto se asume la resolución del sistema Ax=b.
Índice:
Reescalado
de sistemas lineales:
A=[1,2;1.1,2], b=[10;10.4], x=A\b
AA=[1,2;1.05,2], b=[10;10.4], xx=AA\b
Eavec=abs(x-xx), Ervec=Eavec./abs(x)
Ea=norm(Eavec,inf), Er=norm(Ervec,inf)
D1=diag([10,10]), D2=diag([1,1])
Are=D1*A*D2, bre=D1*b, xre=D2*(Are\bre)
AAre=Are;AAre(2,1)=AAre(2,1)-0.05, xxre=D2*(AAre\bre)
Se resuelve Ax=b y se descompone
A=LU siempre que sea posible.
n=length(b);x=zeros(n,1);Aa=[A,b];L=eye(n);
for k=1:n-1, %bucle por columnas
m=Aa(k+1:n,k)/Aa(k,k);L(k+1:n,k)=m; %columna de multiplicadores
Aa(k+1:n,k:n+1)=Aa(k+1:n,k:n+1)-m*Aa(k,k:n+1); %eliminacion
end, U=Aa(:,1:n); L,U %matrices L,U
x=zeros(n,1);x(n)=Aa(n,n+1)/Aa(n,n);
for i=n-1:-1:1, x(i)=(Aa(i,n+1)-Aa(i,i+1:n)*x(i+1:n))/Aa(i,i);
end, x
Se
resuelve (P*A)x=(P*b) y se descompone P*A=L*U siempre que A sea
regular.
n=length(b);x=zeros(n,1);Aa=[A,b];P=eye(n);L=eye(n);
for k=1:n-1, %bucle por columnas
[pivot,r]=max(abs(Aa(k:n,k)));r=r+k-1;%fila pivote
aux=Aa(k,:);Aa(k,:)=Aa(r,:);Aa(r,:)=aux; %pivotar filas en A
aux=P(k,:);P(k,:)=P(r,:);P(r,:)=aux; %pivotar filas en P
m=Aa(k+1:n,k)/Aa(k,k);L(k+1:n,k)=m; %columna de multiplicadores
Aa(k+1:n,k:n+1)=Aa(k+1:n,k:n+1)-m*Aa(k,k:n+1); %eliminacion
end, U=Aa(:,1:n); L,U,P %matrices L,U,P
x=zeros(n,1);x(n)=Aa(n,n+1)/Aa(n,n);
for i=n-1:-1:1, x(i)=(Aa(i,n+1)-Aa(i,i+1:n)*x(i+1:n))/Aa(i,i);
end, x
Resolución alternativa del sistema
Ax=b para la descomposición PA=LU:
y=L\(P*b), x=U\y
Ir al inicio
Gauss
con pivotaje completo. Sustitución regresiva:
Se resuelve (P*A*Q)(Qx)=(P*b) siempre que A sea
regular.
n=length(b);x=zeros(n,1);Aa=[A,b];Q=eye(n);P=eye(n);
for k=1:n-1, %bucle por columnas
[pivot,ii]=max(abs(Aa(k:n,k:n)));[pivot,jj]=max(pivot);
r=ii(jj)+k-1;s=jj+k-1; %fila y columna del pivote
aux=Aa(k,:);Aa(k,:)=Aa(r,:);Aa(r,:)=aux; %pivotar filas en Aa
aux=P(k,:);P(k,:)=P(r,:);P(r,:)=aux; %pivotar filas en P
aux=Aa(:,k);Aa(:,k)=Aa(:,s);Aa(:,s)=aux; %pivotar columnas en Aa
aux=Q(:,k);Q(:,k)=Q(:,s);Q(:,s)=aux; %pivotar columnas en Q
m=Aa(k+1:n,k)/Aa(k,k); %columna de multiplicadores
Aa(k+1:n,k:n+1)=Aa(k+1:n,k:n+1)-m*Aa(k,k:n+1); %eliminacion
end,U=Aa(:,1:n),Q,P %matrices U, Q y P
x=zeros(n,1);x(n)=Aa(n,n+1)/Aa(n,n);
for i=n-1:-1:1, x(i)=(Aa(i,n+1)-Aa(i,i+1:n)*x(i+1:n))/Aa(i,i);
end, x=Q*x
A=[3 1 1;0 8 3;1 2 9];b=[1;1;1];x0=[0;0;0];tol=1e-3;itermax=100;
L=tril(A,-1); U=triu(A,1); D=diag(diag(A));
MJ=-D\(L+U); vJ=D\b;
iter=1; x=MJ*x0+vJ; xJ=[x0,x];
while (norm(x-x0)>tol) & (iter<itermax)
x0=x; iter=iter+1; x=MJ*x0+vJ; xJ=[xJ,x];
end, xJ
A=[3 1 1;0 8 3;1 2
9];b=[1;1;1];x0=[0;0;0];tol=1e-3;itermax=100;
L=tril(A,-1); U=triu(A,1); D=diag(diag(A));
MGS=-(D+L)\U; vGS=(D+L)\b;
iter=1; x=MGS*x0+vGS; xGS=[x0,x];
while (norm(x-x0)>tol) & (iter<itermax)
x0=x; iter=iter+1; x=MGS*x0+vGS; xGS=[xGS,x];
end, xGS
A=[3 1 1;0 8 3;1 2 9];b=[1;1;1];x0=[0;0;0];tol=1e-3;itermax=100;
D=diag(diag(A));L=tril(A,-1);U=triu(A,1);MJ=-D\(L+U);
w=2/(1+sqrt(1-max(abs(eig(MJ)))^2) ) %valor de w elegido
iter=1; x=x0+w*D\(b-A*x0); xJw=[x0,x];
while (norm(x-x0)>tol) & (iter<itermax)
x0=x; iter=iter+1; x=x0+w*D\(b-A*x0); xJw=[xJw,x];
end, xJw
Barrido en búsqueda del w óptimo:
w=0.1:0.0001:1.9; N=length(w);R=zeros(N,1);L=length(diag(A));
for i=1:N
MJw=eye(L)-w(i)*D\A;R(i)=max(abs(eig(MJw)));
end, plot(w,R),[RR,pos]=min(R),w(pos)
A=[3 1 1;0 8 3;1 2
9];b=[1;1;1];x0=[0;0;0];tol=1e-3;itermax=100;
D=diag(diag(A));L=tril(A,-1);U=triu(A,1);MGS=-(D+L)\U;
w=2/(1+sqrt(1-max(abs(eig(MGS)))^2) ) %valor de w elegido
iter=1; x=(D+w*L)\(((1-w)*D-w*U)*x0+w*b); xGSw=[x0,x];
while (norm(x-x0)>tol) & (iter<itermax)
x0=x; iter=iter+1; x=(D+w*L)\(((1-w)*D-w*U)*x0+w*b);
xGSw=[xGSw,x]; end, xGSw
Barrido en búsqueda del w óptimo:
w=0.1:0.0001:1.9;
N=length(w);R=zeros(N,1);
for i=1:N
MGSw=(D+w(i)*L)\((1-w(i))*D-w(i)*U);R(i)=max(abs(eig(MGSw)));
end, plot(w,R),[RR,pos]=min(R),w(pos)
Se presupone que se han ejecutado todos los métodos iterativos anteriores para realizar un estudio de la convergencia de éstos.
xJn=sqrt(sum((diff(xJ'))'.^2));
NJ1=length(xJ)-1; %norma 2
semilogy(1:NJ1,xJn);hold on,semilogy(1:4:NJ1,xJn(1:4:NJ1),'o')
xGSn=sqrt(sum((diff(xGS'))'.^2)); NGS1=length(xGS)-1; grid
semilogy(1:NGS1,xGSn);semilogy(1:4:NGS1,xGSn(1:4:NGS1),'+')
xJwn=sqrt(sum((diff(xJw'))'.^2)); NJw1=length(xJw)-1;
semilogy(1:NJw1,xJwn);semilogy(1:4:NJw1,xJwn(1:4:NJw1),'s')
xGSwn=sqrt(sum((diff(xGSw'))'.^2)); NGSw1=length(xGSw)-1;
semilogy(1:NGSw1,xGSwn);semilogy(1:4:NGSw1,xGSwn(1:4:NGSw1),'x')
Ir al inicio
Estudio
de la rapidez de convergencia:
Se presupone
que se han ejecutado todos los métodos iterativos anteriores.
Rapidez
(orden) de convergencia para el método de Jacobi sin y con
relajación:
mJ=mean(diff(log10(xJn))),mJ=10^mJ %Alternativa por Pitagoras
rho_mJ=max(abs(eig(MJ)))
mJw=polyfit(1:NJw1-1,log10(xJwn(1:NJw1-1)),1),mJw=10^mJw(1)
w=1.0683;rho_MJw=max(abs(eig(eye(3)-w*D\A)))
Rapidez (orden) de convergencia
para el método de Gauss-Seidel sin y con sobrerelajación:
mGS=mean(diff(log10(xGSn(1:NGS1-1)))),mGS=10^mGS %Alternativa
por Pitagoras
rho_mGS=max(abs(eig(MGS)))
mGSw=polyfit(1:NGSw1-1,log10(xGSwn(1:NGSw1-1)),1),mGSw=10^mGSw(1)
w=1.0035;rho_MGSw=max(abs(eig((D+w*L)\((1-w)*D-w*U))))
Ir al inicio