Determining the right relative and absolute tolerances

Question

bayern science 2011-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/12853-determining-the-right-relative-and-absolute-tolerances

编辑： Torsten 2025-1-5

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?

11 个评论
显示 9更早的评论隐藏 9更早的评论

sajjad 2024-12-26

移动：John D'Errico 2024-12-26

在 MATLAB Online 中打开

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 2025-1-5

I am not receiving more output yet.

Torsten 2025-1-5

编辑：Torsten 2025-1-5

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 ?

请先登录，再进行评论。

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

Determining the right relative and absolute tolerances

11 个评论
显示 9更早的评论隐藏 9更早的评论

回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Determining the right relative and absolute tolerances

11 个评论 显示 9更早的评论隐藏 9更早的评论

回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

11 个评论
显示 9更早的评论隐藏 9更早的评论