Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Array Parsing

Subject: Array Parsing

From: Scott

Date: 20 Jun, 2013 02:41:07

Message: 1 of 7

I have an array of increasing numbers that represent points in time. I then create a separate "clock" process that ticks along and at each point checks the array to determine if the next desired point in time has been reached. If it has the clock time is recorded. I then print out both the original array and the identified clock points. I was expecting the two outputs to be identical but that is not the case. Would someone mind having a look at my code please (see below)? The example is quite artificial because I have boiled the issue down to its essentials (I am trying to use this method to model jump-diffusion processes).

Thanks,
Scott


T = 1.0; %final time
steps = 20; %numberof time steps for clock
dt = T/steps; %clock time step
N = 10;
times = T/N : T/N : T %create time points
output = zeros(1,N);
t = 0;
n = 1;
for i=1:steps %clock process
    t = t + dt;
    if times(n) <= t %check to see if we have reached jump time
        output(1,n) = t; %record clock time
        n = n + 1;
    end
end
output %print identified clock times

Subject: Array Parsing

From: Josh Meyer

Date: 20 Jun, 2013 13:13:31

Message: 2 of 7

I was expecting the two outputs to be identical but that is not the case.
Would someone mind having a look at my code please (see below)? The example
is quite artificial because I have boiled the issue down to its essentials
(I am trying to use this method to model jump-diffusion processes).
>
> T = 1.0; %final time
> steps = 20; %numberof time steps for clock
> dt = T/steps; %clock time step
> N = 10;
> times = T/N : T/N : T %create time points
> output = zeros(1,N);
> t = 0;
> n = 1;
> for i=1:steps %clock process
> t = t + dt;
> if times(n) <= t %check to see if we have reached jump time
> output(1,n) = t; %record clock time
> n = n + 1;
> end
> end
> output %print identified clock times

The problem is that you are accumulating small errors that are causing the
"if times(n) <= t" statement to evaluate to false when you want it to be
true. They arise because of the line "t = t +dt", where dt = 0.05.

For example, in the iteration corresponding to i=6 and n=3, this happens:

 t = t+0.05

t =

    0.3000

times(3) <= t

ans =

     0

The reason is because:

times(3) - t

ans =

   5.5511e-17

To correct this problem, compare the difference of times(n) and t to a
tolerance instead of to each other:

T = 1.0; %final time
steps = 20; %numberof time steps for clock
dt = T/steps; %clock time step
N = 10;
times = T/N : T/N : T %create time points
output = zeros(1,N);
t = 0;
n = 1;
tol = eps;
for i=1:steps %clock process
    t = t + dt;
    if abs(times(n)-t) <= tol %check to see if we have reached jump time
        output(1,n) = t; %record clock time
        n = n + 1;
    end
end
output %print identified clock times
output %print identified clock times

Subject: Array Parsing

From: dpb

Date: 20 Jun, 2013 13:37:46

Message: 3 of 7

On 6/19/2013 9:41 PM, Scott wrote:
> I have an array of increasing numbers that represent points in time. I
> then create a separate "clock" process that ticks along and at each
> point checks the array to determine if the next desired point in time
> has been reached. If it has the clock time is recorded. I then print out
> both the original array and the identified clock points. I was expecting
> the two outputs to be identical but that is not the case. ...
...
>
> T = 1.0; %final time
> steps = 20; %numberof time steps for clock
> dt = T/steps; %clock time step
> N = 10;
> times = T/N : T/N : T %create time points
> output = zeros(1,N);
> t = 0;
> n = 1;
> for i=1:steps %clock process
> t = t + dt;
> if times(n) <= t %check to see if we have reached jump time

...

As another poster noted the above is fraught w/ floating point rounding
errors. Another solution would be to use the integer step numbers as
comparison instead of the actual estimated times.

Another possible way would be to use a datenum and use integer steps to
create them from an input vector. I'm presuming from the above that the
timesteps are fairly coarse. If one uses a consistent technique of 1:N
(say msecs) in the datenum argument, then the internals of datenum will
generate consistent ML datenums and the comparison of them will work
because the rounding is the same having used the same generation method.

--

Subject: Array Parsing

From: Steven_Lord

Date: 20 Jun, 2013 13:48:42

Message: 4 of 7



"Scott " <itsme@somewhere.com> wrote in message
news:kptq43$g96$1@newscl01ah.mathworks.com...
> I have an array of increasing numbers that represent points in time. I
> then create a separate "clock" process that ticks along and at each point
> checks the array to determine if the next desired point in time has been
> reached. If it has the clock time is recorded. I then print out both the
> original array and the identified clock points. I was expecting the two
> outputs to be identical but that is not the case. Would someone mind
> having a look at my code please (see below)? The example is quite
> artificial because I have boiled the issue down to its essentials (I am
> trying to use this method to model jump-diffusion processes).

See question 1 in the Math/Algorithms section of the newsgroup FAQ.

http://matlab.wikia.com/wiki/FAQ

Quoting from the Cleve's Corner article linked to by that question:

"So, the value stored in t is very close to, but not exactly equal to, 0.1.
The distinction is occasionally important. For example, the quantity 0.3/0.1
is not exactly equal to 3 because the actual numerator is a little less than
0.3 and the actual denominator is a little greater than 0.1.

Ten steps of length t are not precisely the same as one step of length 1.
MATLAB is careful to arrange that the last element of the vector 0:0.1:1 is
exactly equal to 1, but if you form this vector yourself by repeated
additions of 0.1, you will miss hitting the final 1 exactly"

