How to make a FOR loop to be efficient

Hi,

I have a code that looks something like this.

tol = 1e-16;
a=1; b=500*pi; % b can be larger
pts = a:pi:b;
err=realmax;
fx = @(x) besselj(0,x); % or some other function
for i=1:size(pts)-1
   [result,err] = somefun(fx,pts(i),pts(i+1)); %or some other function with the sub-intervals as scalar inputs.
end %  for loop

Problem is that this type of logic takes a long time to converge. Is there a way to vectorize this loop? Modifying somefun.m itself is not really an option since they could be built-in functions.

Thanks, Jung.

[Nirmal]

what is pi? is it calculated by the function quad or is it predetermined?

[Jung]

pi = 3.1416....

Answer 1

[Sean de Wolski]

eps(pi)

ans = 4.4409e-16

1e-16 is going to be a hard target to hit... Why not use a realistic tolerance?

[Geoff]

Hmmm let me count the number of times over the past 20 years that I've needed a double-precision accuracy of one in ten quadrillion........... I would say it's approximately floor(1e-16).

Answer 2

[Mike Hosea]

There is no way to do what you ask, and it would not help you much if there were. What is it that that you're really trying to do?

Answer 3

[Jung]

Hi Mike,

My ultimate goal is to run a loop more efficiently - i.e. less number of iterations. The code is actually overly simplified since the original code is a bit too long. But essentially, I need to run quad (or some other function with the sub-interval endpoints as scalars) for each of the interval in 1:pi:500*pi (interval could be shorter or significantly longer). I've edited the original posting above to hopefully make it clearer.

Thanks, Jung

[Sean de Wolski]

Well then make the tolerances realistic so that you stand a chance of hitting them after a finite number of iterations. For example: 10^-8

Also, make sure it's in a function so the JIT will pay attention to it.

[Jung]

problem with reducing the tolerance is that this is part of a larger application in physics labs. so the loop in itself will be run many 10s of thousands of times (or whatever it takes to get to a stable solution), so the final final solution may have rounding errors of much greater magnitude than 1e-8. and besides physicists would just laugh in my face if i told them my tolerance was that high.

[Sean de Wolski]

working in the double precision floating point world you have to have realistic expectations for this kind fo thing. Looking for zero error (which is what you're doing above) is improbable if not impossible to attain. If you need that kind of precision to appease your physics buddies, use the Sympolic Math Toolbox and VPA along with INT or whatever other operations you need to hit that kind of precision.

http://scisyhp.deviantart.com/art/Physicists-vs-Engineers-179857993

[Mike Hosea]

The nature of "somefun" is critical. If it's numerical integration and the integrand is smooth, you could use a high-order Gaussian quadrature rule, evaluate the function at each node with column j corresponding to pts(j) <= x <= pts(j+1), and do the bulk of the work in one matrix multiplication. This is obviously not an adaptive scheme, but there is a symbolic toolbox function that will calculate Gauss-Legendre nodes and weights, so you can go as high in order as you want (provided you have the memory for the resulting matrix of function evaluations). There are also various strategies you could employ to estimate the error.

Answer 4

[newbie traveller]

if you have GPU (graphics Processing Unit ) then you can use them to compute the calculation in the for loop....check that whether the work that you are doing in the for loop is independent data or not!!i mean independent means one data is not dependent on the calculation of other data.....check that or you can also use the cores of the CPU or labs ......you can see on the command "open matlabpool" to get idea from this this actually helps you to use all of your cores or lab in the CPU........so the computation time will be reduced to a great extent and remember the command for PCT(parallel computing toolbox) the for will be changed and you will have to use 'parfor'..........clear enough?

[Jung]

we do have parallel toolbox and have tried parfor. but problem here is that not all machines have the toolbox so we can't distribute it to others. i've also tried building a stand'alone .exe, but people need to tinker with the code so that's a no-go either.

Answer 5

[Ryan]

Perhaps break this into two separate iterations where the first optimization happens more crudely using parameters spaced more greatly apart, then perform another optimization using more closely spaced parameter values based upon the results of the first loop.