MATLAB: Help solving an equation using fzero

fzerotraps

I need to solve an equation 
CT = (1/2)*sigma*r.^2*a.*(theta_75+(theta_tw*(r-0.75))-(lambda_i/r'))
by first using the trapz function to integrate this equation, and then using the fzero function to solve CT-CT_pres = 0, using theta_75 as an independent variable.
This is my code, I'm not sure if i'm using trapz correctly, and I can't get fzero to work at all, it's just one error after the other. Also i'm basically solving all this three times simultaniously by declaring theta_tw = [0,-10,-20] so I've tried to solve each theta_tw individually but still no luck.
Any help would be appreciated.
clear all; close all; clc;
Weight              = 20000;    % lbs
T                   = Weight;   % Thrust = Weight (lbs)
rho                 = 0.00237;  % Density (slug/ft^3)
R                   = 30;       % Radius (ft)
A                   = 2827;     % Area (ft^2)
v_tip               = 650;      % Tip Speed (ft/sec)
kappa               = 1.15;     % Induced Power Factor
Nb                  = 4;        % Number of Blades
c                   = 2;        % Chord (ft)
CD                  = 0.01;     % Drag Coefficient
sigma               = 0.085;    % Solidity
a                   = 6;        % Lift Curve Slope
%% Problem 2
% Rotor Disk Loading
DL = T/A;                   % lbs/ft^2
% Induced Velocity
v = sqrt(DL/(2*rho));       % ft/sec
% Ideal Induced Power
Pi_Ideal = (T*v)/550;       % hp
% Ideal Power Loading
PL_Ideal = T/Pi_Ideal;      % lbs/hp

% Non-Dimensional Coefficients at Hover
CT_pres = T/(rho*A*v_tip^2);            % Thrust
CPi = kappa*(CT_pres*sqrt(CT_pres/2));  % Induced Power
CP0 = (1/8)*sigma*CD;                   % Rotor Profile Power
CQ0 = CP0;                              % Torque 
CP = CPi+CP0;                           % Rotor Total Power
CQ = CP;                                % Torque
% Dimensional Coefficients at Hover
Pi = (CPi*rho*A*v_tip^3)/550;   % Induced Power (hp)
P0 = (CP0*rho*A*v_tip^3)/550;   % Profile Power (hp)
P = Pi+P0;                      % Rotor Total Power (hp)
% Figure of Merit
FM = Pi/P;
% Actual Power Loading
PL = T/P;                       % lbs/hp
fprintf(['----------- Problem 2 Answers -----------\n\n'...
    'Rotor Disk Loading = %f lbs/ft^2 \n'...
        'Ideal Induced Power = %f hp \nIdeal Power Loading = %f hp \n'...
            'Thrust Coefficient = %f \nTorque Coefficient = %f \n'...
                'Figure of Merit = %f \nProfile Power = %f hp \n'...
                    'Actual Power Loading = %f lbs/hp \n\n'],...
    DL, Pi_Ideal,PL_Ideal, CT_pres, CQ, FM, P0, PL);
%% Problem 3
lambda_i = sqrt(CT_pres/2); % Uniform Inflow
theta_tw = [0;-10;-20];     % Twist Rate (Degrees)
r = 0:0.1:1;                % Rotor section from root(0) to tip(1).
% Estimated Twist angle at 75% Raidus. (Degrees)
theta_75 = (((6*CT_pres)/(sigma*a))+((3/2)*sqrt(CT_pres/2)))
% Integrating for CT
dct = dCT(sigma,r,a,theta_75,theta_tw,lambda_i)'
ct = trapz(dct)
% Solving for theta_.75
fun = @(sigma,r,a,X,theta_tw,lambda_i)dCT;
X = 0;
fzero = (fun-Ct_pres,X); 
%theta_r = (((6*CT)/(sigma*a))+((3/2)*sqrt(CT/2)))+(theta_tw*(r-0.75))
function lambda = lambda_r(sigma,a,theta_r)
lambda = (1/16)*(sqrt(1+(32/(sigma*a))*theta_r)-1);
end
function CT = dCT(sigma,r,a,theta_75,theta_tw,lambda_i)
CT = (1/2)*sigma*r.^2*a.*(theta_75+(theta_tw*(r-0.75))-(lambda_i/r'));
end

Best Answer

Since you specifically asked for help on fzero, this works: 
X = 1;
for k = 1:numel(theta_tw)
    T_75(k) = fzero(@(theta_75)dCT(sigma,R,a,theta_75,theta_tw(k),lambda_i)-CT_pres, X) 
end
producing: 
T_75 =
       0.00201198547991055           292.50201198548           585.00201198548
Note that since ‘theta_tw’ is a vector, it’s necessary to iterate through its elements. The values for ‘T_75’ correspond to those eleements. Also, it is never appropriate to begin with an initial parameter estimate of 0, especially in a root-finding function such as fzero. I set ‘X’ to 1 for this reason.
The rest of your code appears to run without error with the fzero change I am posting here. You will need to determine if it produces the correct results. Note that it may be necessary to iterate through ‘T_75’ (name it whatever you want) if you are using non-vectorised code with its results.
Related Question