Error: Warning: Matrix is singular to working precision

Hello...I am currently working on a project for university and I need to create an A and b matrix to solve for Ax=b. When I run my code I get the error: "Warning: Matrix is singular to working precision." which I know means my A matrix is non-invertible, correct? would someone be able to please help me understand how to fix this?

Thank you!

I believe all the code you'll need is below...please let me know if you need more. I've been testing using an n = 20 value.

CODE:

%Create Mesh based on n
meshsize = n*n/2;
nodes = 1:meshsize;

%Make Matricies
nodetype = 16*ones(n/2, n);
A = zeros(n*(n/2),n*(n/2));
b = zeros(n*(n/2),1);

%Make Meshgrid
d_x = 2/n;
[X,Y] = meshgrid(d_x/2:d_x:2,1-d_x/2:-d_x:d_x/2); 

%Problem Givens
T1 = 28;                %C
T2 = 50;                %C
T3 = @(y) 10 + 14*y;    %C
T_inf = 32;             %C
h = 28;                 %W/m^2*K
k = 12;                 %W/m*K
q_dot = 47000;          %W/m^3

nodetype = nodetype';
nodetype_vector = reshape(nodetype,n*(n/2),1);

for i = 1:length(nodetype_vector)
    if nodetype(i) == 0
        A(i,i) = 0;
        b(i) = 0;
    elseif nodetype(i) == 1
        %Make A matrix
        A(i,i) = -6;
        A(i,i+1) = 1;
        A(i,i+n) = 1;
        %Make B Matrix
        b(i) = -2 * (T1 + T2);
    elseif nodetype(i) == 2
        %Make A matrix
        A(i,i) = -6;
        A(i,i-1) = 1;
        A(i,i+n) = 1;
        %Make B Matrix
        b(i) = -2 * (T2 + T3(Y(i)));
    elseif nodetype(i) == 3
        %Make A Matrix
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,i+n) = 1;
        %Make B Matrix
        b(i) = -2 * T1;
    elseif nodetype(i) == 4
        %Make A Matrix
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,in) = 1;
        %Make B Matrix
        b(i) = 2 * T3(Y(i));
    elseif nodetype(i) == 5
        %Make A Matrix
        A(i,i) = -5;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i+n) = 1;       
        %Make B Matrix
        b(i) = 2 * T2;
    elseif nodetype(i) == 6
        %Make A Matrix
        A(i,i) = -5;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = 2 * T1;
    elseif nodetype(i) == 7
        %Make A Matrix
        A(i,i) = -5;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = 2 * T3(Y(i));
    elseif nodetype(i) == 8
        %Make A Matrix
        A(i,i) = -3;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = 0;
    elseif nodetype(i) == 9    
        %Make A Matrix
        A(i,i) = -(3 + ((h*Y(i))/2));
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = -((h*Y(i))/2)*T_inf;
    elseif nodetype(i) == 10    
        %Make A Matrix
        A(i,i) = -(3 + ((h*Y(i))/2));
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = -((h*Y(i))/2)*T_inf;
    elseif nodetype(i) == 11    
        %Make A Matrix
        A(i,i) = -(2 + ((h*Y(i)))/2);
        A(i,i-1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = -((h*Y(i))/2)*T_inf;
    elseif nodetype(i) == 12
        %Make A Matrix
        A(i,i) = -(2 + ((h*Y(i)))/2);
        A(i,i-1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = -((h*Y(i))/2)*T_inf;
    elseif nodetype(i) == 13
        A(i,i) = -(3 + ((h*X(i))/2));
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i-n) = 1;
        %Make B Matrix
        b(i) = -((h*X(i))/2)*T_inf;
    elseif nodetype(i) == 14
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i+n) = 1;       
        A(i,i-n) = 1;   
        %Make B Matrix
        b(i) = - ((q_dot * X(i)^(2))/(2*k));
    elseif nodetype(i) == 15
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i+n) = 1;       
        A(i,i-n) = 1;   
        %Make B Matrix
        b(i) = - ((q_dot * X(i)^(2))/(k));
    elseif nodetype(i) == 16  
        %Make A Matrix
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i+n) = 1;       
        A(i,i-n) = 1;       
        %Make B Matrix
        b(i) = 0;
    end
