fplot glitch when plotting a square function

Problem description

When plotting a custom function (based on the square() function) using the fplot() command, the plot leaves some gaps where none should be and when panning the plot, these gaps appear at different positions. To see this, just count the gaps (marked red) between 0 and 50 in the plots below - in the first plot, there are 4 gaps and in the second, there are 7 gaps.

I see 3 different sources of error for this:

1. A bug in my implementation
2. A hardware problem
3. A MATLAB problem

MATLAB code:

% Time-conversion factors
seconds_per_minute = 60;
minutes_per_hour   = 60;
hours_per_day      = 24;
seconds_per_hour   = seconds_per_minute*minutes_per_hour;
seconds_per_day    = seconds_per_minute*minutes_per_hour*hours_per_day;

% Custom function
photoperiod = 16;                 % [hours/day]
duty        = photoperiod/24*100; % [percent]
day_start   = 6;
I_A_max     = 150;
I_A = @(t) I_A_max/2*square(((t-day_start)*(2*pi)/24)/seconds_per_hour, duty)+I_A_max/2;

% Make the plot
fplot(@(t) I_A(t*seconds_per_day))
xlabel('time (days)')
ylabel({'Light';'[W/m^2]'})
xlim([0,50])

Can anyone comment on what is causing this behavior?

Answer 1

[Stephen23]

"Can anyone comment on what is causing this behavior?"

The main cause is missing from your list: mathematics. It essentially comes down to this:

https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

You seem to expect MATLAB to take infinite samples within that range. Such a thing is not possible.

But you could specify the number of samples (following the nyquist-shannon sampling theorem, at more than double your data frequency):

% Time-conversion factors
seconds_per_minute = 60;
minutes_per_hour   = 60;
hours_per_day      = 24;
seconds_per_hour   = seconds_per_minute*minutes_per_hour;
seconds_per_day    = seconds_per_minute*minutes_per_hour*hours_per_day;
% Custom function
photoperiod = 16;                 % [hours/day]
duty        = photoperiod/24*100; % [percent]
day_start   = 6;
I_A_max     = 150;
I_A = @(t) I_A_max/2*square(((t-day_start)*(2*pi)/24)/seconds_per_hour, duty)+I_A_max/2;
% Make the plot
md = 50*16; % a wild guess (you need to specify this)
fplot(@(t) I_A(t*seconds_per_day), 'MeshDensity',md)
xlim([0,50])