What do I have to fix the plot to make it work?

1 回表示 (過去 30 日間)
Maria Ruiz
Maria Ruiz 2021 年 11 月 5 日
編集済み: Alberto Cuadra Lara 2021 年 11 月 5 日
%Euler's Method range from t=0 a 50 hrs
%consider x(0) = 1.25, S(0) = 86.63 and P(0) = 0.
%t = 0:25:50;
t=linspace(0,25,50);
%h =t(2)-t(2);
h=25;
%Initial values
xI(1) =1.25;
SI(1) =86.63;
PI(1) =0;
Kd = 0.0032;
a = 0.6212;
Yps = 0.5820;
Ki = 243.95
%Eq for mu
MU = 0.5557;
Ks = 0.0191;
p = 30.7600
%Aplicattion of eulers explicit method
for i = 1:length(t)
mu = (MU*SI(i)/Ks+SI(i)+(SI(i)^2/Ki))*((1-PI(i)/p)^1);
xI(i+1) =xI(i) +h*(mu*xI(i) - Kd*xI(i));
SI(i+1) =SI(i) +h*((-a/Yps)*(mu*xI(i)));
PI(i+1) =PI(i) +h*(a*mu*xI(i));
end
plot(t,xI(1:end-1))
hold on
plot(t,SI(1:end-1))
hold on
plot(t,PI(1:end-1))
Error using plot
Vectors must be the same length.

採用された回答

Walter Roberson
Walter Roberson 2021 年 11 月 5 日
The code works for me.
My guess is that you had an existing larger xI, SI,or PI in your workspace. You did not use functions for your code, so any existing variables in your workspace would have been inherited.
%Euler's Method range from t=0 a 50 hrs
%consider x(0) = 1.25, S(0) = 86.63 and P(0) = 0.
%t = 0:25:50;
t=linspace(0,25,50);
%h =t(2)-t(2);
h=25;
%Initial values
xI(1) =1.25;
SI(1) =86.63;
PI(1) =0;
Kd = 0.0032;
a = 0.6212;
Yps = 0.5820;
Ki = 243.95
Ki = 243.9500
%Eq for mu
MU = 0.5557;
Ks = 0.0191;
p = 30.7600
p = 30.7600
%Aplicattion of eulers explicit method
for i = 1:length(t)
mu = (MU*SI(i)/Ks+SI(i)+(SI(i)^2/Ki))*((1-PI(i)/p)^1);
xI(i+1) =xI(i) +h*(mu*xI(i) - Kd*xI(i));
SI(i+1) =SI(i) +h*((-a/Yps)*(mu*xI(i)));
PI(i+1) =PI(i) +h*(a*mu*xI(i));
end
size(t)
ans = 1×2
1 50
size(xI)
ans = 1×2
1 51
size(SI)
ans = 1×2
1 51
size(PI)
ans = 1×2
1 51
plot(t,xI(1:end-1))
hold on
plot(t,SI(1:end-1))
hold on
plot(t,PI(1:end-1))

その他の回答 (2 件)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021 年 11 月 5 日
It is working ok now, but your formulations have some problems. Thus, the calculated results (xI, SI, PI) are not OK.
clc; clearvars; close all
%Euler's Method range from t=0 a 50 hrs
%consider x(0) = 1.25, S(0) = 86.63 and P(0) = 0.
t=linspace(0,25,50);
h=25;
%Initial values
xI =1.25;
SI =86.63;
PI =0;
Kd = 0.0032;
a = 0.6212;
Yps = 0.5820;
Ki = 243.95;
%Eq for mu
MU = 0.5557;
Ks = 0.0191;
p = 30.7600;
%Aplicattion of eulers explicit method
for i = 1:length(t)
mu = (MU*SI(i)/Ks+SI(i)+(SI(i)^2/Ki))*((1-PI(i)/p)^1);
xI(i+1) =xI(i) +h*(mu*xI(i) - Kd*xI(i));
SI(i+1) =SI(i) +h*((-a/Yps)*(mu*xI(i)));
PI(i+1) =PI(i) +h*(a*mu*xI(i));
end
plot(t,xI(1:end-1), 'r-', 'linewidth', 1.5, 'DisplayName', 'xI')
hold on
plot(t,SI(1:end-1), 'g', 'linewidth', 1.5, 'DisplayName', 'SI')
plot(t,PI(1:end-1), 'b', 'linewidth', 1.5, 'DisplayName', 'PI')
grid on; legend
hold off

Alberto Cuadra Lara
Alberto Cuadra Lara 2021 年 11 月 5 日
編集済み: Alberto Cuadra Lara 2021 年 11 月 5 日
As @Sulaymon Eshkabilov has said, the formulation has some kind of typo. Check the definition of mu first, e.g., is the last term to the power of one?. The value of mu is what is blowing up the solution.
Here is the code with some comments. I encourage you to always try to include them for better readability.
% Euler's Explicit Method to solve numerically a Initial Value Problem
% Range from t= 0 to 50 hrs
% Initial conditions:
% * x(0) = 1.25,
% * S(0) = 86.63,
% * P(0) = 0.
% 1. Create mesh
Npoints = 200;
a = 0; % Initial mesh value
b = 50; % Final mesh value
t = linspace(a, b, Npoints);
% 2. Set step size
h = (b-a)/(Npoints - 1);
% 3. Define constants of functions
Kd = 0.0032;
a = 0.6212;
Yps = 0.5820;
Ki = 243.95;
MU = 0.5557;
Ks = 0.0191;
p = 30.7600;
% 4. Preallocate variables so as not to increase their size in each iteration
xI = zeros(1, Npoints);
SI = zeros(1, Npoints);
PI = zeros(1, Npoints);
% 5. Set the initial conditions for the Initial Value Problem
xI(1) = 1.25;
SI(1) = 86.63;
PI(1) = 0;
% 6. Use the Euler explicit/forward method to solve numerically
for i = 1:Npoints-1
mu = (MU * SI(i)/Ks + SI(i) + (SI(i)^2 / Ki)) * (1 - PI(i)/p)^-1;
xI(i+1) = xI(i) + h * (mu * xI(i) - Kd * xI(i));
SI(i+1) = SI(i) + h * (-a/Yps * (mu * xI(i)));
PI(i+1) = PI(i) + h * (a * mu * xI(i));
end
% 7. Plot results
figure; hold on;
plot(t, xI, 'LineWidth', 1.5)
plot(t, SI, 'LineWidth', 1.5)
plot(t, PI, 'LineWidth', 1.5)
legend({'xI', 'SI', 'PI'}, 'location', 'northeastoutside')

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

タグ

製品


リリース

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by