why am I getting this difference in the plotting?

3 次查看（过去 30 天）
Desiree
Desiree 2024-7-14
编辑： Torsten 2024-7-19
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 个评论
显示 1更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-7-14
编辑：Torsten 2024-7-14
Ran in:
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 个评论
显示 7更早的评论
Desiree
Desiree 2024-7-19
编辑：Desiree 2024-7-19
@Torsten So the full code looks like this:
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)./abs(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)/abs(a*z(1)-b)));
tspan = [0,1500];
z0 = x0;
[T2,Z2] = ode45(fun,tspan,z0,odeset('RelTol',1e-12,'AbsTol',1e-14));
Z3 = 1-Z2-(1-b)/a*log(abs((a*x0-b)./(a*Z2-b)));
figure
plot(T1,Z1(:,2),T2,Z3)
xlim([0,500])
% figure
% plot(Z3-Z1(1:length(Z3),2))
[idx, ~]=knnsearch(T1,T2);
TT1=zeros(length(idx),1);
Z1_new_1=zeros(length(idx),1);
Z1_new_2=zeros(length(idx),1);
for i=1:length(idx)
TT1(i)=T1(idx(i));
Z1_new_1(i)=Z1(idx(i),1);
Z1_new_2(i)=Z1(idx(i),2);
end
figure
plot(TT1,Z1_new_2-Z3,'LineWidth',2)
title('Difference of Y(t) for the two methods')
xlabel('t')
figure
plot(TT1,Z1_new_1-Z2)
title('Difference of X(t) for the two methods')
xlabel('t')
Torsten
Torsten 2024-7-19
编辑：Torsten 2024-7-19
Ran in:
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')

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by