Problem with Singular Matrices
显示 更早的评论
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)
回答（1 个）
Walter Roberson
2016-2-11
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 个评论
Jamil Villarreal
2016-2-11
Walter Roberson
2016-2-11
I am not surprised; if I remember my algebra correctly, matrix multiplication cannot increase the rank.
Jamil Villarreal
2016-2-12
类别
在 帮助中心 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
2016-2-11
2016-2-12
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
