Warning: Matrix is singular to working precision. Help

So when the code gets to the line "uhat=A\b" it says Warning: Matrix is singular to working precision and every entry is NaN. I have a feeling it is because b is mostly zeros, which is probably wrong and I'll be asking my professor. I am using m=8. What am I doing wrong?

function finvolHW(m)

 A=zeros(2*m^2);
A(1,1)=-3;
A(1,2)=1;
A(1,9)=1;
for i=2:(2*m)-1
    A(i,i-1)=1;
    A(i,i)=-4;
    A(i,i+1)=1;
    A(i,i+8)=1;
end
for i=2*m
    A(i,i-1)=1;
    A(i,i)=-3;
    A(i,i+8)=1;
end
for i=1:m-2
    b=((i*(2*m))+1);
    A(b,b-8)=1;
    A(b,b)=-3;
    A(b,b+1)=1;
    A(b,b+8)=1;
end
for i=0:m-3
    for j=((m*2)+2)+(2*i*m):1:((4*m)-1)+((2*i*m))
        A(j,j-8)=1;
        A(j,j-1)=1;
        A(j,j)=-4;
        A(j,j+1)=1;
        A(j,j+8)=1;
    end
end
for i=2:1:m-2
    for j=(m*i*2)
        A(j,j-8)=1;
        A(j,j-1)=1;
        A(j,j)=-3;
        A(j,j+8)=1;      
    end
end
for i=2*(m^2)-(2*m)+1
    A(i,i-8)=1;
    A(i,i)=-3;
    A(i,i+1)=1;
end
for i=2*(m^2)-(2*m)+2:1:2*(m^2)-(2*m)+(m)
    A(i,i-8)=1;
    A(i,i-1)=1;
    A(i,i)=-4;
    A(i,i+1)=1;   
end
for i=(2*m^2)-m+1:1:2*m^2-1
    A(i,i-1)=1;
    A(i,i)=-3;
    A(i,i+1)=1;
end
for i=2*m^2
    A(i,i-1)=1;
    A(i,i)=-2;
end
%A has been created

 b=zeros((2*m^2),1);
for i=1:1:2*m
    b(i,1)=-1;
end
uhat=A\b;

Answer 1

[Jan]

The created matrix A is singular or almost singular. "Singular" means, that the matrix does not have an inverse, like e.g. the scalar 0 does not have an inverse value also. Therefore either the problem cannot be solved or you create the wrong matrix. Because we see only the code you use to create A, there is no chance that we could guess, if there is anything wrong.

Lines like this:

for j=(m*i*2)

are confusing only. This is a loop over 1 element, such that it would be smarter to write:

j = m*i*2;

[Kyle]

Here is the problem statement I have

-Write a MATLAB function that solves the problem of two-dimensional heat ow in a rectangular bar with no heat sources or sinks. The bar is two units long and one unit tall. The entire top boundary of the bar is xed at a normalized temperature of u = 1. Exactly one half of the bottom boundary is fixed at u = 0. The remainder of the boundary is insulated,with no heat flow in or out.Use finite difference, finite volume or finite elements to solve the problem. The function should take the total number of intervals m along one vertical column as input. All intervals will be square, meaning that the total number of discrete squares on the domain will be 2m^2. Write a script to call the function for m = 4, 8 and 16. The script will also plot the solution, by putting it into a rectangular matrix and using the MATLAB function contourf().

[Jan]

What kind of help do you expect now? I cannot see the correlation between parts of the assignment text and you code, because the code is free of comments. Although I could rewrite the solution from scratch and compare it with your code until I find the difference. But this would be your job, actually.

[Serap Karahanogullari]

Hello I have the same warning and i cannot find the solution why its so. I would be appereciated if you could help me. here you can find my all code with all comments.

many thank in advance

