Determining the right relative and absolute tolerances

My program generates an essentially sinusoidal output. Using the various relative and absolute tolerances I obtained the following fits to a sine function (with about 2e4 data points):

ReTol = 1e-5, AbTol = 1e-7, Sinefit amplitude = 22.1383, phase = -0.10532, fminsearch err = 71211.296;

ReTol = 1e-5, AbTol = 1e-8, Sinefit amplitude = 22.1383, phase = -0.10532, fminsearch err = 71211.296;

ReTol = 1e-6, AbTol = 1e-7, Sinefit amplitude = 24.5425, phase = +0.22298, fminsearch err = 61670.879;

ReTol = 1e-6, AbTol = 1e-8, Sinefit amplitude = 24.5335, phase = +0.22210, fminsearch err = 106022.2668;

I tend to choose ReTol = 1e-5, AbTol = 1e-7, as reducing absolute tolerance doesn't change the result. However if I reduce relative tolerance I obtain different results (especially the phase). Thus what should I choose? ReTol = 1e-5, AbTol = 1e-7, or ReTol = 1e-6, AbTol = 1e-7?

[sajjad]

How to get more outputs of chaotic maps by adjusting Relative Tolerance (ReTo1) Absolute Tolerance (AbsT1)? I am running a code i need more points for my futher work. Please guide. Thank you.

[Walter Roberson]

You should ask this as a new question.

It would be a good idea to share the actual code instead of just an image of the code.

[sajjad]

This is first file of code:

SAVE THIS PART WITH NAME funn

function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);
%persistent iter
%if isempty(iter)
%  iter = 0;
%end
%iter = iter + 1;
%if mod(iter,1000)==0
%  t
%  iter = 0;
%end
G = y(1:nx);
F = y(nx+1:2*nx);
% F
% Compute spatial derivatives
dFdx = zeros(nx,1);
d2Fdx2 = zeros(nx,1);
dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);
dFdx(nx) = Hdot(t);           % This corresponds to F'(1,t) = H'(t)
d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;
%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;
Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);
d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;
% Compute temporal derivatives
dFdt = zeros(nx,1);
dFdt(1) = F(1);
dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...
    F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);
dFdt(nx) = F(nx) + Hdot(t);     % This corresponds to F(1,t) =- H'(t)
% G
% Compute spatial derivatives
dGdx = zeros(nx,1);
d2Gdx2 = zeros(nx,1);
dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);
d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;
% Compute temporal derivatives
dGdt = zeros(nx,1);
dGdt(1) = G(1);
dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...
    1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...
    1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...
    2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...
    (d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...
    (H(t)^2*Reynolds);
dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);
%Taken from the article, should be the same as set
%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));
dydt = [dGdt;dFdt];
end

Following is the second part of code:

format long;
clear all;
close all;
clc;
nx = 101;        %The number of points or grid cells used to discretize the spetial domain
xstart = 0.0;
xend = 1.0;
x = linspace(xstart,xend,nx).';
dx = (xend-xstart)/(nx-1);
tstart = 0;
tend =8000;
tspan = [tstart,tend];
G0 = zeros(nx,1);        % corresponding to F(etta,0)=0
F0 = zeros(nx,1);        % corresponding to G(etta,0)=0
y0 = [G0;F0];
delta = 0.3;
Reynolds = 1000;
H = @(t) 1 + delta*cos(2*t);
Hdot = @(t) -delta*2*sin(2*t);
M = [eye(nx),zeros(nx);zeros(nx,2*nx)];
M(1,1) = 0;
M(nx,nx) = 0;
options = odeset('Mass',M, 'RelTol', 1e-7, 'AbsTol', 1e-8);
[T,Y] = ode23t(@(t,y)funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot),tspan,y0, options);
G = Y(:,1:nx);
F = Y(:,nx+1:2*nx);
idx = T >= 0;
figure(1)
plot(T(idx),G(:,nx));

% Assuming G and T have already been defined and filtered for t <= 800
G1 = G (: , nx); % Values of G at eta = 1
t_range = T (idx); % Corresponding time points

% Find maxima
[maxima, locs_max] = findpeaks (G1 (idx)); % Find peaks and their locations
max_times = t_range (locs_max); % Times at which maxima occur
% Find minima
[minima, locs_min] = findpeaks (-G1 (idx)); % Negate to find minima
min_times = t_range (locs_min); % Times at which minima occur