end

T = A\b;

Answer 1

[Walter Roberson]

You have

nodetype_vector = reshape(nodetype,n*(n/2),1);

so nodetype_vector is certain to be length n*n/2 .

Then you have

for i = 1:length(nodetype_vector)

so that is for i = 1 : n*n/2 . In the case of n=20 that would be for i = 1 : 200

With the way you initialize, all of your nodetype_vector entries come out equal to 16. So let us look at the code for that:

    elseif nodetype(i) == 16  
        %Make A Matrix
        A(i,i) = -4;
        A(i,i-1) = 1;
        A(i,i+1) = 1;
        A(i,i+n) = 1;       
        A(i,i-n) = 1;       
        %Make B Matrix
        b(i) = 0;

On the first iteration, i is 1. A(1,1) = -4 -- okay. A(1,1-1) = 1 is A(1,0)=1 and we have a problem. If we got past that then A(i,i-n) would be A(1,1-20)=1 which would be A(1,-19)=1 which is a big problem.

You can alter the for to

for i = n+1:length(nodetype_vector)-n

but when you do, you never store anything in columns 1:n or end-n+1:end. Whenever you have columns that are all-identical that reduces the rank of the matrix by 1, so you would end up with an A matrix which is rank at most n*n/2 - 2*n = n*(n-4)/2 . That will never be intervertable.

[Josh Glenn]

so would it just be my A(i,i-1) that is "broken"? Should I change it to something different such as A(i,i+2)? Would that make sense?

Thank you!

[Walter Roberson]

A(i,i-n) is a problem as well.

Potential revised version enclosed. Please pay special attention to the code for nodetype 4, which had a typing mistake in your original.

With your initialization of everything to nodetype 16, you initialize b to all 0, and the result of the \ operator turns out all 0...

[Josh Glenn]

Okay...I understand what you're saying here. However, for nodes such as node 16 where the A matrix has 5 inputs, it must have those 5 inputs so that the x vector (T) can be solved for correctly so I'm not sure an if statement would work in this sense. Does that make sense?

[Walter Roberson]

For any one given instance, will the nodetypes all be the same, or could there be a mix of node types?

If the node types are all the same, then your code is construction of a band matrix, with the exact number of bands depending on the node type. The question then becomes "what if you reach the edge and there is not room for a band of full size?" -- which is a situation that is certain to happen for every node type except type 0. If it makes sense mathematically to just let the band end at the edge, then the code I supplied implements that. If it would instead make sense to "wrap", then a bit more work would be required. "Wrapping" would be like where when an edge is hit then it continues in the same direction starting from the other side

    0 0 0 0 0 1 0 0
    0 0 0 0 0 0 1 0
    0 0 0 0 0 0 0 1
    1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0
    0 0 1 0 0 0 0 0

If all node types are the same, then a more efficient implementation would by using diag() -- but it is probably worth getting the implementation correct first before worrying about the efficiency in building it.

[Josh Glenn]

Here's my full code...It is essentially creating a heat map due to heat transfer conditions and plotting the T matrix which (as of right now, in non-invertible). Does the whole file help?

[Walter Roberson]

Use this code, but put in a breakpoint at the last executable line T = A\b;

Invoke with

newProject(20)

When the code stops executing, command

figure(); spy(A)

Notice the horizontal gap at about 167. Now command

rank(A)

Notice it says 190, not 200.

So your initialization of node types is leaving a gap, probably because of the lines

    nodetype(0.8*(n/2):(n/2),0.3*n:0.6*n) = 0;
    %EMPTY SPACES

That gap is causing the A matrix to be singular. All-zero rows or all-zero columns cause singularity.