Matrix is singular to working precision.

3 views (last 30 days)
SEVEN JHUANG
SEVEN JHUANG on 30 Nov 2022
Answered: Image Analyst on 8 Dec 2022
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
  2 Comments
Jan
Jan on 8 Dec 2022
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.

Sign in to comment.

Answers (1)

Image Analyst
Image Analyst on 8 Dec 2022

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by