Filter löschen
Filter löschen

Matrix is singular to working precision.

1 Ansicht (letzte 30 Tage)
SEVEN JHUANG
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];
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 Kommentare
Sudarshan
Sudarshan am 8 Dez. 2022
Could you elaborate on the exact issue you are facing?
Jan
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.

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Image Analyst
Image Analyst am 8 Dez. 2022

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by