Problem dividing matrix by vector

2 次查看（过去 30 天）

显示更早的评论

Maaz Madha 2020-3-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513561-problem-dividing-matrix-by-vector

评论： Maaz Madha 2020-3-28

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

16 个评论
显示 14更早的评论隐藏 14更早的评论

Maaz Madha 2020-3-28

编辑：Maaz Madha 2020-3-28

在 MATLAB Online 中打开

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;