# Inverse of C matrix is Inf when it shouldn't be. How do I fix the matrix? The formula to calculate C is correct. May be an issue with input c.

2 views (last 30 days)
Shawn Alcorn on 5 Mar 2022
Answered: Simon Chan on 5 Mar 2022
This code gives me an INV(C) as INF. I believe the error is within the input of c into the calculation for C matrix. When I look at the C matrix the numbers do no look unusual enough to cause the inverse to be INF.
clc; clear all; close all
N=10;
V_inf = 50;
lambda = .5; %linspace(0,1,N)
S = 39.2699;
AR = 10.1859;
AoA = 8; % Angle of Attack in degrees
b = 20;
alpha_0 = 0; % Zero lift angle of attack is 0 due to airfoil being symmetric
cr = (2*b)/(AR*(lambda+1));
ct = cr*lambda;
i = 1;
alpha = AoA * (pi/180);
alpha_0 = alpha_0 * ones(N,1);
for theta = linspace(.01,pi-.01,20);
y = (b/2)*cos(theta);
if theta <= pi/2
c = (2/b)*cr*(1-lambda)*y+cr;
else
c = -(2/b)*cr*(1-lambda)*abs(y)+cr;
end
if length(alpha) == 1
alpha = alpha * ones(N,1);
end
if length(c) == 1
c = c * ones(N,1);
end
C = zeros(N,N);
B = zeros(N,1);
for j = 1:N
for n = 1:N
C(j,n) = 2*b/(pi*c(j)) * sin(n*theta) + n * sin(n*theta)/sin(theta);
end
B(j,1) = alpha(j,1) - alpha_0(j,1);
end
A = inv(C)*B;
Gamma = zeros(N,1);
alpha_i = zeros(N,1);
for j = 1:N
for n = 1:N
Gamma(j,1) = Gamma(j,1) + 2*b*V_inf * A(n) * sin(n*theta);
alpha_i(j,1) = alpha_i(j,1) + n * A(n) * sin(n*theta)/sin(theta);
end
end
end
A1(i) =A(1);
Cl(i) = AR*pi*A1;
Stephen23 on 5 Mar 2022
"Inverse of C matrix is Inf when it shouldn't be"
What do you expect the inverse of a non-invertible matrix to be?
"When I look at the C matrix the numbers do no look unusual enough to cause the inverse to be INF."
When I look at the C matrix all of the rows are copies of each other, which means that all of the rows are linearly dependent, which of course means that the matrix has no inverse, i.e. is not invertible.
Perhaps you looked at another matrix C.
So far neither INV nor MATLAB seem to be doing anything unexpected.
N=10;
V_inf = 50;
lambda = .5; %linspace(0,1,N)
S = 39.2699;
AR = 10.1859;
AoA = 8; % Angle of Attack in degrees
b = 20;
alpha_0 = 0; % Zero lift angle of attack is 0 due to airfoil being symmetric
cr = (2*b)/(AR*(lambda+1));
ct = cr*lambda;
i = 1;
alpha = AoA * (pi/180);
alpha_0 = alpha_0 * ones(N,1);
for theta = 0.01 % linspace(.01,pi-.01,20);
y = (b/2)*cos(theta);
if theta <= pi/2
c = (2/b)*cr*(1-lambda)*y+cr;
else
c = -(2/b)*cr*(1-lambda)*abs(y)+cr;
end
if length(alpha) == 1
alpha = alpha * ones(N,1);
end
if length(c) == 1
c = c * ones(N,1);
end
display(c)
C = zeros(N,N);
B = zeros(N,1);
for j = 1:N
for n = 1:N
C(j,n) = 2*b/(pi*c(j)) * sin(n*theta) + n * sin(n*theta)/sin(theta);
end
B(j,1) = alpha(j,1) - alpha_0(j,1);
end
display(C)
display(inv(C))
A = inv(C)*B
...
end
c = 10×1
3.9269 3.9269 3.9269 3.9269 3.9269 3.9269 3.9269 3.9269 3.9269 3.9269
C = 10×10
1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588 1.0324 4.0646 9.0961 16.1257 25.1520 36.1734 49.1876 64.1919 81.1835 100.1588
Warning: Matrix is singular to working precision.
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
Warning: Matrix is singular to working precision.
A = 10×1
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
"The formula to calculate C is correct"
And yet clearly it isn't. Rather than relying on just believing that your code is correct, now is a good time to actually learn to look at what your code is really doing, not just what you imagine/hope/want/believe it to be doing. Actually looking at what code is really doing is an important step in debugging your code.
For example, your code defines c to be a column vector of one repeated value:
if length(c) == 1
c = c * ones(N,1);
end
and there is nothing else that depends on the loop iterator j when you calculate C, so clearly every row of C will also be exactly the same, based on that repeated value in c.
Note: of course it won't fix C's linearly dependent rows (and hence being non-invertible), but rather than using INV and MTIMES you should use the more robust MLDIVIDE operator:
A = C\B
Warning: Matrix is singular to working precision.
A = 10×1
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Simon Chan on 5 Mar 2022
I have no idea about the details of the calculation but spot some of the issues stated as follows:
clc; clear all; close all
N=20; % Set it to 20, equals to number of theta
V_inf = 50;
lambda = .5; %linspace(0,1,N)
S = 39.2699;
AR = 10.1859;
AoA = 8; % Angle of Attack in degrees
b = 20;
alpha_0 = 0; % Zero lift angle of attack is 0 due to airfoil being symmetric
cr = (2*b)/(AR*(lambda+1));
ct = cr*lambda;
i = 1;
alpha = AoA * (pi/180);
alpha_0 = alpha_0 * ones(N,1);
Ntheta = 20; % Number of theta = 20
theta = linspace(.01,pi-.01,Ntheta); % Set the value of theta
for k = 1:Ntheta % Calculate each theta
y = (b/2)*cos(theta(k));
if theta(k) <= pi/2
c = (2/b)*cr*(1-lambda)*y+cr;
else
c = -(2/b)*cr*(1-lambda)*abs(y)+cr;
end
if length(alpha) == 1
alpha = alpha * ones(N,1);
end
if length(c) == 1
c = c * ones(N,1);
end
C = zeros(Ntheta,N); % Should be Ntheta instead of N
B = zeros(Ntheta,1); % Should be Ntheta instead of N
for j = 1:Ntheta
for n = 1:N
C(j,n) = 2*b/(pi*c(j)) * sin(n*theta(j)) + n * sin(n*theta(j))/sin(theta(j));
end
B(j,1) = alpha(j,1) - alpha_0(j,1);
end
A = inv(C)*B;
Gamma = zeros(Ntheta,1); % Should be Ntheta instead of N
alpha_i = zeros(Ntheta,1); % Should be Ntheta instead of N
for j = 1:Ntheta
for n = 1:N
Gamma(j,1) = Gamma(j,1) + 2*b*V_inf * A(n) * sin(n*theta(j));
alpha_i(j,1) = alpha_i(j,1) + n * A(n) * sin(n*theta(j))/sin(theta(j));
end
end
end
plot(theta,Gamma) 