{%********************************************************* % The namespaces have been normalized. The following % table shows the attribuation. % Normalized Name --> Full Name ---> User-defined Name % =================================== % e0 --> e[0]97618 --> %*********************************************************

%********************************************************* % The variables are named according to the notation % provided in the Mosaic model. % % The variable names can be read as follows: % ========================================== % e0_greek_alpha_aus % α: % Superscripts % aus: null % % e0_greek_alpha_in % α: % Superscripts % in: null % % e0_ep % ep: % % e0_l % l: % % e0_lambda % lambda: % % e0_ro % ro: % % e0_sig % sig: % % e0_T_gas % T: % Superscripts % gas: null % % e0_T_wand % T: % Superscripts % wand: null % % e0_cp_luft % cp: % Superscripts % luft: null % % e0_d % d: % % e0_flaeche % flaeche: % % e0_Q_kon % Q: % Superscripts % kon: null % % e0_Q_st % Q: % Superscripts % st: null % % e0_U % U: % % e0_deltat % deltat: % %*********************************************************

function[X,Y]=solveEquationSystem()

% declare interval of independent variable X_START=0.0; X_END=1000; X_INTERVAL=[X_START X_END]; % e0_t

% specify the initial values for the state variables Y_INIT(1) = 100.0; % e0_deltat Y_INIT(2) = 100.0; % e0_Q_kon Y_INIT(3) = 50.0; % e0_Q_st Y_INIT(4) = 10.0; % e0_U Y_INIT(5) = 5.0; % e0_flaeche

% declare parameters PARAMS(1) = 23.0; % e0_greek_alpha_aus PARAMS(2) = 8.1; % e0_greek_alpha_in PARAMS(3) = 395.0; % e0_T_gas PARAMS(4) = 300.0; % e0_T_wand PARAMS(5) = 1005.0; % e0_cp_luft PARAMS(6) = 0.5; % e0_d PARAMS(7) = 0.04; % e0_ep PARAMS(8) = 0.5; % e0_l PARAMS(9) = 150.0; % e0_lambda PARAMS(10) = 2.7; % e0_ro PARAMS(11) = 0.0; % e0_sig

M = daeSystemMM();

OPTIONS = odeset('Mass',M);

[X,Y] = ode45(@(X,Y)daeSystemLHS(X,Y,PARAMS),X_INTERVAL,Y_INIT',OPTIONS);

displayResults(X,Y);

end

function MASS = daeSystemMM()

MASS = zeros(5, 5); %total number of equations MASS(1:1, 1:1)=eye(1,1); % number of odes

end

% evaluate the differential function. function[DYDX] = daeSystemLHS(X,Y,PARAMS)

% read out variables e0_deltat = Y(1); e0_Q_kon = Y(2); e0_Q_st = Y(3); e0_U = Y(4); e0_flaeche = Y(5);

% read out differential variable e0_t = X;

% read out parameters e0_greek_alpha_aus = PARAMS(1); e0_greek_alpha_in = PARAMS(2); e0_T_gas = PARAMS(3); e0_T_wand = PARAMS(4); e0_cp_luft = PARAMS(5); e0_d = PARAMS(6); e0_ep = PARAMS(7); e0_l = PARAMS(8); e0_lambda = PARAMS(9); e0_ro = PARAMS(10); e0_sig = PARAMS(11);

% evaluate the function values DYDX(1) = 4* (e0_Q_kon + e0_Q_st )/(e0_ro * e0_l * e0_cp_luft * ( e0_d )^2); DYDX(2) = e0_deltat - ( e0_T_gas - e0_T_wand ); DYDX(3) = e0_U - ( ( ( 1/( e0_greek_alpha_in ) ) + ( 1/( e0_greek_alpha_aus ) ) + ( e0_d/e0_lambda ) )^- 1 ); DYDX(4) = e0_Q_kon - e0_U * e0_flaeche * e0_deltat; DYDX(5) = e0_Q_st - (e0_ep * e0_sig * e0_flaeche *( e0_T_wand )^4);

DYDX=DYDX';

end

function[] = displayResults(X,Y)

% decide for a plot type: % 0 -> Plot the variables into individual figures % 1 -> Plot into sub figures % 2 -> Plot all selected into one figure % other -> Do not plot plotType = 2;

% set a line width: linewidth = 1.5;

% define which dependent variables should be plotted % 1 -> Plot. % other -> Do not plot. plotControl=[ 1 % e0_deltat 1 % e0_Q_kon 1 % e0_Q_st 1 % e0_U 1 % e0_flaeche ];

% decide wether to normalize the y axis % 1 -> Normalized % other -> Individual maximum scale axisControl = 2;

%====================================================

% labels of the dependent variables yAxisLabels=[ 'e0.deltat ' % e0_deltat 'e0.Q^{kon}' % e0_Q_kon 'e0.Q^{st} ' % e0_Q_st 'e0.U ' % e0_U 'e0.flaeche' % e0_flaeche ]; xAxisLabel = 't';

% plot the variables figureIndex=1; xMinVal = min(X); xMaxVal = max(X); yMinVal = min(min(Y)); yMaxVal = max(max(Y)); if (plotType==0) % create a plot for each state variable individually for i=1:length(Y(1,:)) if (plotControl(i)==1) figure(i) plot(X,Y(:,i),'LineWidth',linewidth) title('Solution of the equation system'); xlabel(xAxisLabel); ylabel(yAxisLabels(i,:)); if (axisControl==1) axis([xMinVal xMaxVal yMinVal yMaxVal]); end legend(yAxisLabels(i,:)); figureIndex=figureIndex+1; end end elseif (plotType==1) % use a subplot environment firstTime = 1; numberOfPlots = sum(plotControl); figure(1) for i=1:length(Y(1,:)) if (plotControl(i)==1) subplot(numberOfPlots,1,figureIndex); plot(X,Y(:,i),'LineWidth',linewidth) ylabel(yAxisLabels(i,:)); legend(yAxisLabels(i,:)); if (axisControl==1) axis([xMinVal xMaxVal yMinVal yMaxVal]); end figureIndex=figureIndex+1; if (firstTime) title('Solution of the equation system'); firstTime = 0; end end end xlabel(xAxisLabel); elseif (plotType==2) % plot in one figure colors = [ 'r' % -> Red 'g' % -> Green 'b' % -> Blue 'c' % -> Cyan 'm' % -> Magenta 'y' % -> Yellow 'k' % -> Black ]; colorCtr=1; maxColors=7; figure(1) hold on; for i=1:length(Y(1,:)) if (plotControl(i)==1) linespec = colors(colorCtr); if (colorCtr<=maxColors) colorCtr = colorCtr+1; end plot(X,Y(:,i),linespec,'LineWidth',linewidth); end end title('Solution of the equation system'); xlabel(xAxisLabel); ylabel(yAxisLabels(i,:)); legend(yAxisLabels(:,:)); hold off; end

end}

[Serap Karahanogullari]

It semmst too bad therefore i added .m file.

thanks one more time.