How to fix error with while loop?

1 visualizzazione (ultimi 30 giorni)
Katherine
Katherine il 21 Mar 2023
Modificato: Stephen23 il 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);

Risposte (1)

Stephen23
Stephen23 il 21 Mar 2023
Modificato: Stephen23 il 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

Categorie

Scopri di più su Elementary Math in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by