getting rid of for loop in a function with long arrays
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Here is a small function to calculate the neff of SM fiber. To this function I insert three long arrays: n_core, n_clad and lamda (each can be even 2^17 long array - will need it for FFT later on) If the arrays are very long, the computer get stuck and even if not it's a long processing time. Was wondering if someone knows how to make this more efficient since matlab is more efficient on matrix calculations. Basically get rid of the for loop. Note that this loop exists to solve eq. numerically something I wasn't manage to do with the "solve" function.
Thanks.
Here is the code:
function [neff,betta]=neff_fiber2(n_core,n_clad,core_diameter,lamda)
%a is the radius in um;
%lamda is the wavelenght in um;
a = core_diameter/2; %radius
k0 = 2*pi./lamda;
V = k0.*a.*sqrt(n_core.^2 - n_clad.^2);%V number
betta = zeros(1,length(V));
%%now we want to solve the following
% J0(x) / ( x*J1(x))=K0(sqrt(V^2-x^2) )/(sqrt(V^2-x^2)*K1(sqrt(V^2-x^2)))
% where x=Ka
for ii=1:length(V)
x=linspace(0,V(ii),1e6);
Diff=abs(besselj(0,x)./(x.*besselj(1,x))- ...
besselk(0,sqrt(V(ii).^2-x.^2))./(sqrt(V(ii)^2- ...
x.^2).*besselk(1,sqrt(V(ii)^2-x.^2))));
[~, ind1] = min(Diff(1,:)); % using the index of the minimum value in the row
Kappa = x(ind1)./a;
betta(ii) = sqrt(k0(ii).^2 * n_core(ii)^2 - Kappa.^2);
end
neff=betta./k0;
end
13 comentarios
dpb
el 5 de Ag. de 2017
Yet again, what are some typical input values???
Give us something that actually will run.
Respuesta aceptada
dpb
el 5 de Ag. de 2017
Editada: dpb
el 6 de Ag. de 2017
OK, once have some data to work with it's not too bad at all...it is virtually always really, really helpful to be able to have actual problem to see; particularly with nonlinear functions.
Try the following
function [neff,betta]=neff_fiber2(n_core,n_clad,core_diameter,lamda)
%a is the radius in um;
%lamda is the wavelenght in um;
a = core_diameter/2; %radius
k0 = 2*pi./lamda;
V = k0.*a.*sqrt(n_core.^2 - n_clad.^2);%V number
betta = zeros(1,length(V));
%%now we want to solve the following
% J0(x) / ( x*J1(x))=K0(sqrt(V^2-x^2) )/(sqrt(V^2-x^2)*K1(sqrt(V^2-x^2)))
% where x=Ka
[x,y]=arrayfun(@solveFn,V);
Kappa = x/a;
betta = sqrt(k0.^2.*n_core.^2-Kappa.^2);
neff=betta./k0;
function [x,y]=solveFn(V)
x=0.8*V;
b=sqrt(V*V-x*x);
[x,y]=fzero(@nFn,x);
% solution function nFn
function y=nFn(x)
y=besselj(0,x)/(x*besselj(1,x))-besselk(0,b)/(b*besselk(1,b));
end
end
This started with investigating the functional form of the function you're trying to solve as:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/177693/image.jpeg)
Here's the two pieces that make up the Diff vector for V(1) -- note particularly the humongously large values for small values of x; there's absolutely no sense in starting from 0 or even anywhere near there so a huge fraction of your computational overhead is just totally wasted looking in the wrong part of the world...as had suggested, even something as simple as a binary search to bracket the root would have helped immensely.
Blowing this up into an area of interest and adding the difference, one can see what the functional form is in the area of interest:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/177694/image.jpeg)
This shows the functions are quite smooth and so fzero given at least a reasonable starting point should have no trouble in finding the solution.
Note that the slopes are pretty sizable in the neighborhood of the solution; hence it was observed that the solution to near machine precision differs from your approximate solution by as much as 20% or so in the x value; this is observed in your Diff vector being very jaggedy if plotted. The above solution is quite smooth, in contrast.
I just used 0.8*V for the starting guess as a first cut; it ran through the entire vector with no problems that way; you could perhaps make it a little quicker if used the previous result as the subsequent guess as it changes very little between observations; some, in fact, were the same to machine precision between successive solutions altho the error changed noticeably.
Part of debugging session; just used the first two values for brevity--
>> nf2=neff2(n_core(1:2),n_clad(1:2),6,lamda(1:2));
K>> x(1)
ans =
1.4559
K>> y(1)
ans =
1.1102e-16
K>>
The solution using the min() was
K>> x=linspace(0,V(1),1E6);
K>> b=sqrt(V(1).^2-x.^2);
K>> F1=besselj(0,x)./(x.*besselj(1,x));
K>> F2=besselk(0,b)./(b.*besselk(1,b));
K>> [mn,ix]=min(abs(F1-F2))
mn =
8.7172e-07
ix =
813897
K>> x(ix)
ans =
1.4441
K>>
which is roughly 1% of magnitude but interestingly some
K>> x(ix)-x(ix+1)
ans =
-1.7742e-06
K>> 0.0149/ans
ans =
-8.3979e+03
K>>
steps removed from the next value computed. That's the result of the slopes being sizable so the difference even on such a small scale isn't that precise.
The errors in the neighborhood; you see it was indeed the smallest computed but fzero got quite a lot more accurate solution.
K>> [F1(ix-1:ix+1)-F2(ix-1:ix+1)]
ans =
1.0e-05 *
0.4511 0.0872 -0.2768
K>>
3 comentarios
dpb
el 6 de Ag. de 2017
Editada: dpb
el 7 de Ag. de 2017
y is saved solely in case it is ever wanted for debugging or further analysis; it isn't required or presently used; just kept because seemed quite possible it would be wanted in the future; doesn't cost anything else except a little memory which is cheap these days...
The "trick" using nested function to pass additional parameters to the minimization routines is documented here: <Parameterizing Functions>. Scope and details on nested functions is at <Nested Functions>
While I didn't check on it, you should be able to eliminate the for loop entirely with arrayfun
ADDENDUM
Well, just couldn't leave well enough alone... :)
See the updated Answer for using arrayfun() thereby eliminating the explicit loop and need for preallocation. I didn't compare run time so don't know if it makes any real difference or not.
ADDENDUM 2
I used the nested function alternative as I think it's a little easier to follow for the inexperienced user. As an "exercise for the student", your mission, should you choose to accept it, is to rewrite the solution using the anonymous function form which can reduce the source form to only the single line (I think w/o having actually done this case). :)
Más respuestas (1)
Walter Roberson
el 5 de Ag. de 2017
x=linspace(0,V(ii),1e6);
can be replaced by
x = V(ii) * linspace(0, 1, 1e6);
That in turn leads you to the vectorized
xG = bsxfun( V(:), linspace(0, 1, 1e6) );
and then
temp = sqrt( bsxfun( @minus, V(:).^2, x.^2) );
DiffG = abs(besselj(0, xG) ./ (xG .* besselj(1,xG)) - ...
besselk(0, temp) ./ temp .* besselk(1, temp)));
[~, ind1] = min(DiffG, [], 2); %per_row minima
kappa_ind = sub2ind(size(xG), 1:size(xG,1), ind1); %each ind1 value refers to a different row
kappa = xG(kappa_ind) ./ a; %his is a vector, one for each V
betta = sqrt(k0.^2 .* n_core.^2 - Kappa.^2);
should be the vector solution.
4 comentarios
Walter Roberson
el 5 de Ag. de 2017
Ah. In that case you will not be able to vectorize the calculation.
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!