I've modified your example slightly to add some extra display:

T = 1.0; %final time
steps = 20; %numberof time steps for clock
dt = T/steps; %clock time step
N = 10;
times = T/N : T/N : T; %create time points
output = zeros(1,N);
t = 0;
n = 1;
fprintf('n\ttimes\t\t\t\tt\t\t\t\t\tt-times\n');
for i=1:steps %clock process
    t = t + dt;
    fprintf('%d\t%bx\t%bx\t%g\n', n, times(n), t, t-times(n))
    if times(n) <= t %check to see if we have reached jump time
        output(1,n) = t; %record clock time
        n = n + 1;
    end
end
output %print identified clock times

You will see that at the steps where the output results are different from
what you expect, t is smaller than times by a very small amount. If you look
at the hex representations of times and t you'll see that the difference is
in the least significant bit. That's still enough to fail your "check to see
if we have reached jump time" condition.

Since this is an artificial example I can't offer suggestions on how to
resolve the problem in your original code. In this example I would rescale
the problem; change the units of T, dt, times, etc. from seconds to
(seconds/20). This would make times:

2*T:2*T:20*T % All integer values

and dt would now be 1. This avoids roundoff error entirely, since integers
in the range with which you're working are exactly representable in double
precision.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Array Parsing

From: Scott

Date: 21 Jun, 2013 03:03:09

Message: 5 of 7

"Steven_Lord" <slord@mathworks.com> wrote in message <kpv17p$ef3$1@newscl01ah.mathworks.com>...
>
>
> "Scott " <itsme@somewhere.com> wrote in message
> news:kptq43$g96$1@newscl01ah.mathworks.com...
> > I have an array of increasing numbers that represent points in time. I
> > then create a separate "clock" process that ticks along and at each point
> > checks the array to determine if the next desired point in time has been
> > reached. If it has the clock time is recorded. I then print out both the
> > original array and the identified clock points. I was expecting the two
> > outputs to be identical but that is not the case. Would someone mind
> > having a look at my code please (see below)? The example is quite
> > artificial because I have boiled the issue down to its essentials (I am
> > trying to use this method to model jump-diffusion processes).
>
> See question 1 in the Math/Algorithms section of the newsgroup FAQ.
>
> http://matlab.wikia.com/wiki/FAQ
>
> Quoting from the Cleve's Corner article linked to by that question:
>
> "So, the value stored in t is very close to, but not exactly equal to, 0.1.
> The distinction is occasionally important. For example, the quantity 0.3/0.1
> is not exactly equal to 3 because the actual numerator is a little less than
> 0.3 and the actual denominator is a little greater than 0.1.
>
> Ten steps of length t are not precisely the same as one step of length 1.
> MATLAB is careful to arrange that the last element of the vector 0:0.1:1 is
> exactly equal to 1, but if you form this vector yourself by repeated
> additions of 0.1, you will miss hitting the final 1 exactly"
>
> I've modified your example slightly to add some extra display:
>
> T = 1.0; %final time
> steps = 20; %numberof time steps for clock
> dt = T/steps; %clock time step
> N = 10;
> times = T/N : T/N : T; %create time points
> output = zeros(1,N);
> t = 0;
> n = 1;
> fprintf('n\ttimes\t\t\t\tt\t\t\t\t\tt-times\n');
> for i=1:steps %clock process
> t = t + dt;
> fprintf('%d\t%bx\t%bx\t%g\n', n, times(n), t, t-times(n))
> if times(n) <= t %check to see if we have reached jump time
> output(1,n) = t; %record clock time
> n = n + 1;
> end
> end
> output %print identified clock times
>
> You will see that at the steps where the output results are different from
> what you expect, t is smaller than times by a very small amount. If you look
> at the hex representations of times and t you'll see that the difference is
> in the least significant bit. That's still enough to fail your "check to see
> if we have reached jump time" condition.
>
> Since this is an artificial example I can't offer suggestions on how to
> resolve the problem in your original code. In this example I would rescale
> the problem; change the units of T, dt, times, etc. from seconds to
> (seconds/20). This would make times:
>
> 2*T:2*T:20*T % All integer values
>
> and dt would now be 1. This avoids roundoff error entirely, since integers
> in the range with which you're working are exactly representable in double
> precision.
>
> --
> Steve Lord
> slord@mathworks.com
> To contact Technical Support use the Contact Us link on
> http://www.mathworks.com

Thanks for your help Steve,
Scott

Subject: Array Parsing

From: Scott

Date: 21 Jun, 2013 03:43:08

Message: 6 of 7

Thanks to all who took the time to answer my question. Your comments were very helpful.

Scott

Subject: Array Parsing

From: Scott

Date: 21 Jun, 2013 04:38:07

Message: 7 of 7

I thought I would post the solution to this issue in case anyone else reading this thread is interested. Motivated by one of the respondents, I converted the array "times" to integer points via "int32(steps*times)" and then compared the looping variable "i" to this converted array. The code now reads:

T = 1.0; %final time
steps = 836492; %numberof time steps for clock
dt = T/steps; %clock time step
N = 10;
times = T/N : T/N : T %create time points
times = int32(steps*times); %convert to integers

output = zeros(1,N);
t = 0;
n = 1;
for i=1:steps %clock process
    t = t + dt;
    if times(n) <= i %check to see if we have reached jump time
        output(1,n) = t; %record clock time
        n = n + 1;
    end
end
output %print identified clock times

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us