Hi

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)
    local theta_n_array:=Array(
        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)]
        , real)
        ,datatype=float[8]);
    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
    end;
    theta_n_array;
end;

And call it like so (for say n=1000)

st:=time();
result:=get_theta_n_array(1000, 100e9, 1e8, 0.3, 20e-14, 50.0);
time()-st;

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...

thanks!

I have a system of linear differential equations and am trying to solve them using Fourier transforms.

I can reduce the system to a final result (for the variable of interest) in Fourier space as this (note the frequency variable is 'w'):

restart:

with(inttrans):

vout_fourier_num := fourier(phi[3](t), t, w) = 6.63569999999998*10^(-15)*w^2*fourier(V(t), t, w)/(-5.69875218358308*10^(-40)*w^4+(9.19473390627057*10^(-29)*I)*w^3+2.15219369729956*10^(-18)*w^2-(4.14691648617110*10^(-8)*I)*w-700.8);

#the drive can be defined as:
drive:=5.70000000000000*10^(-6)*exp(-3.18877551020408*10^18*(t-2.0*10^(-9))^2)*cos(4.8*10^9*Pi*(t-2.0*10^(-9)))*10^(19/20);

#substitute the drive in - this is not necessary, but it should work!...
vout_fourier_num2:=subs(V(t)=drive, vout_fourier_num);

#now take the inverse... note this gives 0!...
invfourier(vout_fourier_num2, w, t);

the final results calculated is zero. It is wrong... it seems like an accuracy issue, but increasing the digits does not help. I should note that I can calculate the solution directly via dsolve, and get completely reasonable answer.

any ideas how to get the Fourier method to work?

thanks!

Say I have some function that cannot be changed and it returns this:

all_plots:=display([plot(sin(10*x+0.2),x=0..1, thickness=10, color=blue), plot(1-sin(10*x),x=0..1, thickness=10, color=red), plot(sin(10*x),x=0..1, thickness=10, color=green)]):


Now can plot it as:
plots:-display( all_plots );


...but what if I need the (say) red curve to be "on top" (fully visible), green curve in the middle, and blue curve at the bottom (as is now).

what is the most effective way to do this?

EDIT:

i guess this works:

display([convert(all_plots, list)[1], convert(all_plots, list)[3],  convert(all_plots, list)[2]]);

EDIT 2:

... lookst like above doesn't keep all the extra options that a plot can have... here is a version that seems to do it.

rearrangeCurves:=proc(v_items, v_reorder:=[])
  #Reorder should be a list of index pairs like
  #[[3,5], [-1, 1], ...]
  local p, temp, curves:=[], rest:=[]:
  #separate curves from rest (there must be a prettier way to do this)
  for p in convert(v_items, list) do
    if type(p,'function') and op(0,p) = ('CURVES') then 
      curves:=[curves[], p]:
    else
      rest:=[rest[], p]:
    end;
  end;
  #now reorder
  for p in v_reorder do
    temp:=curves[p[1]]:
    curves[p[1]]:=curves[p[2]]:
    curves[p[2]]:=temp:
  end;
  PLOT(curves[], rest[]):
end:

#change say second last with last curve:

rearrangeCurves(all_plots, [[-2, -1]]);

is there a built in command to get rid of terms of some order or greater in a given variable?

say i have:

myeq:=diff(mu[m](t), t, t) = (-2*d^2*phi[0]^2-2*phi[0]^2)*mu[m](t)/(2*C_J*phi[0]^2*L)+(J*L*sin(mu[m](t)/phi[0])*d^2+J*sin(mu[m](t)/phi[0])*L)*mu[p](t)^2/(2*C_J*phi[0]^2*L)+J*cos(mu[m](t)/phi[0])*d^3*mu[p](t)/(phi[0]*C_J)+(I_b*L*d*phi[0]^2-2*J*L*sin(mu[m](t)/phi[0])*phi[0]^2-d^2*Phi[xx]*phi[0]^2-phi[0]^2*Phi[xx])/(2*C_J*phi[0]^2*L);

and i want to get rid of all terms of O(d^2).

i know i can do this

lhs(myeq)=algsubs(d^2=0, rhs(myeq));

or

lhs(myeq)=convert(series(rhs(myeq), d=0, 2), polynom);

but I wonder if there is a built-in command. 

thanks

I have a complicated function that I call on a elements of a list. This function inside calls various other functions and in particular fsolve, which is not thread safe. In order to speed things up I tried using Grid:-Map, but have found a slight drop in performance.

So say a run with map takes 31 seconds, a run with Grid:-Map will take 32 seconds. By looking at the cpu activity, I can tell that only a single core is used at any one time (this is also evident from timing). By looking at the processes, I do see "mserver -gridnode 1/4", etc. present however, so new mserves are certainly spawned - they are just just not used concurrently.

Grid:-NumNodes() gives 4 for my CPU.

Any ideas how to force my code to run on all four nodes??

I am on ubuntu 12.04 x86_64.

thanks.

 ** EDIT **

In case someone finds this useful. Here is a quick-hack that fixes this problem for me.

Look at the current implementation of Grid:-Map with

showstat(Grid:-Map):

Define a new function called say hackMap that is identical to Grid:-Map, except, change the second line from this

     var := `union`(`union`(map(proc (x) op(0,x) end proc,indets(f,function)),indets(f,name)),{'`grid/mapcmd`'});

to this:

     var := {anames('user'), '`grid/mapcmd`'};

now use hackMap as you would Grid:-Map.

It basically forces maple to copy all the 'user' variables over to the new processes. On a machine with 4 nodes i get a ~4 seedup with my particular code.

NOTE: I am not certain that there aren't some side effects as I dont' know all the code behind the scenes, but I seem to get correct results with my particular code.