Problem with Singular Matrices
3 views (last 30 days)
Show older comments
I am working on a crank-nicolson scheme for solving the Allen-Cahn equation. Below is the code that I have written, but I keep on getting the warning of matrix is singular to working precision. I have also attached the paper scheme on which the code is based on for reference.
Many thanks in advance, any help on getting this program running will be greatly appreciated.
Here is my code:
clear all
epsilon=10^(-3);
xL=-5;xR=-xL;
yD=-5;yU=-yD;
L=xR-xL; %the length of domain
tau=0.02; %time step "k"
h=0.2; %spatial step
%x=(xL+h):h:(xR-h); y=(yD+h):h:(yU-h);
x=xL:h:xR; y=yD:h:yU;
r=tau/h^2;
N=L/h; T=10;
Utemp=zeros(N+1,N+1);
x0=floor((N+1)/2);y0=x0;
for i=1:N+1
for j=1:N+1
if ((i-x0)^2+(j-y0)^2)<=5^2
Utemp(i,j)=1;
end
end
end
%mesh(Utemp)
U=reshape(Utemp',(N+1)^2,1);
%U=sparse(U);
%mesh(U)
%B1=diag((2+4*r)*ones(N-1,1))+diag(-r*ones(N-2,1),1)+diag(-r*ones(N-2,1),-1);
%B2=diag((2-4*r)*ones(N-1,1))+diag( r*ones(N-2,1),1)+diag( r*ones(N-2,1),-1);
K1 = diag([0,-r*ones(1,N-1),0]);
K11 = diag([0,r*ones(1,N-1),0]);
for t=0:tau:T
% for the left parts
A = kron(diag([ones(1,N-1),0],-1),K1) + kron(diag([0,ones(1,N-1)],1),K1);
% for the right parts
A1 = kron(diag([ones(1,N-1),0],-1),K11) + kron(diag([0,ones(1,N-1)],1),K11);
for i=2:N
Ui = Utemp(i,:);
K2 = diag([1,(2+4*r-tau*epsilon^(-2)*(1-Ui(2:N).^2)),1]) + diag([-r*ones(1,N-1),0],-1) + diag([0,-r*ones(1,N-1)],1);
A = A + kron(diag([zeros(1,i-1),1,zeros(1,N+1-i)]),K2);
K22 = diag([1,(2-4*r+tau*epsilon^(-2)*(1-Ui(2:N).^2)),1]) + diag([r*ones(1,N-1),0],-1) + diag([0,r*ones(1,N-1)],1);
A1 = A1 + kron(diag([zeros(1,i-1),1,zeros(1,N+1-i)]),K22);
end
A = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
A1 = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
U=(A\A1)*U;
Utemp=reshape(U,N+1,N+1);
Utemp=Utemp';
pause
mesh(Utemp)
end
%U=reshape(Utemp,N+1,N+1);
mesh(Utemp)
2 Comments
Answers (1)
Walter Roberson
on 11 Feb 2016
You have
A = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
A1 = kron(diag([1,zeros(1,N-1),1]),eye(N+1));
U=A\A1*U;
Your A and A1 are the same and are rank 2*(N+1) but size (N+1)^2 x (N+1)^2. You then \ those identical low-rank matrices. That is the cause of the singularity warning.
Remember that A\A1*U is (A\A1)*U not A\(A1*U)
If somehow A and A1 were not singular, then because they are identical, A\A1 would be the numeric approximation of the identity matrix.
3 Comments
Walter Roberson
on 11 Feb 2016
I am not surprised; if I remember my algebra correctly, matrix multiplication cannot increase the rank.
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!