How to fix error with while loop?

1 Ansicht (letzte 30 Tage)
Katherine
Katherine am 21 Mär. 2023
Bearbeitet: Stephen23 am 21 Mär. 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);

Antworten (1)

Stephen23
Stephen23 am 21 Mär. 2023
Bearbeitet: Stephen23 am 21 Mär. 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

Kategorien

Mehr zu Elementary Math finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by