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

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

[Jan]

What is your question?

[Naveen Krish]

using Euler’s Method instead of ODE89. Use a timestep of 3600 sec. Plot the concentrations of A, B, and P on the same plot generated from ODE89 to compare the results.

[Naveen Krish]

I already have my ODE89 results i just wanna compare it in the same graph with euler's method

[Torsten]

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]

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

Answer 1

[Torsten]

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}))

[Naveen Krish]

t{j}{1} = ti; Unable to perform assignment because brace indexing is not supported for variables of this type.