Hi Doug thanks for your help In the meantime I've found a way, using for cycles (the reason I want to use for is that my array's entries are summations of terms, and I want to check for each of this addend whether it is positive or not): par_dim:=4: a:=proc() global aaa; for i from 1 to par_dim do for j from 1 to i do for k from 1 to j do for m from 1 to k do aaa[i,j,k,m]:= 2^i+3^j+4^k+5^m; for ind in combinat:-permute([i,j,k,m]) do aaa[op(4,ind),op(3,ind),op(2,ind),op(1,ind)]:=aaa[i,j,k,m]; end do end do end do end do end do; end proc; and it works well! :-) S.

To be clearer, I know that I can initialize the independent entries of the array with something like for i from 1 to 4 do for j from 1 to i do for k from 1 to j do for m from 1 to k do V[ i, j, k, m]:= something: end do end do end do end do; but how to tell maple that it can initialize the other entries by simply commutating the index of these ones? S.

Hi I tried to use interface(displayprecision=2) ; -1 in both document and worksheet mode (maple 13) but I alway obtain a wrong behavior a:=evalf(1/3); a := .33333333333333333333 with ten digits How so? Thanks S.

I solved myself

I shoud write

if is(a>0) then etc...

S.

Hi

 

can I use as label the solar mass symbol (see http://en.wikipedia.org/wiki/Solar_mass) ?

Thanks.

S.

Additionally, how can I tell to Maple to write the numbers on the axis much bigger?

 

thanks

 

it works, thanks!

I use maple 13.01 on Linux Fedora 11, using worksheet with maple notation input, and 2D output

 

S.

Ok, I discovered labelfont=[SYMBOL] and that solve for eta. I could try the italic for M, but how to assign different labelfonts for different axis?

 

Thaks

S.

hi

 

I tried with loglogplot but I still get zeroes in the vertical axis :

 

 

 can't I have exponential notation in the vertical axis??

Thanks

Hi

I've thought about my problem, and I can very well understand and accept that results will always have a certain degree of error  (even an important one) due to round off approximations. What I can't accept is that the approximations done will always be different, and will lead to different final results. Is Maple similar to quantum mechanics??

If Maple chooses for ex. 2/3 = 0.6666666667  the first time I execute the program, it shouldn't assume 2/3 = 0.6666666666 the second time, so the errors introduced by treating 2/3 like if it was truly equal to 0.6666666667 should be always the same.

The exemple of Robert Istrael illuminating: same input, three outputs. Why? Are we joking?

In the example I talked about before http://www.mapleprimes.com/forum/sameinputsamecomputerdifferentresults#comment-29414  there aren't sums (so: no possible catastrofic effects), only a singe integral. You could tell to me that an integration is actually a sum (true) but I checked and the integral gives always the same result (same 10 or 11 digits) in the main scope of the program, is only when it is calculated inside an array that it gives different answers every time.

 

Hi

I've found the source of the problem. They are the V's.

I've taken as exemple V31 defined as follows.

 

v31:=proc()
global h1,h2,h3,param_dim, V31;
V31:=Array(1..param_dim,1..param_dim,1..param_dim,1..param_dim, (ii,jj,kk,mm)->
 scalar_pr(h3[ii,jj,kk],h1[mm])
):
end proc:
 

(param_dim is 5)

if I run the program twice I get two different aswers for (say) V31[5,5,5,5]. The same happens if I give v31() twice in the same run.

Now, I've checked that the h3 and h1 arrays have the same values at every run of the program (I tried ten times, always same values).

The same is true for scalar_pr(h31[5,5,5],h1[5]): always the same value at every run. But this (should be) is exactly V31[5,5,5,5]!!

Then the problem must be the Array call inside v31, the costructor for V31.

Resuming: if I call scalar_pr(h31[5,5,5],h1[5]): in the main scope I always get the same result, but inside v31, the 5,5,5,5 component (I focused on this particular component, I guess the same is true for all others) change every run.

I think this is the source of the error in the procedure try() of my second post, as there the V's are summed few hunred of times.

Any idea?

Thanks

S.

Ok

Here I copy a modified version of my input file. I launch it via ssh as was suggested to me in the forum

nohup maple <input.mpl> out.out &

and then I renice it +19. The system is a Linux Fedora 11 64bit, Maple 13.

 

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

parameters:= [Enn,t,phi,Emm,eta]:
param_dim:= nops(parameters):
f_low:=50.0:
f_lso:= (6^(3/2)*Pi*M_tot)^(-1):
f_high:=f_lso:

M_test:=7*10^(-6):

m1:= M_test*1.4 :  #########change it
m2:= M_test*1.4 :   ######### change it

M_tot:=m1+m2:
eta_eff:=(m1*m2)/M_tot^2:
v_lso:=(Pi*M_tot*f_lso)^(1/3):

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

Digits:=15:

p:= proc()
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;
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
,3*eta
,-16*Pi];
Psi:= (2*Pi*f*t-phi-Pi/4+ 3/(128*eta*v^5)*add(alpha_par[i+1]*v^i, i=0..3));
Enn*f^(-7/6)*exp(I*Psi):
 end proc:

