Solving a set of complex trig equation

Hello,

I have a very complex set of trig equations:

eqns1 = (-rake + (R_values(i)*(d/2))*deg2rad(skew)*tan(p) + (0.5 * chord_len(i) - x_c(1)) * sin(p) + y_u_L(1) * cos(p)) == xp(i*200);
eqns2 = R_values(i)*(d/2) * sind(skew -(180 * ((0.5 * chord_len(i) - x_c(1)) * cos(p) - y_u_L(1) * sin(p))) /(pi * R_values(i)*(d/2))) == yp(i*200);
eqns3 = R_values(i)*(d/2) * cosd(skew -(180 * ((0.5 * chord_len(i) - x_c(1)) * cos(p) - y_u_L(1) * sin(p))) /(pi * R_values(i)*(d/2))) == zp(i*200);

Attaching a picture of the equations for convineince of viewing.

Here I have named iG => rake, θs => skew, and θnt => p in my MATLAB script, The unknowns are: rake, p, and skew. All the other values are known.

I have tried using solve, vpasolve, and fsolve. solve and vpasolve generates an emptry structure for S.

Code of the loop that solves the equations:

for i = 1:9
    % Calculate airfoil properties
    chord_len = chord_length_at_R * d;
    max_c = max_camber_at_R .* chord_len;
    max_th = max_thickness_at_R * d;
    cam = -cam_d * max_c(i);
    th = th_d * max_th(i);
    x_c = X_c * chord_len(i);
    theta1 = atan2(der_y, 1);
    y_u_L = cam + th .* cos(theta1);
    % Define equations
    eqns1 = (-rake + (R_values(i)*(d/2))*deg2rad(skew)*tan(p) + ...
        (0.5 * chord_len(i) - x_c(1)) * sin(p) + y_u_L(1) * cos(p)) == xp(i*200);
    eqns2 = R_values(i)*(d/2) * sind(skew - ...
        (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(p) - y_u_L(1) * sin(p))) / ...
        (pi * R_values(i)*(d/2))) == yp(i*200);
    eqns3 = R_values(i)*(d/2) * cosd(skew - ...
        (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(p) - y_u_L(1) * sin(p))) / ...
        (pi * R_values(i)*(d/2))) == zp(i*200);
    % Solve equations
    try
        S = vpasolve(simplify([eqns1, eqns2, eqns3]), [p, skew, rake], [pitch(i),0, skew_angle_at_R(i)]);  % With initial guess
        if isempty(S)
            disp(['No solution for section ', num2str(i)]);
        else
            rake_sol = double(S.rake);
            p_sol = double(S.p);
            skew_sol = double(S.skew);
            fprintf('Section %d: Rake = %.4f, P = %.4f, Skew = %.4f\n', i, rake_sol, p_sol, skew_sol);
        end
    catch ME
        disp(['Error solving section ', num2str(i)]);
        disp(ME.message);
    end
end

And with fsolve I can only converge the first two sections.

Code of the loop that solves the equations:

for i = 1:n
    % Scale airfoil data for current section
    chord_len = chord_length_at_R * d;
    max_c = max_camber_at_R .* chord_len;
    max_th = max_thickness_at_R * d;
    cam = cam_d * max_c(i);
    th = th_d * max_th(i);
    x_c = X_c * chord_len(i);
    theta1 = atan2(der_y, 1);
    y_u_L = cam + th .* cos(theta1);
    % Define equations for fsolve
    fun = @(vars) [
        (-vars(2) + (R_values(i) * (d / 2)) * vars(3) * tan(vars(1)) + (0.5 * chord_len(i) - x_c(1)) * sin(vars(1)) + y_u_L(1) * cos(vars(1))) - xp(i * 200);
        R_values(i) * (d / 2) * sind(vars(3) - (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(vars(1)) - y_u_L(1) * sin(vars(1)))) / (pi * R_values(i) * (d / 2))) - yp(i * 200);
        R_values(i) * (d / 2) * cosd(vars(3) - (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(vars(1)) - y_u_L(1) * sin(vars(1)))) / (pi * R_values(i) * (d / 2))) - zp(i * 200)
        ];
    % Use pitch, skew_angle_at_R, and rake as initial guesses
    initial_guess = [pitch(i),0, skew_angle_at_R(i)];
    % Solve using fsolve
    options = optimoptions('fsolve', 'Display', 'none', 'MaxIterations', 10000, 'MaxFunctionEvaluations', 20000);
    [solution, fval, exitflag] = fsolve(fun, initial_guess, options);
    % Handle non-convergence
    if exitflag <= 0
        fprintf('fsolve did not converge for section %d. Trying alternative guesses...\n', i);
        alt_guesses = [rand(1, 3); [pitch(i),0, skew_angle_at_R(i)]]; % Use random or skew-related guesses
        for j = 1:size(alt_guesses, 1)
            [solution, fval, exitflag] = fsolve(fun, alt_guesses(j, :), options);
            if exitflag > 0
                break;
            end
        end
    end
    % Store solutions
    if exitflag > 0
        solutions(i, :) = solution;
    else
        fprintf('Failed to converge for section %d even with alternative guesses.\n', i);
    end