% Plot Maxima as points only
figure(2);
plot(max_times, maxima, 'ro', 'MarkerFaceColor', 'r','MarkerSize', 3)  % Plot maxima as red points
xlabel('Time (t)');
ylabel('Maxima of G(1,t)');
title('Maxima of G(1,t)');
grid on;


% Plot Minima as points only
figure(3);
plot (min_times, -minima, 'go', 'MarkerFaceColor', 'g','MarkerSize', 3) % Plot minima as green points
xlabel (' Time (t)');
ylabel (' Minima of G (1, t)');
title (' Minima of G (1, t)');
grid on;

[pks1,locs] = findpeaks(G(:,nx));
figure(4)
plot(T(locs),pks1,'o')

[pks2,locs] = findpeaks(-G(:,nx));
%pks = -pks
figure(5)
plot(T(locs),pks2,'o')


% Create return map for maxima: plot max(i+1) vs max(i)
figure (6)
plot(maxima(1:end-1), maxima(2:end), 'ro', 'MarkerFaceColor', 'r','MarkerSize', 3);  % Return map for maxima
xlabel('Max(i)');
ylabel('Max(i+1)');
title('Return Map of Maxima');
grid on;

% Display data for return map
return_map_maxima = [maxima(1:end-1), maxima(2:end)];  % Prepare data for return map
disp('Return Map Data for Maxima (Max(i), Max(i+1)):');
disp(return_map_maxima);  % Display the consecutive maxima pairs in the Command Window


% Create return map for minima: plot min(i+1) vs min(i)
figure(7)
plot(-minima(1:end-1), -minima(2:end), 'go', 'MarkerFaceColor', 'g','MarkerSize', 3);  % Return map for minima
xlabel('Min(i)');
ylabel('Min(i+1)');
title('Return Map of Minima');
grid on;

% Display data for return map
return_map_minima = [minima(1:end-1), minima(2:end)];  % Prepare data for return map
disp('Return Map Data for Maxima (Max(i), Max(i+1)):');
disp(return_map_minima);  % Display the consecutive maxima pairs in the Command Window

[sajjad]

in above code i ma getting 2000 plus maximas i want to in crease these maximas. How to increase?

[Torsten]

I don't know exactly what you mean, but "tspan" is the vector of times when you want to output the solution.

You only define start and end in your code. This means that the solution is returned at times determined by the solver.

If you set "tspan" as a vector with more than two elements, the solution will be supplied exactly at these specified times.

[Walter Roberson]

... This should all have gone as a new Question.

[Torsten]

This should all have gone as a new Question.

Yes, best would have been here:

https://uk.mathworks.com/matlabcentral/answers/2154545-how-to-solve-the-system-of-time-dependent-coupled-pde-s?s_tid=srchtitle

[sajjad]

How to adjust the relative and absolute tolerances in ODE23t for more outputs.?

Actully i need 600000 outpouts of maxima.

How to acheive?

at RelTol', 1e-3, 'AbsTol', 1e-6, i got the 2000 plus points.

Or i have to descritize tspan = linspace(tstart,tend,10000); like this??

[Torsten]

Or i have to descritize tspan = linspace(tstart,tend,10000); like this??

Yes.

If the results are stable even if you strengthen the tolerances, you don't need to change them for receiving more outputs. Just redefine the vector of output times in "tspan" as you suggested.

The test whether the results remain stable even if you choose smaller tolerances should always be done before trusting a solution - especially for chaotic systems. I had done this for your system, but couldn't find an acceptable tolerance level - the results didn't stabilize. Maybe you were more successful.

[sajjad]

I am not receiving more output yet.

[Torsten]

As I answered in another question of yours

https://uk.mathworks.com/matlabcentral/answers/491227-ode45-solver-results-changes-how-to-select-optimal-reltol-and-abstol?s_tid=srchtitle

, it's impossible to get more maxima because the number of maxima if given by the solution of the PDE system. E.g. if your solution were y(x)=sin(x), you would get maxima and minima at odd multiples of pi/2, and no modification of the code could give you more of them. Or do I get wrong what you are asking for ?