# Matrix is singular to working precision.

1 Ansicht (letzte 30 Tage)
SEVEN JHUANG am 30 Nov. 2022
Beantwortet: Image Analyst am 8 Dez. 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];
disp(M)
disp(n)
end
disp(binvl)
TTo= Aa(end);
disp(TTo);
Tt =[Tt , TTo];
end
xlswrite('C:\Users\NTU\Desktop\new data.xlsx',Tt','betaa')
end
##### 2 KommentareKeine anzeigenKeine ausblenden
Sudarshan am 8 Dez. 2022
Could you elaborate on the exact issue you are facing?
Jan am 8 Dez. 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.

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

Image Analyst am 8 Dez. 2022
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by