Digits:=10:

h1:=Array(1..param_dim, i->
'eval'(
(diff(p(), parameters[i]))
,[Emm=(M_tot*(eta_eff)^(3/5)),eta=eta_eff] )
):

scalar_pr:= proc(A,B)
global variables, parameters,f_low, f_high,S_h,M_tot,eta_eff;
local num;
num:=
simplify(
A*conjugate(B)+conjugate(A)*B
)
;
evalf(
IntegrationTools:-Expand(
2*Int(num/S_h(f), f=f_low..f_high)
)
)
;
#simplify(%);
end proc:

Digits:=10:

AMatrix:=proc()
global low,f_low,f_high,param_dim,M_tot,eta_eff:
local gg,a_1,b_1:
gg:=(a_1,b_1)->
scalar_pr(h1[a_1],h1[b_1])
;
low:= Matrix(param_dim,param_dim, shape=symmetric,gg);
end proc:

AMatrix():

save low, "./low.txt":

Inverse:=proc()
global up,M_tot,param_dim:
local a_3,b_3,gg;
gg:=(a_3,b_3)-> (MatrixInverse(low))[a_3,b_3];
up:= Matrix(param_dim,param_dim,shape=symmetric,gg);
end proc:

Inverse():

h2:=Array(1..param_dim, 1..param_dim, (i,j)->
'eval'(
(diff(p(), parameters[i],parameters[j]))
,[Emm=(M_tot*(eta_eff)^(3/5)),eta=eta_eff])
):
h3:=Array(1..param_dim, 1..param_dim,1..param_dim, (i,j,k)->
'eval'(
(diff(p(), parameters[i],parameters[j],parameters[k]))
,[Emm=(M_tot*(eta_eff)^(3/5)),eta=eta_eff])
):

