Unrecognized function or variable 'A'. Error in ode89vsrk (line 37) plot (t,A,x,Y(:,1))
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Naveen Krish
el 17 de Mzo. de 2022
Comentada: Naveen Krish
el 18 de Mzo. de 2022
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 comentarios
Respuesta aceptada
Torsten
el 17 de Mzo. de 2022
Editada: Torsten
el 17 de Mzo. de 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')
Más respuestas (0)
Ver también
Categorías
Más información sobre Calendar en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!