Clear Filters
Clear Filters

How to fix error with while loop?

4 views (last 30 days)
Katherine
Katherine on 21 Mar 2023
Edited: Stephen23 on 21 Mar 2023
I am working on creating a code to calculate the catalyst weight necessary to achieve 60% conversion and will plot the X,y,f and reaction rate as a function of catalyst weight. I keep getting the error shown below. How can I fix this? Thank you!
k = 0.00392; % mol/(atm kgcat sec)
FA0 = 0.1362; % mol/sec
FB0 = 0.068; % mol/sec
P0 = 10; % atm
alpha = 0.0367; % kg
X_upper = 1; % upper bound for X
X_lower = 0; % lower bound for X
epsilon = 0.4; % given constant
conversion = 0.6; % desired conversion
tolerance = 1e-6; % tolerance for bisection method
function res = ethylene_oxide_residual(X)
y = P0 / (P0 * (1 - 0.5*X)^(2/3) * (1 + epsilon*X));
r = k * FA0^(1/3) * (y/(1 + epsilon*X))^2/3 * (1 - X);
conversion_calculated = (1 - X) / FA0;
res = conversion_calculated - conversion;
end
while (X_upper - X_lower) > tolerance
Function definitions in a script must appear at the end of the file.
Move all statements after the "ethylene_oxide_residual" function definition to before the first local function definition.
X_guess = (X_upper + X_lower) / 2;
if ethylene_oxide_residual(X_guess) > 0
X_upper = X_guess;
else
X_lower = X_guess;
end
end
catalyst_weight = alpha * X_guess;
fprintf('Catalyst weight necessary for 60%% conversion: %.2f kg\n', catalyst_weight);
catalyst_weights = linspace(0, alpha, 100);
X_values = zeros(size(catalyst_weights));
y_values = zeros(size(catalyst_weights));
f_values = zeros(size(catalyst_weights));
reaction_rates = zeros(size(catalyst_weights));
for i = 1:length(catalyst_weights)
X = catalyst_weights(i) / alpha;
y = P0 / (P0 * (1 - 0.5*X)^(2/3) * (1 + epsilon*X));
f = (1 + epsilon*X) / y;
r = k * FA0^(1/3) * (y/(1 + epsilon*X))^2/3 * (1 - X);
X_values(i) = X;
y_values(i) = y;
f_values(i) = f;
reaction_rates(i) = r;
end
We can then plot these values using the following code.
subplot(2, 2, 1);
plot(catalyst_weights, X_values);
xlabel('Catalyst weight (kg)');
ylabel('Conversion (X)');
subplot(2, 2, 2);
plot(catalyst_weights, y_values);

Answers (1)

Stephen23
Stephen23 on 21 Mar 2023
Edited: Stephen23 on 21 Mar 2023
"I keep getting the error shown below. How can I fix this? "
By doing exactly what the error message tells you: move the script code to before the function definition:
k = 0.00392; % mol/(atm kgcat sec)
FA0 = 0.1362; % mol/sec
FB0 = 0.068; % mol/sec
P0 = 10; % atm
alpha = 0.0367; % kg
X_upper = 1; % upper bound for X
X_lower = 0; % lower bound for X
epsilon = 0.4; % given constant
conversion = 0.6; % desired conversion
tolerance = 1e-6; % tolerance for bisection method
while (X_upper - X_lower) > tolerance
X_guess = (X_upper + X_lower) / 2;
if ethylene_oxide_residual(X_guess,P0,epsilon,k,FA0,conversion) > 0
X_upper = X_guess;
else
X_lower = X_guess;
end
end
catalyst_weight = alpha * X_guess;
fprintf('Catalyst weight necessary for 60%% conversion: %.2f kg\n', catalyst_weight);
Catalyst weight necessary for 60% conversion: 0.00 kg
catalyst_weights = linspace(0, alpha, 100);
X_values = zeros(size(catalyst_weights));
y_values = zeros(size(catalyst_weights));
f_values = zeros(size(catalyst_weights));
reaction_rates = zeros(size(catalyst_weights));
for i = 1:length(catalyst_weights)
X = catalyst_weights(i) / alpha;
y = P0 / (P0 * (1 - 0.5*X)^(2/3) * (1 + epsilon*X));
f = (1 + epsilon*X) / y;
r = k * FA0^(1/3) * (y/(1 + epsilon*X))^2/3 * (1 - X);
X_values(i) = X;
y_values(i) = y;
f_values(i) = f;
reaction_rates(i) = r;
end
subplot(2, 2, 1);
plot(catalyst_weights, X_values);
xlabel('Catalyst weight (kg)');
ylabel('Conversion (X)');
subplot(2, 2, 2);
plot(catalyst_weights, y_values);
% here is the function, at the end:
function res = ethylene_oxide_residual(X,P0,epsilon,k,FA0,conversion)
y = P0 / (P0 * (1 - 0.5*X)^(2/3) * (1 + epsilon*X));
r = k * FA0^(1/3) * (y/(1 + epsilon*X))^2/3 * (1 - X);
conversion_calculated = (1 - X) / FA0;
res = conversion_calculated - conversion;
end

Categories

Find more on Elementary Math in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!