Matrix is singular to working precision.

function circular_moving_porous

Tt=[];    
clc;
format long


for  binvl=[0.1,0.2,0.3,0.4,0.5,0.6] %b/c
    lamda = 10;
    b =binvl;
    c = 1;
    a= 0;
    n=10;
    Count = 1; Aa=[0];
    
    
    
    
    while Count > 0.0001
        n = n + 2 ;
        [theta2] = collition2(n); % 0 ~ pi/82
        %disp(theta2);
        Vcount = zeros(4*n,1);
        fluidcal = zeros(4*n,4*n);
        
        k = 1;
        p=1;
        m=1;
        for i = 1:n  % i is angle 
            q=1;
            
            f=1;
            h=1;
            %disp(i)
            %disp('**')
            theta = theta2(:,i);
            
            %------------------------撘?13)蝑??喲?----------------------------------------
            
            %V = lamda*b*sinh(lamda*b)-cosh(lamda*b);
            W = 2*(lamda^3)*(b^3)+(lamda^3)*(a^3)+3*lamda*a;
            F1 = W*(sinh(lamda*b)-cosh(lamda*b));
            F2 = -(3*sinh(lamda*b)+(W+3*lamda*b*cosh(lamda*b)));
            
            Vcount(4*i-3,1) = 0;
            Vcount(4*i-2,1) = -(2/3)*lamda*(F2/F1)/(10^(-10));
            Vcount(4*i-1,1) = 0;
            Vcount(4*i,1) = 0;
            
            for j=2:2:2*n  % j is n 
                %disp("j=");
                %disp(j);
                
                
                
                %------------------------摰儔Polynomial------------------------------------
                
                %-----Gegenbauer and legendre---------------------------------
                P =  legendre(j,cos(theta)) ; 
                P1_subtra = legendre(j-1,cos(theta)) ;
                P2 = legendre(j-2,cos(theta));
                P1_plus= legendre(j+1,cos(theta)) ;
                G1_plus= (P1_subtra(1,:)-P1_plus(1,:))/(2*j+1); 
                G  = (P2(1,:)-P(1,:))/(2*j-1) ;
                %-----A--------------------------------
                A11 = @(w) -((-1)^(j/2))*( (2*(w^j))/(pi*factorial(j)) )*besselk(1,c*w);   %OK
                A21 = @(w) ((-1)^(j/2))*( (2*(w^(j-2)))/(pi*factorial(j)) )*((j-2)*(j-3)*besselk(1,c*w)-(2*j-3)*c*w*besselk(0,c*w));   %OK
                A12 = @(w) ((-1)^(j/2))*( (2*(w^j))/(pi*factorial(j)) )*besselk(0,c*w);    %OK
                A22 = @(w) ((-1)^(j/2))*( (2*(w^(j-2)))/(pi*factorial(j)) )*((2*j-3)*c*w*besselk(1,c*w)-j*(j-1)*besselk(0,c*w));    %OK
                %-----S---------------------------------
                S11 = @(w) (A11(w)*besseli(0,c*w)-A12(w)*besseli(1,c*w))/(c*w*(besseli(0,c*w)^2)-2*besseli(0,c*w)*besseli(1,c*w)-c*w*(besseli(1,c*w)^2));   %OK
                S21 = @(w) (A21(w)*besseli(0,c*w)-A22(w)*besseli(1,c*w))/(c*w*(besseli(0,c*w)^2)-2*besseli(0,c*w)*besseli(1,c*w)-c*w*(besseli(1,c*w)^2));  %OK
                S12 = @(w) (A11(w)-c*w*besseli(0,c*w)*S11(w))/(besseli(1,c*w)) ;   %OK
                S22 = @(w) (A21(w)-c*w*besseli(0,c*w)*S21(w))/(besseli(1,c*w)) ;   %OK
                %----------------gamma(b,thata)--------------------
                r11_n = @(w) (S11(w)*besseli(0,w*b*sin(theta))*w*b*sin(theta)+S12(w)*besseli(1,w*b*sin(theta)))*sin(w*b*cos(theta));
                %r11_n = @(w) S11(w)*besseli(0,w*b*sin(theta))*w*b*sin(theta)+S12(w)*besseli(1,w*b*sin(theta));
                %R11_n = @(w) r11_n(w)*sin(w*b*cos(theta));
                R11_n = @(w) r11_n(w);
                R11_int = gauss_laguerre180(R11_n);
                
                r11 = R11_int-(b^(-j-1))*((j+1)*G1_plus*csc(theta)) ;  
                
                %r21_n = @(w) S21(w)*besseli(0,w*b*sin(theta))*w*b*sin(theta)+S22(w)*besseli(1,w*b*sin(theta));
                %R21_n = @(w) r21_n(w)*sin(w*b*cos(theta));
                r21_n = @(w) (S21(w)*besseli(0,w*b*sin(theta))*w*b*sin(theta)+S22(w)*besseli(1,w*b*sin(theta)))*sin(w*b*cos(theta));
                R21_n = @(w) r21_n(w);
                
                R21_int = gauss_laguerre180(R21_n);
                %R21_int = integral(R21_n,0,Inf);
                r21 = R21_int - (b^(-j+1))*((j+1)*G1_plus*csc(theta)-2*G*cot(theta));  
                
                r12_n = @(w) S21(w)*(2*besseli(0,w*b*sin(theta))+besseli(1,w*b*sin(theta))*w*b*sin(theta))+S22(w)*besseli(0,w*b*sin(theta));
                R12_n = @(w) r12_n(w)*cos(w*b*cos(theta)) ;         
                R12_int = gauss_laguerre180(R12_n);
                r12 = R12_int - (b^(-j-1))*P(1,:); 
                
                r22_n = @(w) S21(w)*(2*besseli(0,w*b*sin(theta))+besseli(1,w*b*sin(theta))*w*b*sin(theta))+S22(w)*besseli(0,w*b*sin(theta));
                R22_n = @(w) r22_n(w)*cos(w*b*cos(theta)) ;
                R22_int = gauss_laguerre180(R22_n);
                r22 = R22_int - (b^(-j+1))*(P(1,:)+2*G*cot(theta));       
                
                %----------------gamma(b,thata)-star--------------------
                C_star = @(w) -(w)*(cos(2*theta)*cos(w*b*cos(theta))*besseli(1,w*b*sin(theta))-(1/4)*w*sin(2*theta)*sin(w*b*cos(theta))*(3*besseli(0,w*b*sin(theta))+besseli(2,w*b*sin(theta))) );
                D_star = @(w) 2*w*((sin(theta))^2)*besseli(1,w*b*sin(theta))*(cos(w*b*cos(theta))-w*b*cos(theta)*sin(w*b*cos(theta)))-(w*b*cos(2*theta)*cos(w*b*cos(theta))+3*cos(theta)*sin(w*b*cos(theta)))*besseli(0,w*b*sin(theta))*w*sin(theta);
                C_double_star = @(w) -(1/4)*(w)*sin(w*b*cos(theta))*besseli(0,w*b*sin(theta))*(1+3*cos(2*theta))+2*(w)*sin(theta)*cos(theta)*cos(w*b*cos(theta))*besseli(1,w*b*sin(theta))+(1/2)*(w)*((sin(theta))^2)*sin(w*b*cos(theta))*besseli(2,w*b*sin(theta));
                D_double_star = @(w) (sin(w*b*cos(theta))+3*cos(2*theta)*sin(w*b*cos(theta))-4*(sin(theta)^2)*cos(theta)*cos(w*b*cos(theta)) )*besseli(0,w*b*sin(theta))*(w/2) + ...
                    (cos(2*theta)*sin(w*b*cos(theta))*w*b-2*cos(w*b*cos(theta))*cos(theta))*besseli(1,w*b*sin(theta))*w*sin(theta);
                r1_star = @(w) S12(w)*C_star(w)+S11(w)*D_star(w);
                R1_star_int1 = gauss_laguerre180(r1_star);
                R1_star_int = R1_star_int1 +(1-j)*(-1-j)*(b^(-2-j))*G*csc(theta);
                r2_star = @(w) S22(w)*C_star(w)+S21(w)*D_star(w);
                R2_star_int1 = gauss_laguerre180(r2_star);
                R2_star_int = R2_star_int1 +(3-j)*(1-j)*(b^(-j))*G*csc(theta);
                
                r1_double_star = @(w) S12(w)*C_double_star(w)*S11(w)*D_double_star(w);
                R1_double_star1= gauss_laguerre180(r1_double_star );
                R1_double_star = R1_double_star1 -(-1-j)*(b^(-2-j))*P1_subtra(1,:);
                r2_double_star = @(w) S22(w)*C_double_star(w)*S21(w)*D_double_star(w);
                R2_double_star1= gauss_laguerre180(r2_double_star );
                R2_double_star = R2_double_star1 -(1-j)*(b^(-j))*P1_subtra(1,:);
                
                
                %------------------------摰儔?寧?蝯???-------------------------------------
                
                
                B3 = r11;
                
                D3 = r21;
                
                A3 = -(b^(j-2))*((j+1)*G1_plus*csc(theta)-(2*j-1)*G*cot(theta));
                
                
                
                C3 = -(lamda^(1/2))*(b^(-3/2))*( ((j+1)*csc(theta)*G1_plus-(2*j-1)*cot(theta)*G)*besseli(j-1/2,lamda*b)-lamda*b*cot(theta)*G*besseli(j+1/2,lamda*b) );
                
                
                
                B4 = r12;
                
                D4 = r22;
                
                A4 = -(b^(j-2))*((2*j-1)*G+P(1,:));
                
                
                
                C4 = -(lamda^(1/2))*(b^(-3/2))*( (P1_plus(1,:)+(j-1)*G1_plus+j*G)*besseli(j-1/2,lamda*b)+lamda*b*G*besseli(j+1/2,lamda*b) );
                
                
                B5 = R1_star_int;
                
                D5 = R2_star_int;
                
                A5 = (b^(j-3))*j*(j-2)*csc(theta)*G;
                
                
                
                C5 = (lamda^(-1/2))*(b^(-7/2))*( ((j-2)*(2*j+1)+2*(lamda*b)^2)*j*besseli(j+1/2,lamda*b)+lamda*b*(j^2-2*j+(lamda*b)^2)*besseli(j+3/2,lamda*b) )*G*csc(theta);  %ok
                
                
                
                B6 = R1_double_star;
                
                D6 = R2_double_star;
                
                A6 = -(b^(j-1))*(j-2)*P1_subtra(1,:);  
                
                
                
                C6 = (lamda^(-1/2))*(b^(-7/2))*( ((-j+2)*(2*j+1)-2*(lamda*b)^2)*j*besseli(j+1/2,lamda*b)+lamda*b*(j^2-2*j-(lamda*b)^2)*besseli(j+3/2,lamda*b) )*G*csc(theta);  %ok
                
                %Dd6 = (lamda^(1/2))*(b^(-5/2))*( (2-j)*besselk(j-1/2,lamda*b)-lamda*b*besselk(j+1/2,lamda*b) )*P1_subtra(1,:);
                
                %------------------------撘?13)蝑?撌阡?------------------------------------------
                if j==2
                    fluidcal(4*i-3,4*j-7) = B3/(10^(-10));
                    fluidcal(4*i-3,4*j-6) = D3/(10^(-10));
                    fluidcal(4*i-3,4*j-5) = -A3/(10^(-10));
                    
                    fluidcal(4*i-3,4*j-4) = -C3/(10^(-10));
                    
                    fluidcal(4*i-2,4*j-7) = B4/(10^(-10));
                    fluidcal(4*i-2,4*j-6) = D4/(10^(-10));
                    fluidcal(4*i-2,4*j-5) = -A4/(10^(-10));
                    
                    fluidcal(4*i-2,4*j-4) = -C4/(10^(-10));
                    
                    fluidcal(4*i-1,4*j-7) = B5/(10^(-10));
                    fluidcal(4*i-1,4*j-6) = D5/(10^(-10));
                    fluidcal(4*i-1,4*j-5) = -A5/(10^(-10));
                    
                    fluidcal(4*i-1,4*j-4) = -C5/(10^(-10));
                    
                    fluidcal(4*i,4*j-7) = B6/(10^(-10));
                    fluidcal(4*i,4*j-6) = D6/(10^(-10));
                    fluidcal(4*i,4*j-5) = -A6/(10^(-10));
                    
                    fluidcal(4*i,4*j-4) = -C6/(10^(-10));
                    
                else
                    fluidcal(4*i-3,4*j-7-4*q) = B3/(10^(-10));
                    fluidcal(4*i-3,4*j-6-4*q) = D3/(10^(-10));
                    fluidcal(4*i-3,4*j-5-4*q) = -A3/(10^(-10));
                    
                    fluidcal(4*i-3,4*j-4-4*q) = -C3/(10^(-10));
                    
                    fluidcal(4*i-2,4*j-7-4*q) = B4/(10^(-10));
                    fluidcal(4*i-2,4*j-6-4*q) = D4/(10^(-10));
                    fluidcal(4*i-2,4*j-5-4*q) = -A4/(10^(-10));
                    
                    fluidcal(4*i-2,4*j-4-4*q) = -C4/(10^(-10));
                    
                    fluidcal(4*i-1,4*j-7-4*q) = B5/(10^(-10));
                    fluidcal(4*i-1,4*j-6-4*q) = D5/(10^(-10));
                    fluidcal(4*i-1,4*j-5-4*q) = -A5/(10^(-10));
                    
                    fluidcal(4*i-1,4*j-4-4*q) = -C5/(10^(-10));
                    
                    fluidcal(4*i,4*j-7-4*q) = B6/(10^(-10));
                    fluidcal(4*i,4*j-6-4*q) = D6/(10^(-10));
                    fluidcal(4*i,4*j-5-4*q) = -A6/(10^(-10));
                    
                    fluidcal(4*i,4*j-4-4*q) = -C6/(10^(-10));
                    
                    q=q+1;
                end
                %              
            end
        end
        
        Torque=linsolve(fluidcal,Vcount);
        
        M=Torque(2,1);
        
        Aa =[Aa , M]; 
        adj=fix(log10(abs(Aa(:,length(Aa)))));
        
        disp(M) 
        disp(n) 
        Count = abs((Aa(:,length(Aa))-Aa(:,length(Aa)-1))/(10^(adj)));
    end
    
    disp(binvl)
    TTo= Aa(end);
    disp(TTo);
    Tt =[Tt , TTo]; 
    
    
    
    
end  

xlswrite('C:\Users\NTU\Desktop\new data.xlsx',Tt','betaa')

end

[Sudarshan]

Could you elaborate on the exact issue you are facing?

[Jan]

This is a huge code. You should not expect the readers to understand it directly. I cannot run your code, because the function collition2 is missing.

Please post at least the complete error message to show the readers, which line is failing.

By the way, 10^(-10) is an expensive power operation, while 1e-10 is a cheap constant.

Answer 1

[Image Analyst]

FAQ:

https://matlab.fandom.com/wiki/FAQ#Does_MATLAB_only_calculate_to_4_significant_digits?