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.

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]

"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

See the "Tips" section here: https://www.mathworks.com/help/matlab/ref/mldivide.html#bthh_f9

Answer 1

[Simon Chan]

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)