Hi,

If you have a single optiona you can use Optional Ordererd Parameter (see ?parameter_classes)

something like that will do the job:

test:=proc(a,b,c:=NULL)
local d;
d:=a*b;
if c<>NULL then
print(c);
d:=d*c;
end if;
print(d);
end proc:

test(3,4,6);
                                      6
                                     72
test(3,5);
                                     15

 Otherwise, if you need something more complicated, see ?overload

S.

Hi all,

 

thanks for your replies. The error with Eigenvalues disappeared once I restarted maple (It was on from several days). Btw, it was M13 on a 64bit Linux pc.

I'm trying to use logarithms to avoid having large and small numbers, but the problem with the inverse seems to be still there when the dimension of the matrix is N>=7.

Are you aware of any reason that could explain a degradation of the inversion for N>=7?

 

Tkhs

Salvo

It seems that the threads in the "questions" page are ordered following the creation time, and not the last-reply time. This would explain why my answer at the older thread is not there. It seems an illogical way to order the threads, btw. Is there a way to change this behaviour?

Tks

Salvo

Alshom, thank you, but I don't feel comfortable with modules. I don't like that part of the maple advanced programming guide. Alec, than you for your suggestion. I'm going to think about it. Salvo

I forgot the ^2's in the sqrt, sorry

Hi Robert the arrays V's ARE global variables, the procedures' only meaning is to indexing these arrays, making some scalar products and put the results in the array entries. Actually, before the assignment I use a local variable to make the actual calculation, something like v5:=proc() description "calculates the tensor v_abcde"; local temp: global h41,h32,param_dim, V5, variables, parameters; V5:= Array(commutingindexes,1..param_dim,1..param_dim,1..param_dim,1..param_dim,1..param_dim); for ii from 1 to param_dim do for jj from 1 to ii do for kk from 1 to jj do for mm from 1 to kk do for nn from 1 to mm do temp:= h41[ii,jj,kk,mm,nn] + h41[ii,jj,kk,nn,mm] + h41[ii,jj,mm,nn,kk] + h41[ii,kk,mm,nn,jj] + h41[jj,kk,mm,nn,ii] + h32[ii,jj,kk,nn,mm] + h32[ii,jj,mm,nn,kk] + h32[ii,kk,mm,nn,jj] + h32[jj,kk,mm,nn,ii] + h32[ii,jj,nn,mm,kk] + h32[ii,kk,nn,mm,jj] + h32[jj,nn,kk,mm,ii] + h32[ii,mm,nn,kk,jj] + h32[mm,nn,jj,kk,ii] + h32[kk,mm,nn,jj,ii]; provv:= - simplify(temp); V5[ii,jj,kk,mm,nn]:=provv; end do end do end do end do; end do; end proc: is the use of temp unnecessary ? I use the FORs because the array is totally symmetric (communtingindex is the indexing function someone suggested to me here in the forum) and so I don't really need to assign all the values (like I'd be obliged to do using the assignment at the moment of the creation of the array). Is there any other way to do this, not using the FORs? Here the indexing function: `index/commutingindexes`:=proc(idx,M,val) description "indexing function for completely symmetric tensors"; local idx1: idx1:=op(sort(idx)); if nargs = 2 then # retrieving M[idx1]; else # storing M[idx1]:=op(val); end if: end proc: Thanks Salvo

Hi Scott

thanks for you reply. I feel like the program is too long to be published here, and additionally as it is a research work I'd prefer not to make it public.

I can sketch the main parts of the program, if it could be useful

 

1) a lot of self built procedures: mainly a weighted scalar product and the costructors for several arrays (let me call them V[ .. ]  )  from 5x5 dims up to 5x5x5x5x5 dims.

2) beginning of a for cycle

    3) I define a function "H" of 5 variables, that depends also on the value of the for index.
         I call the array V constructors (whose element (i,j ... k) is tipically a combination of derivatives of the function H with respect to the i-th, j-th ... k-th variable)
         The array are calculated (some three hours needed). Their values are saved or appended in external files

4) end of the for

