why am I getting this difference in the plotting?
92 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Desiree
le 14 Juil 2024 à 7:42
Modifié(e) : Torsten
le 19 Juil 2024 à 13:01
I have a system of equations:
This system satisfies a relation
where is an initial condition.
For I plotted the relation in two ways:
First way by solving the system numerically using ode45 with RelTol 1e-12 and AbsTol 1e-15 and by having the solution for x I plotted y.
Second way by solving the equation again with ode45 with Reltol 1e-9 and Abstol 1e-12 and then defining y=... and plotting it.
However, the plots are very different and I don't understand the reason why:
Here is the plot by first method:
Here is the plot by second method:
Help is appreciated!
Codes:
%--------------------------------- second method
function[Y] = S_with_I_defined(a,b,x0)
% a function to define for ode45
d=abs(a*x0-b);
function dS = SIpS1_pR(t,y)
dS = -(a*y-b)*(1-y-((1-b)/a)*log(d/abs(a*y-b)));
end
% solving the system and sketching the curves S,I,R
options = odeset('Refine',6,'RelTol',1e-9,'AbsTol',1e-12);
[t,y] = ode45(@SIpS1_pR, [0 1500], x0, options);
Y=y;
end
%-------------------------------------------- first method
function[X,Y] = RK_SI_pS_1_minus_pR7(a,b,x0,y0)
% a function to define for ode45
function dy = SIpS1_pR(t,y)
dy = zeros(2,1);
dy(1) = - a*y(1)*y(2)+b*y(2);
dy(2) = a*y(1)*y(2) -y(2);
end
% solving the system and sketching the curves
options = odeset('Refine',1,'RelTol',1e-12,'AbsTol',1e-18);
[t,y] = ode45(@SIpS1_pR, [0 1500], [x0 y0], options);
X=y(:,1);
Y=y(:,2);
end
3 commentaires
Réponse acceptée
Torsten
le 14 Juil 2024 à 8:48
Modifié(e) : Torsten
le 14 Juil 2024 à 9:43
Your choice of a and b must be different from the values you posted:
% First method
a = 3;
b = 0.9;
x0 = 0.9;
y0 = 1-x0;
fun = @(t,z)[-a*z(1)*z(2)+b*z(2);a*z(1)*z(2)-z(2)];
tspan = [0,1500];
z0 = [x0;y0];
[T1,Z1] = ode45(fun,tspan,z0);
% Second method
fun = @(t,z) -(a*z(1)-b)*(1-z(1)-(1-b)/a*log((a*x0-b)/(a*z(1)-b)));
tspan = [0,500];
z0 = x0;
[T2,Z2] = ode45(fun,tspan,z0,odeset('RelTol',1e-12,'AbsTol',1e-12));
Z2 = 1-Z2-(1-b)/a*log(abs((a*x0-b)./(a*Z2-b)));
plot(T1,Z1(:,2),T2,Z2)
xlim([0,100])
9 commentaires
Torsten
le 19 Juil 2024 à 12:59
Modifié(e) : Torsten
le 19 Juil 2024 à 13:01
Why the difference for x(t) becomes much smaller than for y(t) ? Regarding y(t) I suspect it is because ode45 isn't designed to preserve quantities but still can't explain it to myself the real reason.
I don't know, but a difference in the results in the order of 1e-5 for y is not that bad. Maybe it's because the second method uses ode45 only to solve for x - so you don't have a control over the error in y.
a = 3;
b = 0.9;
x0 = 0.9;
y0 = 1-x0;
fun = @(t,z)[-a*z(1)*z(2)+b*z(2);a*z(1)*z(2)-z(2)];
tspan = 0:1500;
z0 = [x0;y0];
[T1,Z1] = ode45(fun,tspan,z0,odeset('RelTol',1e-12,'AbsTol',1e-18));
Z1_tilde = 1-Z1(:,1)-(1-b)/a*log((a*x0-b)./(a.*Z1(:,1)-b));
figure
plot(Z1_tilde-Z1(:,2))
title('Difference $\hat{Y}(x(t))-Y(t)$', 'Interpreter','latex')
xlabel('t')
%-----------------------------------------
fun = @(t,z) -(a*z(1)-b)*(1-z(1)-(1-b)/a*log((a*x0-b)/(a*z(1)-b)));
tspan = 0:1500;
z0 = x0;
[T2,Z2] = ode45(fun,tspan,z0,odeset('RelTol',1e-12,'AbsTol',1e-18));
Z3 = 1-Z2-(1-b)/a*log((a*x0-b)./(a*Z2-b));
figure
plot(T1,Z1(:,1)-Z2,'LineWidth',2)
title('Difference of X(t) for the two methods')
xlabel('t')
figure
plot(T1,Z1(:,2)-Z3,'LineWidth',2)
title('Difference of Y(t) for the two methods')
xlabel('t')
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!