try:=proc(i)
local temp,MMM,NNN,PPP,QQQ,TTT,ZZZ:
global up,low,V,param_dim;
return
(
+ add(add(add(add(   
up(i,MMM)*up(i,NNN)*up(PPP,QQQ)*
(
+ V4[NNN,MMM,PPP,QQQ]
+ 3*scalar_pr(h2[NNN,QQQ],h2[PPP,MMM])
+ 2*V31[NNN,MMM,PPP,QQQ]
+ V31[MMM,PPP,QQQ,NNN]
)
,MMM=1..param_dim),NNN=1..param_dim),PPP=1..param_dim),QQQ=1..param_dim)
+
add(add(add(add(add(add(   
up(i,MMM)*up(i,NNN)*up(PPP,ZZZ)*up(QQQ,TTT)*
(
+ V3[NNN,PPP,MMM]*V3[QQQ,ZZZ,TTT]
+ 5/2*(-scalar_pr(h2[NNN,PPP],h1[QQQ])*V3[MMM,ZZZ,TTT]-scalar_pr(h2[QQQ,PPP],h1[NNN])*V3[MMM,ZZZ,TTT])
+5/2*( 2*scalar_pr(h2[NNN,QQQ],h1[PPP])*scalar_pr(h2[MMM,ZZZ],h1[TTT]) + scalar_pr(h2[NNN,QQQ],h1[PPP])*scalar_pr(h2[MMM,TTT],h1[ZZZ]) +scalar_pr(h2[NNN,QQQ],h1[PPP])*scalar_pr(h2[TTT,ZZZ],h1[MMM]) )
+2*V21[QQQ,ZZZ,NNN]*V3[MMM,PPP,TTT]
+2* V3[NNN,MMM,TTT]*V21[QQQ,PPP,ZZZ]
+6* V21[NNN,TTT,ZZZ]*V3[MMM,PPP,QQQ]
+V21[MMM,TTT,NNN]*V3[PPP,QQQ,ZZZ]
+2*V21[NNN,QQQ,ZZZ]*V21[PPP,TTT,MMM]
+2*V21[NNN,QQQ,MMM]*V21[PPP,TTT,ZZZ]
)
,MMM=1..param_dim),NNN=1..param_dim),PPP=1..param_dim),QQQ=1..param_dim),TTT=1..param_dim),ZZZ=1..param_dim)
)
end proc:


prova2:= try(2):
prova2:= simplify(Re(try2));#assuming Enn::positive;
save prova2, "./try2.txt":
 

 

where the v's are array of this kind

v3:=proc()
global h1,h2,h3,param_dim, V3;
V3:= Array(1..param_dim,1..param_dim,1..param_dim, (ii,jj,kk)->  
- scalar_pr(h2[ii,jj],h1[kk])
):
end proc:

v3():

and S_h is a rel valued function. (I made changes with respect to the true code, for privacy reasons, so if something is unclear please ask)

 

I launched the imput 5 times, in the same machine, and the file try2.txt contained always a different value:

prova2 := .2050331085e-93/Enn^4;
prova2 := .7892899913e-93/Enn^4;
prova2 := .3927438368e-93/Enn^4;
prova2 := -.2472605212e-93/Enn^4;
prova2 := .6520808383e-93/Enn^4;
prova2 := -.1950615287e-93/Enn^4;

This means that my work of the last three weeks is garbage. I'd love at least to discover why.

Thanks

S.

 

Hi, with respect to the integration we talked about, I've this integral:

 

evalf(IntegrationTools:-Expand(2*Int((h1[5]*conjugate(h1[4]))+(h1[4]*conjugate(h1[5])), f=f_low..f_high))) :
 

where the range of f is real. So, the integrand is of the kind

A B* + A* B

(asterix means complex conjugation), and this, as everybody knows, should be real, but I obtain a complex number ( Re + I Im) as an answer.

Could you help me with that?

Thaks

S.V.

Ok thanks, with the floating range values and evalf it works.

Now I've something more complicated, for which if I use the method you suggested I don't get a results, but a copy of the input. The integral is

Int((2.50*10^6+27.*f^2)*f^(57/50)*(21.+221.*M^(2/3)*f^(2/3))/(2.27*10^18+2.45*10^13*f^2-1.16*10^14*f^(107/50)+5.43*10^10*f^(207/50)-1.20*10^6*f^(307/50)+13.*f^(407/50)), f = 10. .. 2670.)

where M is a positive constant, that maple could put out from the integral. So that if we break the numerator in a sum of therms, each integrale should be of the same kind as before (post #1, the denominator is the same). So I guess maple should be able to give me a numerical answers (with M as constant unknown value), but it doesn't. Where am I wrong?

 

Thanks

S.V.

After some attempts I discovered that are indeed the terms with M in the numerator that maple can't integrate.

Why don't maple just put out M from the integral and give me a results on the form of

M times a number

?