I want to comapre my A,B and P against time and want to compare my ODE89 function with Eulers method with a timestep of 3600

4 次查看（过去 30 天）
Naveen Krish
Naveen Krish 2022-3-14
评论： Naveen Krish 2022-3-14
nsteps = 12;
t = zeros (nsteps,1);
A = zeros (nsteps,1);
B = zeros(nsteps, 1);
P = zeros(nsteps,1);
A0 = 1;
B0 = 3;
P0 = 0;
ti=0;
tf=12*3600;
Yb = 1;
Yp = 0.15;
K = 5*10^-5;
h=3600;
timestep=[3600 1800 900 450 225];
[t,y]= euler (ti tf,(A0;B0;P0),timestep);
A_t=(A0-B0/YB)/(1-(B0/(YB*A0))*exp(-((YB*A0/B0)-1)*K*B0*t));
for k = 2:13
for j=1:length(timestep)
t(k) = t(k-1)+timestep(j)
A(k) = A(k-1)+(-K*A(k-1)*B(k-1))*timestep(j);
B(k) = B(k-1)+(-Yb*(K*A(k-1)*B(k-1)))*timestep(j);
P(k)= P(k-1)+ Yp*(K*A(k-1)*B(k-1))*timestep(j);
end
end
plot (t,A)
figure (1)
plot(t,A(:,1))
plot (t,B)
figure (2)
plot(t,B(:,1))
plot(t,P)
figure (3)
plot(t,P(:,1))
% e = abs(A_t - y(1)/y(1))*100;
% plot (e)
% hold on
%end
% title('Error compared with Analytical A value');
% xlabel ('Time (t)');
% ylabel (Error(%)');
% legend ('A_e 3600','A_e 1800','A_e 900','A_e 450','A_e 225')
% end
  5 个评论
显示 3更早的评论
Torsten
Torsten 2022-3-14
Why does "timestep" not have as many elements as the k-loop requires, namely 13-2+1 ?
Where do you initialize A(1), B(1) and P(1) ?
Why do you have a j-loop ? The k-loop suffices to implement Euler's method.
Why do you call a function named "euler" if you perform the Euler-method just below the call ?
What is A_t ?
So many questions ...
Naveen Krish
Naveen Krish 2022-3-14
nsteps = 12;
t = zeros (nsteps,1);
A = zeros (nsteps,1);
B = zeros(nsteps, 1);
P = zeros(nsteps,1);
A1 = 1;
B1 = 3;
P1 = 0;
ti=0;
tf=12*3600;
Yb = 1;
Yp = 0.15;
K = 5*10^-5;
h=3600;
timestep=[3600 1800 900 450 225];
for k = 2:13
for j=1:length(timestep)
t(k) = t(k-1)+timestep(j)
A(k) = A(k-1)+(-K*A(k-1)*B(k-1))*timestep(j);
B(k) = B(k-1)+(-Yb*(K*A(k-1)*B(k-1)))*timestep(j);
P(k)= P(k-1)+ Yp*(K*A(k-1)*B(k-1))*timestep(j);
end
end
plot (t,A)
figure (1)
plot(t,A(:,1))
plot (t,B)
figure (2)
plot(t,B(:,1))
plot(t,P)
figure (3)
plot(t,P(:,1))
% e = abs(A_t - y(1)/y(1))*100;
% plot (e)
% hold on
%end
% title('Error compared with Analytical A value');
% xlabel ('Time (t)');
% ylabel (Error(%)');
% legend ('A_e 3600','A_e 1800','A_e 900','A_e 450','A_e 225')
% end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2022-3-14
A1 = 1;
B1 = 3;
P1 = 0;
ti = 0;
tf = 12*3600;
Yb = 1;
Yp = 0.15;
K = 5*10^-5;
h = 3600;
timestep = [3600 1800 900 450 225];
for j = 1:length(timestep)
dt = timestep(j);
t{j}{1} = ti;
A{j}{1} = A1;
B{j}{1} = B1;
P{j}{1} = P1;
nsteps = tf/dt + 1;
for k = 2:nsteps
t{j}{k} = t{j}{k-1} + dt;
A{j}{k} = A{j}{k-1} + (-K*A{j}{k-1}*B{j}{k-1})*dt;
B{j}{k} = B{j}{k-1} + (-Yb*(K*A{j}{k-1}*B{j}{k-1}))*dt;
P{j}{k} = P{j}{k-1} + Yp*(K*A{j}{k-1}*B{j}{k-1})*dt;
end
end
plot(cell2mat(t{1}),cell2mat(A{1}))
hold on
plot(cell2mat(t{2}),cell2mat(A{2}))
hold on
plot(cell2mat(t{3}),cell2mat(A{3}))
hold on
plot(cell2mat(t{4}),cell2mat(A{4}))
hold on
plot(cell2mat(t{5}),cell2mat(A{5}))

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by