Thanks!

yes, this can do the job! Thanks S.

yes, this can do the job! Thanks S.

Don't worry about that. I know that the program is resources-consuming, I don't want to destroy your pc :-)

I think that now, expecially because of your help, it runs pretty well.

I think that I  could still spare some minutes (over a several hours running), but nothing sensational, maybe it isn't worth. The big part of the work has been done already.

Thanks!

Salvo

Don't worry about that. I know that the program is resources-consuming, I don't want to destroy your pc :-)

I think that now, expecially because of your help, it runs pretty well.

I think that I  could still spare some minutes (over a several hours running), but nothing sensational, maybe it isn't worth. The big part of the work has been done already.

Thanks!

Salvo

Hi Robert Is it possible to run the program in a "diagnostic way"? I know that in the mw or classic mode there are "stopat" and other tools. Nothing similar for the non-interactive programs? Thanks S.

Hi Robert Is it possible to run the program in a "diagnostic way"? I know that in the mw or classic mode there are "stopat" and other tools. Nothing similar for the non-interactive programs? Thanks S.

I also remarked that sometimes in the directory that contains the .mpl a file called core.ANumber (like core.26907) appears (I don't know if while the program is running or when it ends). It seems to be a binary file. Any idea about that? Thanks S.

Hi Axel I modified the code in order to make it able to deal with integrand of the form: f^a* g(f)/S_h(f) where g(f), as I told you, is typically ln(f), ln(f)^2, or a lorentzian 1/((f- A)^2-B) it turned out that the only modification in my_int_DE that one needs is to multiply the integrand you proposed by the new function evaluated at f0/x: myInt_DE_modif:=proc(a,g,b,c, DigitsForWorking) description "improved numerical integration of the function f^a g(f)/S_h(f) f=b..c"; local x, f0; Digits:=DigitsForWorking; f0:=215; if g<>1 then #### g present f0^(1+a)*intDE('x -> g(x)*K(x,a)',f0/c,f0/b, Digits); else ##### g(f)=1 ie not present, it computes the old integral. f0^(1+a)*intDE('x -> K(x,a)',f0/c,f0/b, Digits); end if; end proc; The reason I'm writing is to have your advice (other advices are welcomed, obviously!), abut the following code. It works well, but I'd like to know if it is well written, or if it could be more efficient, changing it somehow. The first procedure is the scalar product. It take two inputs and the extrema of integration (as a list). If A and B are scalar it calculate the scalar product by creating a combination of A and B, and passing it to the function splitter (see below) . If A and B are list, of the same dimension, the procedure calls himself with A[1],B[1] as scalar parameters, does the same for each A[n],B[n] and sums the results. scalar_pr:= proc(A,B,extrema::list) description "Calculate the scalar product between A and B in the frequency space. It check whether A and B are lists or not and acts smartly": global variables, parameters,M_tot,eta_eff,f_lso,f,Enn,rho,scal_pr_digits; local num, noA,noB,counter,temp,result; Digits:=scal_pr_digits; if (type(A,'scalar') and type(B,'scalar')) then num:= Re( simplify( A*conjugate(B)+conjugate(A)*B ) ); return 2*splitter(num,extrema); elif (type(A, 'list') and type(B,'list')) then #### if A and B are both lists noA:=nops(A); noB:=nops(B); if noA <> noB then error "The two lists should have the same dimension!"; end if; counter:=1; temp:=0; result:=0; while counter<=noA do temp:=scalar_pr(A[counter],B[counter],[extrema[counter],extrema[counter+1]]); result:=result+temp; counter := counter+1; end do; return simplify(result); else error "I don't recognize your inputs. Are they both lists or scalar?"; end if; end proc: The next function is called splitter, it simply checks if the first of his parameters (the second being only the integration extrema) is a sum or not. If it is a sum, it calls himself with the first term of this sum like first parameter, does the same with all the addends, and sums the result. Shortly: it makes a sum of integrals from the integral of a sum! It makes the integration calling the function my_Int_DE. But in order to make it, it needs a decomposition of his argument, that will be on the general form N*f^a* g(f). For this purpose, it calls the procedure atomizer (see below). if g(f) is present, it convert it in a procedure caos:= x-> g(f0/x), if not g(f)=1. splitter:=proc(Num,extrema::list) description "splits the integral of the sum in a sum of integrals"; local temp,Dig,partial_int,result,f0; global myInt_DE,non_f,f_pow,has_numbers,f; temp:=expand(Num); n_temp:=nops(temp): result:=0; partial_int:=0; f0:=215; if hastype(temp,`+`) then for i from 1 to n_temp do loc:=expand(op(i,temp)); partial_int:=splitter(loc,extrema); result:=result+partial_int: end do: return result; else loc:=expand(temp); atomizer(loc); #### the call to atomizer changes the global vars non_f,f_pow, and has_numbers if is(non_f = 1) then caos:=1; else tmp:=subs(f=f0/x,non_f); caos:=unapply(tmp,x); end if; tot_result:=has_numbers*myInt_DE_ln(f_pow,caos,extrema[1],extrema[2],Digits); ### integration call return simplify(tot_result): end if; end proc: The last procedure is called atomizer. It's called by splitter and does the following: it takes his parameter and create three global variables: 1) f_pow ie the power "a" in f^a if f^a is present in the parameter. 2) non_f ie an eventual second function that is not simply on the form f^a (ie g(f)) 3) has_numbers ie the overall factor that can be numerical or symbolic. if one of more of this quantities are non present on the input, it uses default values: 0 for f_pow (f^0=1 ie not present), 1 for non_f, and 1 for has_numbers. atomizer:=proc(ex) description "searches the constituents of his parameter, the power of f^a the constant coeff, and other function (ln,lorentzian). Used in splitter()"; global f_pow,has_numbers,non_f,f; local pattern, result,temp,loc_has_numbers,loc_f_pow,i,temp_has_numbers; pattern:=f^a; result:=1; temp:=1; loc_has_numbers:=1; loc_f_pow:=0; for i in op(..,ex) do if not has(i,f) then temp_numbers:=i; loc_has_numbers:= loc_has_numbers*temp_numbers: end if; if has(i,f) then if match(i = pattern, f, 'S') then loc_f_pow := eval(a, S); else temp:=i; result:=result*temp; end if; end if; end do; non_f:=result; has_numbers:=loc_has_numbers: f_pow:=loc_f_pow: return; end proc: I guess that is not terrific that atomizer modifies global variables instead of passing them to his caller. But to obtain this behavior, it should make atomizer a module and export the three variables. It that worthwhile? Any other remarks? Thank for you comments and help, I'm very grateful Salvo

Hi Axel I modified the code in order to make it able to deal with integrand of the form: f^a* g(f)/S_h(f) where g(f), as I told you, is typically ln(f), ln(f)^2, or a lorentzian 1/((f- A)^2-B) it turned out that the only modification in my_int_DE that one needs is to multiply the integrand you proposed by the new function evaluated at f0/x: myInt_DE_modif:=proc(a,g,b,c, DigitsForWorking) description "improved numerical integration of the function f^a g(f)/S_h(f) f=b..c"; local x, f0; Digits:=DigitsForWorking; f0:=215; if g<>1 then #### g present f0^(1+a)*intDE('x -> g(x)*K(x,a)',f0/c,f0/b, Digits); else ##### g(f)=1 ie not present, it computes the old integral. f0^(1+a)*intDE('x -> K(x,a)',f0/c,f0/b, Digits); end if; end proc; The reason I'm writing is to have your advice (other advices are welcomed, obviously!), abut the following code. It works well, but I'd like to know if it is well written, or if it could be more efficient, changing it somehow. The first procedure is the scalar product. It take two inputs and the extrema of integration (as a list). If A and B are scalar it calculate the scalar product by creating a combination of A and B, and passing it to the function splitter (see below) . If A and B are list, of the same dimension, the procedure calls himself with A[1],B[1] as scalar parameters, does the same for each A[n],B[n] and sums the results. scalar_pr:= proc(A,B,extrema::list) description "Calculate the scalar product between A and B in the frequency space. It check whether A and B are lists or not and acts smartly": global variables, parameters,M_tot,eta_eff,f_lso,f,Enn,rho,scal_pr_digits; local num, noA,noB,counter,temp,result; Digits:=scal_pr_digits; if (type(A,'scalar') and type(B,'scalar')) then num:= Re( simplify( A*conjugate(B)+conjugate(A)*B ) ); return 2*splitter(num,extrema); elif (type(A, 'list') and type(B,'list')) then #### if A and B are both lists noA:=nops(A); noB:=nops(B); if noA <> noB then error "The two lists should have the same dimension!"; end if; counter:=1; temp:=0; result:=0; while counter<=noA do temp:=scalar_pr(A[counter],B[counter],[extrema[counter],extrema[counter+1]]); result:=result+temp; counter := counter+1; end do; return simplify(result); else error "I don't recognize your inputs. Are they both lists or scalar?"; end if; end proc: The next function is called splitter, it simply checks if the first of his parameters (the second being only the integration extrema) is a sum or not. If it is a sum, it calls himself with the first term of this sum like first parameter, does the same with all the addends, and sums the result. Shortly: it makes a sum of integrals from the integral of a sum! It makes the integration calling the function my_Int_DE. But in order to make it, it needs a decomposition of his argument, that will be on the general form N*f^a* g(f). For this purpose, it calls the procedure atomizer (see below). if g(f) is present, it convert it in a procedure caos:= x-> g(f0/x), if not g(f)=1. splitter:=proc(Num,extrema::list) description "splits the integral of the sum in a sum of integrals"; local temp,Dig,partial_int,result,f0; global myInt_DE,non_f,f_pow,has_numbers,f; temp:=expand(Num); n_temp:=nops(temp): result:=0; partial_int:=0; f0:=215; if hastype(temp,`+`) then for i from 1 to n_temp do loc:=expand(op(i,temp)); partial_int:=splitter(loc,extrema); result:=result+partial_int: end do: return result; else loc:=expand(temp); atomizer(loc); #### the call to atomizer changes the global vars non_f,f_pow, and has_numbers if is(non_f = 1) then caos:=1; else tmp:=subs(f=f0/x,non_f); caos:=unapply(tmp,x); end if; tot_result:=has_numbers*myInt_DE_ln(f_pow,caos,extrema[1],extrema[2],Digits); ### integration call return simplify(tot_result): end if; end proc: The last procedure is called atomizer. It's called by splitter and does the following: it takes his parameter and create three global variables: 1) f_pow ie the power "a" in f^a if f^a is present in the parameter. 2) non_f ie an eventual second function that is not simply on the form f^a (ie g(f)) 3) has_numbers ie the overall factor that can be numerical or symbolic. if one of more of this quantities are non present on the input, it uses default values: 0 for f_pow (f^0=1 ie not present), 1 for non_f, and 1 for has_numbers. atomizer:=proc(ex) description "searches the constituents of his parameter, the power of f^a the constant coeff, and other function (ln,lorentzian). Used in splitter()"; global f_pow,has_numbers,non_f,f; local pattern, result,temp,loc_has_numbers,loc_f_pow,i,temp_has_numbers; pattern:=f^a; result:=1; temp:=1; loc_has_numbers:=1; loc_f_pow:=0; for i in op(..,ex) do if not has(i,f) then temp_numbers:=i; loc_has_numbers:= loc_has_numbers*temp_numbers: end if; if has(i,f) then if match(i = pattern, f, 'S') then loc_f_pow := eval(a, S); else temp:=i; result:=result*temp; end if; end if; end do; non_f:=result; has_numbers:=loc_has_numbers: f_pow:=loc_f_pow: return; end proc: I guess that is not terrific that atomizer modifies global variables instead of passing them to his caller. But to obtain this behavior, it should make atomizer a module and export the three variables. It that worthwhile? Any other remarks? Thank for you comments and help, I'm very grateful Salvo

Ok, I've founded something satisfying, that allow me to create three variable, the (eventual) numerical overall coefficent, the power of f and the other eventual function of f

prova:=proc(ex)
global f_pow,has_numbers,non_f;
local pattern, result,temp,loc_has_numbers,loc_f_pow,i;
pattern:=f^a;
result:=1;
temp:=1;
loc_has_numbers:=1;
loc_f_pow:=0;
for i in op(..,ex) do

if not has(i,f) then
  loc_has_numbers:=i;
end if;

if hastype(i,symbol) then
   if match(i = pattern, f, 'S') then   
     loc_f_pow := eval(a, S);
   else
     temp:=i;
     result:=result*temp;
  end if;
end if;
end do;
non_f:=result;
has_numbers:=loc_has_numbers:
f_pow:=loc_f_pow:
return;
end proc:


en:=-4*f^(-3)*ln(f)^2/(f+1);

                                          2
                                   4 ln(f) 
                                - ----------
                                   3       
                                  f  (f + 1)

prova(en);

has_numbers;
f_pow;
non_f;
                                     -4
                                     -3
                                        2
                                   ln(f)
                                   ------
                                   f + 1
 

This is enough for me purposes, I don't need the function non_f to be further decomposed.

Thank to all

 

Salvo

Ok, I've founded something satisfying, that allow me to create three variable, the (eventual) numerical overall coefficent, the power of f and the other eventual function of f

prova:=proc(ex)
global f_pow,has_numbers,non_f;
local pattern, result,temp,loc_has_numbers,loc_f_pow,i;
pattern:=f^a;
result:=1;
temp:=1;
loc_has_numbers:=1;
loc_f_pow:=0;
for i in op(..,ex) do

if not has(i,f) then
  loc_has_numbers:=i;
end if;

if hastype(i,symbol) then
   if match(i = pattern, f, 'S') then   
     loc_f_pow := eval(a, S);
   else
     temp:=i;
     result:=result*temp;
  end if;
end if;
end do;
non_f:=result;
has_numbers:=loc_has_numbers:
f_pow:=loc_f_pow:
return;
end proc:


en:=-4*f^(-3)*ln(f)^2/(f+1);

                                          2
                                   4 ln(f) 
                                - ----------
                                   3       
                                  f  (f + 1)

prova(en);

has_numbers;
f_pow;
non_f;
                                     -4
                                     -3
                                        2
                                   ln(f)
                                   ------
                                   f + 1
 

This is enough for me purposes, I don't need the function non_f to be further decomposed.

Thank to all

 

Salvo

Hi Robert

 

thanks for your suggestion, it seems very useful. But I'm having some difficulties trying to extend it in more general scenarios.

for ex if I expect to have expression of this kind

c*f^a*ln(f)^b /( (f-d)^2 + e)^g

 

where positives and negatives (and zero) values are allowed for a,b,c,d,e,g

if I use your filter in this way:

 

pattern := c*f^a*ln(f)^b /( (f-d)^2 + e)^g:

if match(en = pattern, f, 'S') then
    thec := eval(c, S); 
    thepower := eval(a, S);
the_ln_power:=eval(b,S);
thed:=eval(d,S);
thee:=eval(e,S);
theg:=eval(g,S);
  else
    error "%1 doesn't match the pattern", A
  end if;

I should be able to match this expression

A:=32*f^(2)*(ln(f))/(f-.1)^2

with c=32, a=2, b=1, d=0.1, e=0, g=1

but the if returns the error in this case. It seems to have problems with the value 0.1 for d.

Have you some idea about this behaviour?

 

Salvo

Hi Robert

 

thanks for your suggestion, it seems very useful. But I'm having some difficulties trying to extend it in more general scenarios.

for ex if I expect to have expression of this kind

c*f^a*ln(f)^b /( (f-d)^2 + e)^g

 

where positives and negatives (and zero) values are allowed for a,b,c,d,e,g

if I use your filter in this way:

 

pattern := c*f^a*ln(f)^b /( (f-d)^2 + e)^g:

if match(en = pattern, f, 'S') then
    thec := eval(c, S); 
    thepower := eval(a, S);
the_ln_power:=eval(b,S);
thed:=eval(d,S);
thee:=eval(e,S);
theg:=eval(g,S);
  else
    error "%1 doesn't match the pattern", A
  end if;

I should be able to match this expression

A:=32*f^(2)*(ln(f))/(f-.1)^2

with c=32, a=2, b=1, d=0.1, e=0, g=1

but the if returns the error in this case. It seems to have problems with the value 0.1 for d.

Have you some idea about this behaviour?

 

Salvo

Hi Axel

I've found the problem, the change of variable implies that the argument of log should be

ln(f0/x)

and not

ln(x*f0)
 

as I've wrongly deduced by the fact that the change of variable was x = t/f0 --> t= x*f0

So now it work with the log, and I'm trying with the lorentzian 1/(1+x^2) that should be good behaved.

I'll tell you

 

thanks

Salvo