why am I getting this difference in the plotting?

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

[Sam Chak]

Hi @Desiree

Something in the equation is not being conserved. Please provide the MATLAB code for further investigation. Why? Because the specified initial value but the graph for begins from 0.1.

[Desiree]

These are the plots for and not for .

[Desiree]

@Sam Chak Here are the codes

Answer 1

[Torsten]

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

[Desiree]

I took your code and changed it a bit for for example:

a = 3;
b = 0.95;
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));
%plot(Z1(:,2))
fun = @(t,z) -(a*z(1)-b)*(1-z(1)-(1-b)/a*log((a*x0-b)/abs(a*z(1)-b)));
tspan = [0,500];
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)));
plot(T1,Z1(:,2),T2,Z3)
xlim([0,500])

[Desiree]

But gives an unpleasant picture... why does this happen?

[Torsten]

The solution from the second method for x can only be used as long as the signs of a*x0-b and a*x-b agree. In this case: as long as x>b/a. For this period of time, the solutions of the two methods agree as you can see from the plot.

For greater values of t, the solution for x becomes complex-valued and can no longer be used to recover y.

a = 3;
b = 0.95;
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);
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));
idx = find(Z2<=b/a,1,'first')

idx = 1387

T2=T2(1:idx-1);
Z2=Z2(1:idx-1);
Z2 = 1-Z2-(1-b)/a*log((a*x0-b)./(a*Z2-b));
plot(T1,Z1(:,2),T2,Z2)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

xlim([0,500])

[Desiree]

@Torsten I'm aware of the possibility of imaginary parts which is why in my rewriting of the code, I put absolute value in the log function yielding to the plot I showed earlier today. But that plot doesn't make sense and don't get why it behaves like this...

[Torsten]

I repeat:

The solution of the system

x'=-a*x*y+b*y,  x(0) = x0
y'=a*x*y-y,     y(0) = y0

and of the system

x' = -(a*x-b)*(1-x-(1-b)/a*log((a*x0-b)/(a*x-b))), x(0) = x0
y  = 1-x-(1-b)/a*log((a*x0-b)/(a*x-b))

are only equivalent as long as (a*x0-b)/(a*x-b) > 0.

You cannot simply replace (a*x0-b)/(a*x-b) by its absolute value.

[Desiree]

@Torsten Since there are abvously differences in T1 and T2, by using knnsearch, I plotted the difference between the solutions of for the two methods and the difference between for the two methods. Why the difference for becomes much smaller than for ? Regarding I suspect it is because ode45 isn't designed to preserve quantities but still can't explain it to myself the real reason.

[Torsten]

T1, T2 and knnsearch have not yet been used. So please post the code you are referring to.

[Desiree]

@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]

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