end

I would really appreacite any help on this as I've been stuck on this for some time. I know that these equations are highly non-linaer which is why MATLAB is probably having a hard time solving this, but I am not sure how else to approach this problem.

[Torsten]

Please supply executable code (with all parameters defined).

[Irtiza Masrur Khan]

Attached are the executable code and the excels.

clear;
clc;

% Load airfoil and coordinate data
airfoil_path = 'airfoil_data_fixed.xlsx';
coords_path = 'ORCA_Pointwise.xlsx';

% Read the airfoil data
Airfoil = xlsread(airfoil_path);
Coords = xlsread(coords_path);

X_c = Airfoil(:, 1);  % Chord coordinates (normalized)
y_c = Airfoil(:, 2);  % Camber distribution
der_y = Airfoil(:, 3);  % Slope of camber
th_d = Airfoil(:, 4);  % Thickness distribution

yp = Coords(:, 1);
zp = Coords(:, 2);
xp = Coords(:, 3);

% Normalize camber distribution
cam_d = y_c / max(y_c);

% Define polynomial functions
MaxCamber = @(x) -4448.8369*x.^12 + 30393.6831*x.^11 - 92977.6043*x.^10 + 168066.8833*x.^9 - 199490.6759*x.^8 + 163416.9800*x.^7 - 94493.8152*x.^6 + 38765.7447*x.^5 - 11176.4246*x.^4 + 2207.9559*x.^3 - 285.0846*x.^2 + 21.9443*x - 0.7490;
ChordLength = @(x) -143202.4761*x.^12 + 978274.9902*x.^11 - 2992184.0323*x.^10 + 5408923.9625*x.^9 - 6424276.8851*x.^8 + 5271614.6993*x.^7 - 3058632.5267*x.^6 + 1261908.7720*x.^5 - 366745.1423*x.^4 + 73096.4455*x.^3 - 9470.1908*x.^2 + 716.1619*x - 23.7157;
MaxThickness = @(x) -630.6311*x.^12 + 2601.5312*x.^11 - 2668.3834*x.^10 - 4627.0740*x.^9 + 16182.8469*x.^8 - 21144.9359*x.^7 + 15960.6510*x.^6 - 7582.2691*x.^5 + 2294.7299*x.^4 - 431.4940*x.^3 + 47.9318*x.^2 - 3.0647*x + 0.1614;
Pitch = @(x)  19344.5071*x.^12 - 114044.8587*x.^11 + 280789.2801*x.^10 - 357377.6146*x.^9 + 207947.2705*x.^8 + 43330.8173*x.^7 - 173099.4797*x.^6 + 143116.5772*x.^5 - 65570.9523*x.^4 + 18410.2055*x.^3 - 3128.7530*x.^2 + 294.0671*x - 10.4419;
SkewAngle = @(x) -719361.0309*x.^12 + 5176951.0162*x.^11 - 16746858.1600*x.^10 + 32161071.3897*x.^9 - 40784306.1013*x.^8 + 35930276.7571*x.^7 - 22515881.4272*x.^6 + 10096627.1844*x.^5 - 3209956.6167*x.^4 + 704272.6595*x.^3 - 100947.8092*x.^2 + 8462.8187*x - 312.8100;

% Scaling factor for the geometry
d = 1.4;

% Define the R values where we want to evaluate the polynomials
R_values = [0.1834, 0.2408, 0.2982, 0.3556, 0.413, 0.4704, 0.5278, 0.5852, 0.6426] / 0.7;