At this point the for should begin another time, and (if I'm thinking in the right way) all the Arrays V's recalculated for the new value of the function H, replacing the old values, and so on.

But after the first run I get the stack error.

 

Salvo

Hi Axel

what you have written should work in principle; the only remark is that this line

map( integrator, % ); # apply the integration to each

is not as simple as it looks, because the integration procedure doesn't work on the expression itself,
but rather on the power of f, the coefficient, and the eventual other function. 
This means that I can't simply pass the expression to the integrator with the map procedure. 
I need before to extract the power, etc, before to pass them to the integrator. 
This is why I used three procedures to make this job. 

I'm going to append the whole code, if you want to look at it

Download 4492_mycode.mpl
View file details

Salvo

Hi Axel

this is, more or less what I want:

I've a signal which depends on four or five parameters (symbols), it's on the form p:= A*f^(-7/6) exp( I * (combination of parameters) ).

I've an array, whose entry "i" (let' me call it h[i]) is the derivative of the signal with respect to the i-th parameter,   calculated for numerical values of the parameters (ie: the frequency f is the only symbol left now)

I want scalar product (in the sense I've defined in my earlier post) between h[i] and h[j] for all possible values of i and j (ie 1.. lenght_parameters), the bounds of the integration are defined in a two terms list Extrema:=[20..some hundreds] (the results will be putten in a two-dim array)

The procedure scal prod makes a combination of his two parameters A, B ---> Re(AB* +A*B) (* means conjugation) and should integrate this combination over the normalization function S_h(f) for f within Extrema.
The general terms of this integral will be

C*D* f^a g(f)/S_h(f)

with C real, D symbol.

Now you have suggested a brilliant way to deal with this integral. But normally the real integral will be a sum of terms like this one, for different C,D,a,g(f), depending on what derivative of the signal I'm considering.

This is why I tried to extend you method, making procedures that break the sum inside the integral (splitter) , understand  for each term the values of C,D,a, g(f) (atomizer), and pass them to your integral proc (that works with one monomium at the time)
 

 I hope now is clearer what I need, is difficult to explain it without a blackboard, you're right !

Thanks

Salvo

Hi Axel

I'm trying  to extend the possibilities of the method with a little modification in the procedure:

myInt_DE_g:=proc(a,g,b,c, DigitsForWorking)
local x, f0;
Digits:=DigitsForWorking;
f0:=215;
if g<>1 then  ### g=1 if no supplementary function is present
f0^(1+a)*intDE('x -> eval(g,f=f0*x)*K(x,a)',f0/c,f0/b, Digits);
else
f0^(1+a)*intDE('x -> K(x,a)',f0/c,f0/b, Digits);
end if;
end proc;

where g can be one of these functions:

ln(f), (ln(f))^2, 1/(1+x^2)

but it seems that there is some problem when the logaritm is involved (the result of myInt are smaller than those of direct evalf(Int - ))
Is that in relation with this particular numerical routine?

Thanks

Salvo

Hi Doug

thanks for your suggestions, I'm going to read your post very carefully.
 

S.

ps: if you want to think about it a little more, b =20.0, c is variable with the mass of my system (yes, it's physical problem), it starts from 2500-3000 and decrease till it reaches the value of b (but I'd be happy to neglect this patological point and stop it before)

Hi everyone after some attempts, I think I've understood that my program could me more efficient if there was a way to symbolically solve this integral for real values of a, and positive b and c: S_h:=proc(t) local f_0,T; f_0 :=215: T:=t/f_0: simplify((T^(-414/100) -5*T^(-2)+ 111*(1-T^2+T^4/2)/(1+T^2/2))) end proc: assume(a,'real'); assume(b,'positive'); assume(c,'positive'): result:= int(f^a/S_h(f), f = b .. c); Do you think is possible? Should I use some hints ? Thanks S.

Hi Doug

 

I use IntegrationTools:-Expand because I was told to follow that way here in the forum! This was becaus of the symbols in the integrand.

I can copy a sheet of the code, till the scalar product functions. Thank for any help

restart;

time_start:= time():

Digits:=50:

S_h:=proc(t)
local f_0,T;
global S_0:
#S_0:=1.0*10^(-49):
f_0 :=215.0:
T:=t/f_0:
simplify(S_0*(T^(-4.14) -5*T^(-2)+ 111*(1-T^2+T^4/2)/(1+T^2/2)))
end proc:

assume(eta,'positive');
assume(phi,'positive');
assume(t,'positive');
assume(Enn,'positive');
assume(M,'positive');
assume(f,'positive');
assume(Emm,'positive');
assume(S_0,'positive');
assume(rho,'positive');
with(LinearAlgebra):

h1_digits:=50;
h2_digits:=50;
h3_digits:=50;


parameters:= [t,phi,Emm,eta]:
param_dim:= nops(parameters):
f_low:=20.0:

M_solar:=4.927*10^(-6):
eta_eff:= .25:
PNorder:=2:
M:=0.0;

`index/commutingindexes`:=proc(idx,M,val)
 local idx1:
 idx1:=op(sort(idx));
 if nargs = 2 then
   # retrieving
   M[idx1];
 else
   # storing
   M[idx1]:=op(val);
 end if:
end proc:

for counter from 0 by 1 to 19 do ##################### FOR CYCLE

M_tot:=(2.8 + 219.0*counter/20)*M_solar:
appendto("./M_tot_seq.txt"):
lprint(M_tot);
appendto("./times.txt");
printf("Counter value is %d \n", counter);
appendto("./garbage.txt"):

f_lso:= (6^(3/2)*Pi*M_tot)^(-1):
v_lso:=(Pi*M_tot*f_lso)^(1/3):
f_high:=f_lso:

variables:=[f]:
variab_dim:=nops(variables):

p:= proc()
option hfloat;
local Psi,i,psi,alpha_par,v,theta,euler_gamma,lambda;
global Enn,PNorder,param_dim,f,eta,M,M_tot,f_lso,v_lso,Emm,eta_eff,parameters,variables,S_0,rho,t,phi;
Digits:=50;
M:=Emm*eta^(-3/5);
v:=(Pi*M*f)^(1/3);
theta:=-1.28;
lambda:=-0.6451;
euler_gamma:=0.577215664901532860606512090082:
alpha_par:=[
1
,0
,20/9*(743/336 + 11/4*eta)
,-16*Pi
,10*(3058673/1016064+5429/1008*eta+617/144*eta^2)
, Pi*(38645/756+38645/252*ln(v/v_lso) -65/9*eta*(1+3*ln(v/v_lso)))
, (11583231236531/4694215680-640/3*Pi^2-6848/21*euler_gamma)+eta*(-15335597827/3048192+2255/12*Pi^2-1760/3*theta+12320/9*lambda)+76055/1728*eta^2-127825/1296*eta^3-6848/21*ln(4*v)
,Pi*(77096675/254016+378515/1512*eta-74045/756*eta^2)];
Psi:= (2*Pi*f*t-phi-Pi/4+ 3/(128*eta*v^5)*add(alpha_par[i+1]*v^i, i=0..PNorder));
Enn*f^(-7/6)*exp(I*Psi):
 end proc:

H1:=proc()
global Enn,PNorder,param_dim,f,eta,M,M_tot,f_lso,v_lso,Emm,eta_eff,parameters,variables,S_0,rho,t,phi,h1, counter_h1,h1_digits;
Digits:=h1_digits;
h1:=Array(1..param_dim, i->
'eval'(
(diff(p(), parameters[i]))
,[Emm=(M_tot*(eta_eff)^(3/5)),eta=eta_eff] )
):
end proc;

time_h1:=time(H1());
appendto("./time-h.txt"):
printf("h1 has been calculated in %a seconds with %d digits \n", time_h1, h1_digits);
appendto("./garbage.txt"):

prova:=simplify(h1):
save prova,"./h1":


scalar_pr:= proc(A,B)
option hfloat;
global variables, parameters,f_low, f_high,S_h,M_tot,eta_eff,S_0,f_lso,f,Enn,rho;
local num;
description "Calculate the scalar product between two tensors A and B in the frequency space":
Digits:=15;
num:=
Re(
simplify(
A*conjugate(B)+conjugate(A)*B
)
)
;
evalf(
IntegrationTools:-Expand(
2.0*Int(num/S_h(f), f=f_low..f_high)
)
);
appendto("./time_snr.txt");
lprint(time_snr_end);
appendto("./garbage.txt"):

Low_fisher:=proc()
option hfloat;
description " Calculate the fisher information matrix called low_fisher(a,b)";
global Enn,low_fisher,S_0, f_low,f_high,param_dim,M_tot,eta_eff,parameters,variables,rho,t,phi,M,eta,Emm:
local gg,a_1,b_1:
gg:=(a_1,b_1)->
scalar_pr(h1[a_1],h1[b_1])
;
low_fisher:= Matrix(param_dim,param_dim, shape=symmetric,gg);
end proc:

Low_fisher():

#save low_fisher, "./low_fisher.txt":

Up_fisher:=proc()
option hfloat;
description " Calculate the upper index fisher information matrix called up_fisher(a,b)";
global Enn, up_fisher,M_tot,eta_eff,S_0,param_dim,variables,parameters,rho,t,phi,Emm,eta,M:
local a_3,b_3,gg;
gg:=(a_3,b_3)-> (MatrixInverse(low_fisher))[a_3,b_3];
up_fisher:= Matrix(param_dim,param_dim,shape=symmetric,gg);
end proc:

Up_fisher():
save up_fisher,"./up_fisher":
end proc:

time_snr_in:=time():
snr:=scalar_pr(p(),p()):
time_snr_end:=time()-time_snr_in:
save snr, "./snr_squared.txt":
appendto("./snr-seq.txt"):
snr;
 
end do;

Hi everybody I'm a little bit confuse about how maple works with local values of Digits. Until yesterday my program was schematically like follow: I had a function h(a,b,c,d) with both analytical and numerical coefficients of this kind h = a f^(-7/6) exp( I * (2*pi*b + c + etc) and his derivatives (first to third) evaluated for numerical values of his arguments. a "scalar product" of two derivatives of h that I wrote as a procedure: scalar_pr:= proc(A,B) option hfloat; global variables, parameters,f_low, f_high,S_h,M_tot,eta_eff,S_0,f_lso,f,Enn,rho; local num; num:= Re(simplify(A*conjugate(B)+conjugate(A)*B)); evalf(IntegrationTools:-Expand( 2.0*Int(num/S_h(f), f=f_low..f_high) )); end proc; where now all the quantities (integral extrema, S_h, etc) have numerical values, and the integrand has an overall symbolic coefficient with is pulled out by integrationTools. Then, later in the program, I build tensor with N indexes (N=1..4) which are (roughly speaking) sum of scalar products of derivatives of h with all the possible combination of that particular set of index, for exemple: v21:=proc() option hfloat; global h1,h2,h3,param_dim, V21,variables, parameters,f_low, f_high,S_h,M_tot,eta_eff,S_0,f_lso,f,Enn,rho,Emm,eta,M,t,phi; V21:=Array(1..param_dim,1..param_dim,1..param_dim, (ii,jj,kk)-> simplify( scalar_pr(h2[ii,jj],h1[kk]) )): end proc: Then these tensors were summed in a particular combination, to give a final numerical answer. So it seemed evident to me that the most important part of my code was the scalar product, as all the other function would use it. It turned out that with Digits<=15 I had a good performance, but bad numerical answers (too few digits), and that with (global) Digits>15 bad performance, as I complained before in this thread. I tried to increase digits only in some functions, in particular, in particular in those which calculate the derivatives of h, and those which calculate the tensors. But the result suffered of the same numerical errors. What changed the situation was to put Digits:=30 elsewhere but in the scalar_proc procedure where I left Digits:=15 locally. With this choice the programs didn't loose time-efficiency and gained numerical efficiency. So my first question is: why is my scalar product procedure so sensible to the actual value of digits and not the rest of the code? The second is: are the results that I obtain now (with 30 digits) really more precise, or the fact that Digits in scal_prod is still 15 means that the "true" faithful digits in the final results are still 15? Thank for any hints S.

Thank to both of you. Yes, the efficiency is an issue (a big one!!), so I'll try Pagan's suggestion. S.