Problem dividing matrix by vector

Hi

I have a matrix of 20301*20301

and a vector of 20301*1 ...(after 201 iterations, so at cell 202 the numbers start.Before it was all zeros)

When I try dividing the using A\b, the error "Warning:Singular Matrix .." comes up and my answer is just NaNs. I compared my answer to a friend and while he got the exact values of A and b his division worked and he got actual values. I was wondering if anyone could shed light on this situation

[the cyclist]

Can you upload your data (in a mat file) and the code you used (in either an m file, or by pasting the code)? Please don't paste an screenshot of the code, which would require us to retype it out to use it.

[the cyclist]

Are you sure that your friend has identical inputs? It seems to me that between you and your friend, you must have at least one of the following, if you are getting different results:

• different inputs
• different calculation
• different versions of MATLAB
• different hardware

(or maybe some other difference I am forgetting).

[Torsten]

Do you really want to solve A*x = b for x or do you only want to divide each row (column) of A by the vector b ?

[Maaz Madha]

clear
clc 
clear all

%%code for temperature distribution over a sheet of metal using discertisation..
%uses central differencing method 

%%Initialisation 1 
alpha=0.1;
beta=0.1;
dx=(2*pi)/100;%%increments
dy=dx;
dt=0.01; %%timeskip
nu=0.1;
H=2*pi;%%height of metal sheet
L=4*pi;%length
T_end=5;
x=[0:dx:L];%%setting up x and y values
y=[0:dx:H];



n=round(size(x,2));
m=round(size(y,2));
t=round((T_end/dt)+1);
v=@(x,t) 5*(1-((x-2*pi)^2)/(4*pi^2))*cos(pi*t)*sin(x);%%velocity distribution
T=@(x,y) 20*cos(x)*sin(y);%%initial temperature
A=zeros(n*m);%%preallocation
b=zeros(n*m,1);


%A matrix
for i=1:n
    for k=1:t
        an(i,k)= v((i-1)*dx,dt)*dt/(2*dy) - beta*dt/dy^2;%%discetising (a_north)
        as(i,k)= -v((i-1)*dx,dt)*dt/(2*dy) - beta*dt/dy^2; %a south
    end 
end 

ae = -dt*alpha/dx^2; %a east
ap = 1+(2*alpha*dt/dx^2)+(2*beta*dt/dy^2);%a present
aw=ae; %a west


%T1
for i=2:n-1
    for j=2:m-1
        pointer=(j-1)*n+i;%to map 2d onto 1d
        A(pointer,pointer)=ap;
        A(pointer,pointer-n)=as(i);
        A(pointer,pointer-1)=aw;
        A(pointer,pointer+1)=ae;
        A(pointer,pointer+n)=an(i); 
        b(pointer)=T((i-1)*dx,(j-1)*dy);
      
    end 
     
end 

%T2
 i=1;
    for j=2:m-1
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,pointer+n)=an(i);
        A(pointer,pointer-n)=as(i);
        A(pointer,pointer+i)=aw+ae;
        b(pointer)=T((i-1)*dx,(j-1)*dy);
        
       
        
    end 


%T3

for i=n
    for j=2:m-1
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,pointer-n)=as(i);
        A(pointer,pointer-1)=aw+ae;
        A(pointer,pointer+n)=an(i);
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 

%T4 

for i=2:n-1
    for j=1
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,pointer+n)=an(i);
        A(pointer,pointer-1)=aw;
        A(pointer,pointer+1)=ae;
        A(pointer,(m-2)*n+n)=as(i);
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 
%for T5
for i=2:n-1
    for j=m-1
        pointer=(j-1)*n+i;    
        A(pointer,pointer)=ap;
        A(pointer,pointer+n)=as(i);
        A(pointer,pointer-1)=aw;
        A(pointer,pointer+n)=an(i);
        A(pointer,pointer+1)=ae;
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 

%T6 

for i=1
    for j=1
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,pointer+n)=an(i);
        A(pointer,(m-2)*n+n)=as(i);
        A(pointer,pointer+1)=aw+ae;   
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 
        
%T7

for i=n
    for j=1
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,(m-2)*n+n)=as(i);
        A(pointer,pointer-1)=(aw+ae);
        A(pointer,pointer+n)=an(i);
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 


%T8
for i=1
    for j=m
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,pointer-1)=as(i);
        A(pointer,pointer+1)=(aw+ae);
        A(pointer,n+i)=an(i);
        b(pointer)=T((i-1)*dx,(j-1)*dy);
    end 
end 

%t9

for i=n
    for j=m
        pointer=(j-1)*n+i;
        A(pointer,pointer)=ap;
        A(pointer,(m-2)*n+n)=as(i);
        A(pointer,pointer-1)=aw+ae;
        A(pointer,i-n+2)=an(i);
        b(pointer)=T((i-1)*dx,(j-1)*dy);
        
    end
   
end 

Temp=A\b;

[Maaz Madha]

I want to solve for temperature in the format Ax=b

[Torsten]

The code will neither produce a matrix A nor a vector b.

[Maaz Madha]

it does produce a matrix A and a vector b when I run it

[Torsten]

an and as are matrices of two dimensions - you address them as arrays of one dimension.

T((i-1)*dx,(j-1)*dy) is usually undefined since (i-1)*dx and (j-1)*dy are not integers in general.

[Maaz Madha]

hmmm i see your point about as and an. About the b vector my purpose was to input data in using the initial temperature at k=1 (time zero) which was

T=@(x,y) 20*cos(x)*sin(y);%%initial temperature

So my thought was that by doing (i-1)*dx,(j-1)*dy over the whole grid it would give me all the initial temperatures across those points. What would you recommend

[Torsten]

Sorry, I did not see that T was not a matrix, but a function.

[Maaz Madha]

i've fixed as and an values

% Discretisation
for k=2:t
    for i=1:n       
        v(i)=5*(1-((i-1)*dx-2*pi)^2/(4*pi^2))*cos(pi*(k-1)*dt)*sin((i-1)*dx); 
    end 
    as=-v*dt/(2*dy) - beta*dt/dy^2; % Calculate the coefficient matix of Tn(i,j-1)
    an=v*dt/(2*dy) - beta*dt/dy^2;  % Calculate the coefficient matix of Tn(i,j+1)    
    ae = -dt*alpha/dx^2; %discretised value of 'ae'
    ap = 1+(2*alpha*dt/dx^2)+(2*beta*dt/dy^2); %discretised value of 'ap'
    aw=ae; %discretised value of 'aw' which is the same as 'ae'
end 

but despite this the answer produced is still sinfular

[Torsten]

Equations for T2 - T9 must be derived from the boundary conditions of your problem to get a nonsingular matrix A. I don't see where you incorporate boundary conditions.

Further, beginning with time step 2, T in the equations must be replaced by the Temp vector of the preceeding time step. Thus it will be better to work with Temp instead of T in the definition of b.

[Maaz Madha]

i used the boundary conditions when i derived the eqns for T2-T9. Using the boundary conditions I had to use the pointer function to manipulate the neighbours of ap(the basis of everything).

I am working with temperature as temperature is T is my code. If you look it says T=@(x,y) 20*cos(x)*sin(y) and this is what I have been using.

[Torsten]

Which boundary conditions did you try to incorporate at the four edges ?

T in your code is the initial temperature at time t=0. If you want to stick to T in your b vector, you will have to redefine your function for T since T becomes Temp.

[Maaz Madha]

In the Q, the boundary conditions are and T is periodic

[Maaz Madha]

At the 4 edges basically I had to look at periodicity so and