I have to numerically solve this equation many (tens of thousands) of times for theta_n (as I vary different parameters, especially n):
tan(theta_n) - C_a*Z_0*(-v^2*(Pi*n-theta_n)^2/x_l^2+omega_a^2)*x_l/(v*(Pi*n-theta_n)) = 0
theta_n should be in [-pi/2, pi/2). It seems like solving this is the slowest single component of the chain of calculations (that follow this).
Currently I do it with this function:
get_theta_n_array:=proc(max_n::integer, omega_a::float, v::float, x_l::float, C_a::float, Z_0::float)
select(x->Re(x)<evalf(Pi/2) and Re(x)>=-evalf(Pi/2), [seq(
#NOTE: careful with fsolve - in some cases returns unevaluated equation
fsolve(tan(theta_n) - C_a*Z_0*(-v^2*(Pi*n-theta_n)^2/x_l^2+omega_a^2)*x_l/(v*(Pi*n-theta_n)) = 0, theta_n) , n=1..max_n)]
if ArrayNumElems(theta_n_array) <> max_n then
printf("Bad Array Dimensions! Got too many or not enough solutions.");
theta_n_array:="CHECK: get_theta_n_array()": #dirrrrty hack that will ring an alarm bell if array is not the right size
And call it like so (for say n=1000)
result:=get_theta_n_array(1000, 100e9, 1e8, 0.3, 20e-14, 50.0);
This will take say 3.5s on my PC.
Does anyone have any ideas how to speed this up? I would hope this to take at least an order of magnitude less time. I played with DirectSearch lib but that was not faster.
Also, I should note that this is the only portion of my code that is not thread safe (because of the fsolve call), which leads to "extra" slowdowns because I have to use Grid:-Map, instead of Thread:-Map when parallelizing, and more importantly because I can't compile the rest of the code (Grid:-Map is not compatible with compiled functions).
Let me know if you have any ideas...