I got two different graphs from my code like ODE89 and Eulers method. I need to compare the graphs in it.

ODE89

clc

clear all

A0=1;

B0=3;

P0=0;

K=5*10^-5;

Yb=1;

Yp=0.15;

tspan = [0 43200];

[t,Y] = ode89(@(t,Y) odefun(t,Y, K, Yb, Yp),tspan, [A0;B0;P0]);

figure (1)

plot(t,Y(:,1))

figure (2)

plot(t,Y(:,2))

figure (3)

plot(t,Y(:,3))

function dYdt = odefun(t,Y,K,Yb,Yp)

dYdt = [(-K*Y(1)*Y(2));

(-Yb*(K*Y(1)*Y(2)));

(Yp*(K*Y(1)*Y(2)))];

end

Eulers

nsteps = 12;

t = zeros (nsteps,1);

A = zeros (nsteps,1);

B = zeros(nsteps, 1);

P = zeros(nsteps,1);

A(1) = 1;

B(1) = 3;

C(1) = 0;

K = 5*10^-5

for k = 2:13

t(k) = t(k-1)+3600

A(k) = A(k-1)+(-K*A(k-1)*B(k-1))*3600;

B(k) = B(k-1)+(-Yb*(K*A(k-1)*B(k-1)))*3600;

P(k)= P(k-1)+ Yp*(K*A(k-1)*B(k-1))*3600;

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

I need to compare the ODE89 graph of A,B and P with Euler's A,B and P

[AndresVar]

They look pretty similar, so to make it easier to draw comparisons you can plot the corresponding solutions on the same axis.

Maybe you can plot the percent differce on the right-hand side axis also.

[Naveen Krish]

how can i plot the percent difference

Answer 1

[Davide Masiello]

clc
clear all

%% ODE89
A0=1; 
B0=3; 
P0=0; 
K=5*10^-5; 
Yb=1; 
Yp=0.15;
tspan = [0 43200];
[x,Y] = ode89(@(t,Y) odefun(t,Y, K, Yb, Yp),tspan, [A0;B0;P0]);

%% Euler
nsteps = 12;
t = zeros (nsteps,1);
A = zeros (nsteps,1);
B = zeros(nsteps, 1);
P = zeros(nsteps,1);
A(1) = 1;
B(1) = 3;
C(1) = 0;
K = 5*10^-5;
for k = 2:13
    t(k) = t(k-1)+3600;
    A(k) = A(k-1)+(-K*A(k-1)*B(k-1))*3600;
    B(k) = B(k-1)+(-Yb*(K*A(k-1)*B(k-1)))*3600;
    P(k)= P(k-1)+ Yp*(K*A(k-1)*B(k-1))*3600;
end


figure(1)
subplot(1,3,1)
plot (t,A,x,Y(:,1))
xlabel('t')
ylabel('A')
subplot(1,3,2)
plot (t,B,x,Y(:,2))
xlabel('t')
ylabel('B')
subplot(1,3,3)
plot(t,C,x,Y(:,3))
xlabel('t')
ylabel('P')
legend('ODE89','Euler','Location','Best')

function dYdt = odefun(t,Y,K,Yb,Yp)
dYdt = [(-K*Y(1)*Y(2));
 (-Yb*(K*Y(1)*Y(2)));
  (Yp*(K*Y(1)*Y(2)))];
end

[Naveen Krish]

Unrecognized function or variable 'odefun'.

Error in Eluervsode89 (line 11)

[x,Y] = ode89(@(t,Y) odefun(t,Y, K, Yb, Yp),tspan, [A0;B0;P0]);

[Naveen Krish]

Nevermind got it

[Naveen Krish]

Thank you