How to fix error with while loop?

조회 수: 1 (최근 30일)
Katherine
Katherine 2023년 3월 21일
편집: Stephen23 2023년 3월 21일
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);

답변 (1개)

Stephen23
Stephen23 2023년 3월 21일
편집: Stephen23 2023년 3월 21일
"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

카테고리

Help CenterFile Exchange에서 Elementary Math에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by