% Evaluate the polynomials at R_values
max_camber_at_R = MaxCamber(R_values);
pitch_dia_at_R = Pitch(R_values) * d * 1.4;
chord_length_at_R = ChordLength(R_values) * 0.99;
max_thickness_at_R = MaxThickness(R_values);
skew_angle_at_R = SkewAngle(R_values);
pitch = atan2(pitch_dia_at_R, (2 * pi * R_values)); % Angle in radians

chord_len = chord_length_at_R * d;
max_c = max_camber_at_R .* chord_len;
max_th = max_thickness_at_R * d;

% Number of sections (airfoils)
n = length(R_values);

% Solve for each section
solutions = zeros(n, 3);

for i = 1:n
    % Scale airfoil data for current section
    chord_len = chord_length_at_R * d;
    max_c = max_camber_at_R .* chord_len;
    max_th = max_thickness_at_R * d;

    cam = cam_d * max_c(i);
    th = th_d * max_th(i);
    x_c = X_c * chord_len(i);

    theta1 = atan2(der_y, 1);
    y_u_L = cam + th .* cos(theta1);
R_values(i)^2*d^2/4-yp(i*200)^2-zp(i*200)^2
    % Define equations for fsolve
    fun = @(vars) [
        (-vars(2) + (R_values(i) * (d / 2)) * vars(3) * tan(vars(1)) + (0.5 * chord_len(i) - x_c(1)) * sin(vars(1)) + y_u_L(1) * cos(vars(1))) - xp(i * 200);
        R_values(i) * (d / 2) * sind(vars(3) - (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(vars(1)) - y_u_L(1) * sin(vars(1)))) / (pi * R_values(i) * (d / 2))) - yp(i * 200);
        R_values(i) * (d / 2) * cosd(vars(3) - (180 * ((0.5 * chord_len(i) - x_c(1)) * cos(vars(1)) - y_u_L(1) * sin(vars(1)))) / (pi * R_values(i) * (d / 2))) - zp(i * 200)
    ];

    % Use pitch, skew_angle_at_R, and rake as initial guesses
    initial_guess = [pitch(i), skew_angle_at_R(i), 0];

    % Solve using fsolve
    options = optimoptions('fsolve', 'Display', 'none', 'MaxIterations', 10000, 'MaxFunctionEvaluations', 20000);
    [solution, fval, exitflag] = fsolve(fun, initial_guess, options);

    % Handle non-convergence
    if exitflag <= 0
        fprintf('fsolve did not converge for section %d. Trying alternative guesses...\n', i);
        alt_guesses = [rand(1, 3); [pitch(i), skew_angle_at_R(i), 0]]; % Use random or skew-related guesses
        for j = 1:size(alt_guesses, 1)
            [solution, fval, exitflag] = fsolve(fun, alt_guesses(j, :), options);
            if exitflag > 0
                break;
            end
        end
    end

    % Store solutions
    if exitflag > 0
        solutions(i, :) = solution;
        
    else
        fprintf('Failed to converge for section %d even with alternative guesses.\n', i);
    end
   
end

ans = 8.1922e-10

ans = 4.1339e-11

ans = 0.0309

fsolve did not converge for section 3. Trying alternative guesses...

Failed to converge for section 3 even with alternative guesses.

ans = 0.0375

fsolve did not converge for section 4. Trying alternative guesses...

Failed to converge for section 4 even with alternative guesses.

ans = 0.0816

fsolve did not converge for section 5. Trying alternative guesses...

Failed to converge for section 5 even with alternative guesses.

ans = 0.0948

fsolve did not converge for section 6. Trying alternative guesses...

Failed to converge for section 6 even with alternative guesses.

ans = 0.1521

fsolve did not converge for section 7. Trying alternative guesses...

Failed to converge for section 7 even with alternative guesses.

ans = 0.1719

fsolve did not converge for section 8. Trying alternative guesses...

Failed to converge for section 8 even with alternative guesses.

ans = 0.2424

fsolve did not converge for section 9. Trying alternative guesses...

Failed to converge for section 9 even with alternative guesses.

 disp(solutions);

    2.1293    0.2920    0.2256
    0.8014    0.1956    0.8773
         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0
         0         0         0

Answer 1

[Torsten]

Do you see the problem ? Your R-values and your y- and z-values are not compatible (at least for n = 3).

It must hold that

R_values(i)^2*d^2/4-yp(i*200)^2-zp(i*200)^2 = 0

for all i.

[Irtiza Masrur Khan]

You are right. My previous code was reading the xp,yp,zp values wrong from the excel files.