Unrecognized function or variable 'A'. Error in ode89vsrk (line 37) plot (t,A,x,Y(:,1))
2 views (last 30 days)
Show older comments
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]);
%%RK Method
h=3600;
A0=1;
B0=3;
P0=0;
K=5*10^-5;
Yb=1;
Yp=0.15;
t = 0:h:43200;
fyt = @(t,y) [(-K*y(1)*y(2));
(-Yb*(K*y(1)*y(2)));
(Yp*(K*y(1)*y(2)))];
Y = zeros(3,numel(t))
Y(1,1) = 1.0;
Y(2,1) = 3.0;
Y(3,1) = 0;
for i=1 : numel(t)-1
k1 = fyt(t(i),Y(:,i));
k2 = fyt(t(i)+0.5*h,Y(:,i)+0.5*h*k1);
k3 = fyt(t(i)+0.5*h,Y(:,i)+0.5*h*k2);
k4 = fyt(t(i)+h,Y(:,i)+h*k3);
Y(:,i+1) = Y(:,i) + (h/6)*(k1+2*k2+2*k3+k4);
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','RK Method','Location','Best')
I need to plot ODE89 vs RK Method
0 Comments
Accepted Answer
Torsten
on 17 Mar 2022
Edited: Torsten
on 17 Mar 2022
%% ODE89
A0=1;
B0=3;
P0=0;
K=5*10^-5;
Yb=1;
Yp=0.15;
xM = [0 43200];
fyt = @(t,y,K,Yb,Yp) [(-K*y(1)*y(2));
(-Yb*(K*y(1)*y(2)));
(Yp*(K*y(1)*y(2)))];
[xM,yM] = ode89(@(t,y) fyt(t,y, K, Yb, Yp),xM, [A0;B0;P0]);
%%RK Method
h=3600;
A0=1;
B0=3;
P0=0;
K=5*10^-5;
Yb=1;
Yp=0.15;
xR = (0:h:43200).';
fyt = @(t,y) [(-K*y(1)*y(2)),...
(-Yb*(K*y(1)*y(2))),...
(Yp*(K*y(1)*y(2)))];
yR = zeros(numel(xR),3)
yR(1,1) = A0;
yR(1,2) = B0;
yR(1,3) = P0;
for i=1 : numel(xR)-1
k1 = fyt(xR(i),yR(i,:));
k2 = fyt(xR(i)+0.5*h,yR(i,:)+0.5*h*k1);
k3 = fyt(xR(i)+0.5*h,yR(i,:)+0.5*h*k2);
k4 = fyt(xR(i)+h,yR(i,:)+h*k3);
yR(i+1,:) = yR(i,:) + (h/6)*(k1+2*k2+2*k3+k4);
end
figure(1)
subplot(1,3,1)
plot (xR,yR(:,1),xM,yM(:,1))
xlabel('t')
ylabel('A')
subplot(1,3,2)
plot (xR,yR(:,2),xM,yM(:,2))
xlabel('t')
ylabel('B')
subplot(1,3,3)
plot(xR,yR(:,3),xM,yM(:,3))
xlabel('t')
ylabel('P')
legend('ODE89','RK Method','Location','Best')
More Answers (0)
See Also
Categories
Find more on Calendar in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!