- #1
mfarooq52003
- 22
- 0
Hi There
Can anyone check the following algorithm and tell me what I am doing wrong. If every thing is OK then please give me some feedback to make it robust. Please check both the algorithms
Algorithm-1
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
r(:,1)=b(:,1)-A*x(:,1);
y(:,1)=y;
p=A*r(:,1);
p1=A*p;
c0=y'*r(:,1);
c1=y'*p;
c2=y'*p1;
c3=y'*A*p1;
delta=c1*c3-c2^2;
alpha=(c0*c3-c1*c2)/delta;
beta=(c0*c2-c1^2)/delta;
r(:,2)=r(:,1)-(c0*p)/c1;
x(:,2)=x(:,1)+(c0/c1)*r(:,1);
r(:,3)=r(:,1)-alpha*p+beta*p1;
x(:,3)=x(:,1)+alpha*r(:,1)-beta*p;
y(:,2)=A'*y(:,1);
y(:,3)=A'*y(:,2);
y(:,4)=A'*y(:,3);
eps=0.0001;
tic;
%Step-2
for k=3:m
y(:,k+2)=A'*y(:,k+1);
q1=A*r(:,k-1);
q2=A*q1;% q2=A*q1=A*(A*r(:,k-1))=A^2r(:,k-1)
q3=A*r(:,k-2);
a11=y(:,k-1)'*r(:,k-1);
a13=y(:,k-2)'*r(:,k-2);
a21=y(:,k)'*r(:,k-1);
a22=a11;
a23=y(:,k-1)'*r(:,k-2);
a31=y(:,k+1)'*r(:,k-1);
a32=a21;
a33=y(:,k)'*r(:,k-2);
s=y(:,k+2)'*r(:,k-1);
t=y(:,k+1)'*r(:,k-2);
F_{k+1}=-a11/a13;
b1=-a21-F_{k+1}*a23;
b2=-a31-F_{k+1}*a33;
b3=-s-F_{k+1}*t;
Delta_{k+1}=a11*(a22*a33-a32*a23)+a13*(a21*a32-a31*a22);
B_{k+1}=1/Delta_{k+1}*[b1*(a22*a33-a32*a23)+a13*(b2*a32-b3*a22)];
G_{k+1}=[b1-B_{k+1}*a11]/a13;
C_{k+1}=[b2-B_{k+1}*a21-G_{k+1}*a23]/a22;
A_{k+1}=1/[C_{k+1}+G_{k+1}];
r(:,k+1)=A_{k+1}*[q2+B_{k+1}*q1+C_{k+1}*r(:,k-1)+F_{k+1}*q3+G_{k+1}*r(:,k-2)];
x(:,k+1)=A_{k+1}*[C_{k+1}*x(:,k-1)+G_{k+1}*x(:,k-2)-(q1+B_{k+1}*r(:,k-1)+F_{k+1}*r(:,k-2))];
%Step-3
if (r(:,k+1))~=eps
continue;
else
x=x(:,k+1);
end
end
toc;
x1=A\b;
x2=x(:,k+1);
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)
Algorithm-2
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
y(:,1)=y;
r(:,1)=b-A*x(:,1);
p(:,1)=r(:,1);
C=1;
A_{2}=(-1*y(:,1)'*r(:,1))/(y(:,1)'*A*r(:,1));
x(:,2)=x(:,1)-A_{2}*r(:,1);
r(:,2)=r(:,1)+A_{2}*A*r(:,1);
eps=0.0001;
tic;
%Step-2
for k=2:m
y(:,k)=A'*y(:,k-1);
D_{k+1}=(-1*y(:,k)'*r(:,k))/(C*y(:,k)'*p(:,k-1));
p(:,k)=r(:,k)+D_{k+1}*C*p(:,k-1);
A_{k+1}=(-1*y(:,k)'*r(:,k))/(y(:,k)'*A*p(:,k));
r(:,k+1)=r(:,k)+A_{k+1}*A*p(:,k);
x(:,k+1)=x(:,k)-A_{k+1}*p(:,k);
%Step-3
if r(:,k+1)~=eps,
C=C/A_{k};
fprintf('I am here')
else
x=x(:,k+1);
fprintf('You are here')
end
end
toc;
x2=x(:,k+1);
x1=A\b;
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)
Please let me know of any blunder if I have done there in the loops and I will appreciate if you can correct it and post it in correct form.
I will be grateful for your help and support.
Regards
Can anyone check the following algorithm and tell me what I am doing wrong. If every thing is OK then please give me some feedback to make it robust. Please check both the algorithms
Algorithm-1
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
r(:,1)=b(:,1)-A*x(:,1);
y(:,1)=y;
p=A*r(:,1);
p1=A*p;
c0=y'*r(:,1);
c1=y'*p;
c2=y'*p1;
c3=y'*A*p1;
delta=c1*c3-c2^2;
alpha=(c0*c3-c1*c2)/delta;
beta=(c0*c2-c1^2)/delta;
r(:,2)=r(:,1)-(c0*p)/c1;
x(:,2)=x(:,1)+(c0/c1)*r(:,1);
r(:,3)=r(:,1)-alpha*p+beta*p1;
x(:,3)=x(:,1)+alpha*r(:,1)-beta*p;
y(:,2)=A'*y(:,1);
y(:,3)=A'*y(:,2);
y(:,4)=A'*y(:,3);
eps=0.0001;
tic;
%Step-2
for k=3:m
y(:,k+2)=A'*y(:,k+1);
q1=A*r(:,k-1);
q2=A*q1;% q2=A*q1=A*(A*r(:,k-1))=A^2r(:,k-1)
q3=A*r(:,k-2);
a11=y(:,k-1)'*r(:,k-1);
a13=y(:,k-2)'*r(:,k-2);
a21=y(:,k)'*r(:,k-1);
a22=a11;
a23=y(:,k-1)'*r(:,k-2);
a31=y(:,k+1)'*r(:,k-1);
a32=a21;
a33=y(:,k)'*r(:,k-2);
s=y(:,k+2)'*r(:,k-1);
t=y(:,k+1)'*r(:,k-2);
F_{k+1}=-a11/a13;
b1=-a21-F_{k+1}*a23;
b2=-a31-F_{k+1}*a33;
b3=-s-F_{k+1}*t;
Delta_{k+1}=a11*(a22*a33-a32*a23)+a13*(a21*a32-a31*a22);
B_{k+1}=1/Delta_{k+1}*[b1*(a22*a33-a32*a23)+a13*(b2*a32-b3*a22)];
G_{k+1}=[b1-B_{k+1}*a11]/a13;
C_{k+1}=[b2-B_{k+1}*a21-G_{k+1}*a23]/a22;
A_{k+1}=1/[C_{k+1}+G_{k+1}];
r(:,k+1)=A_{k+1}*[q2+B_{k+1}*q1+C_{k+1}*r(:,k-1)+F_{k+1}*q3+G_{k+1}*r(:,k-2)];
x(:,k+1)=A_{k+1}*[C_{k+1}*x(:,k-1)+G_{k+1}*x(:,k-2)-(q1+B_{k+1}*r(:,k-1)+F_{k+1}*r(:,k-2))];
%Step-3
if (r(:,k+1))~=eps
continue;
else
x=x(:,k+1);
end
end
toc;
x1=A\b;
x2=x(:,k+1);
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)
Algorithm-2
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
y(:,1)=y;
r(:,1)=b-A*x(:,1);
p(:,1)=r(:,1);
C=1;
A_{2}=(-1*y(:,1)'*r(:,1))/(y(:,1)'*A*r(:,1));
x(:,2)=x(:,1)-A_{2}*r(:,1);
r(:,2)=r(:,1)+A_{2}*A*r(:,1);
eps=0.0001;
tic;
%Step-2
for k=2:m
y(:,k)=A'*y(:,k-1);
D_{k+1}=(-1*y(:,k)'*r(:,k))/(C*y(:,k)'*p(:,k-1));
p(:,k)=r(:,k)+D_{k+1}*C*p(:,k-1);
A_{k+1}=(-1*y(:,k)'*r(:,k))/(y(:,k)'*A*p(:,k));
r(:,k+1)=r(:,k)+A_{k+1}*A*p(:,k);
x(:,k+1)=x(:,k)-A_{k+1}*p(:,k);
%Step-3
if r(:,k+1)~=eps,
C=C/A_{k};
fprintf('I am here')
else
x=x(:,k+1);
fprintf('You are here')
end
end
toc;
x2=x(:,k+1);
x1=A\b;
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)
Please let me know of any blunder if I have done there in the loops and I will appreciate if you can correct it and post it in correct form.
I will be grateful for your help and support.
Regards