charlesforgy

30 Reputation

2 Badges

6 years, 61 days

MaplePrimes Activity


These are replies submitted by charlesforgy

@Carl Love Thank you for your feedback.

 

The bn's can take any of the following values:
ax, ay, az bx, b_y, bz, cx, cy, cz, dx, dy, dz, or 0

(b_y used to avoid naming a variable by)

They should all be local to each procedure that uses them, and are not created by cat or the like.  However, inside f3 they sometimes are reassigned as the program changes combinations.  

 

If it helps, here is the relevant code, f1 is L_gen, f2 is F5_gen, and f3 is actually both T_first_deriv_2 and T_sec_deriv_2.

 

L_gen:=proc(ai, aj, ak, bi, bj, bk, ci, cj, ck, di, dj, dk)

#takes angular momenta numbers and returns list for F5_gen, F6_gen, F7_gen, or F8_gen  

local L, ax, ay, az, bx, b_y, bz, cx, cy, cz, dx, dy, dz, n:

if ai+aj+ak+bi+bj+bk+ci+ cj+ ck+ di+dj+ dk=5 then

        L:=[0,0,0,0,0]:

elif ai+aj+ak+ bi+bj+bk+ci+ cj+ ck+ di+dj+ dk=6 then

        L:=[0,0,0,0,0, 0]:

elif ai+aj+ak+ bi+bj+bk+ci+ cj+ ck+ di+dj+ dk=7 then

        L:=[0,0,0,0,0, 0, 0]:

elif ai+aj+ak+bi+bj+bk+ci+ cj+ ck+ di+dj+ dk=8 then

        L:=[0,0,0,0,0, 0, 0, 0]:

fi:

 

 

 

 n:=1;

if ai=1then

        L[n]:=1; n:=n+1;

elif ai =2 then

        L[n]:=1; L[n+1]:=1 ; n:=n+2;

fi;

if aj=1then

        L[n]:=2; n:=n+1;

elif aj =2then

        L[n]:=2; L[n+1]:=2; n:=n+2;

fi;

if ak=1then

        L[n]:=3; n:=n+1;

elif ak =2then

        L[n]:=3; L[n+1]:=3; n:=n+2;

fi;

if bi=1then

        L[n]:=4; n:=n+1;

elif bi=2then

        L[n]:=4; L[n+1]:=4; n:=n+2;

fi;

 ...

...

return L:

end:

 

F5_gen := proc (alpha, beta, gamma, delta, L::list)

local a, b, n, i, j, k, l, L0, L1, V3, V4, V5, V0;

V3 := 0;

V4 := 0;

L0 :=[0,0,0,0,0];

L1 :=[0,0,0];

for i to 4 do

        for j from i+1 to 5 do

                L0 := L;

                L0[i] := 0;

                L0[j] := 0;

                V0 := T_sec_deriv_2(alpha, beta, gamma, delta, L[i], L[j]);

                n := 1;

                for a to 5 do

                        if L0[a] <> 0 then

                                L1[n] := L0[a];

                                n := n+1

                        end if;

                end do;

                b := 3;

                for k to 2 do

                        for l from k+1 to 3 do

                                V3 := V3+V0*T_sec_deriv_2(alpha, beta, gamma, delta, L1[k], L1[l])*T_first_deriv_2(alpha, beta, gamma, delta, L1[b]);

                                b := b-1

                        end do;

                end do;

        end do;

enddo;

 

for i to 4 do

        for j from i+1 to 5 do

                L0 := L;

                L0[i] := 0;

                L0[j] := 0;

                V4 := V4+T_sec_deriv_2(alpha, beta, gamma, delta, L[i], L[j])*T_first_deriv_2(alpha, beta, gamma, delta, L0[1])*T_first_deriv_2(alpha, beta, gamma, delta, L0[2])*T_first_deriv_2(alpha, beta, gamma, delta, L0[3])*T_first_deriv_2(alpha, beta, gamma, delta, L0[4])*T_first_deriv_2(alpha, beta, gamma, delta, L0[5])

        end do;

enddo;

V5 := T_first_deriv_2(alpha, beta, gamma, delta, L[1])*T_first_deriv_2(alpha, beta, gamma, delta, L[2])*T_first_deriv_2(alpha, beta, gamma, delta, L[3])*T_first_deriv_2(alpha, beta, gamma, delta, L[4])*T_first_deriv_2(alpha, beta, gamma, delta, L[5]);

 

return -((1/2)*F3*V3)+F4*V4-F5*V5:

end proc:

 

T_first_deriv_2:=proc(alpha, beta, gamma, delta,  a)

local T1, Ax, Ay, Az, Ta_b, Tc_d, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, L0;

L0:=[0,0,0,  0,0,0,   0,0,0,   0,0,0]:

if a =1  then

        L0[1]:=1;

elif a =2  then

        L0[2]:=1;

elif a =3  then

        L0[3]:=1;

elif a =4  then

        L0[4]:=1;

 ...

...

if L0[1]=1then

        v1 := T1[1] #The symbolic value of the derivative

else

        v1 := 1

endif;

 

if L0[2]=1then

        v2 := T1[2]

else

        v2 := 1

endif;

 ...

...

return v1*v2*v3*v4*v5*v6*v7*v8*v9*v10*v11*v12:

end proc:

 

T_sec_deriv_2:=proc(alpha, beta, gamma, delta,  a, b)

local L0, T1, T2, Ta_b2, Tc_d2, Tmix;

L0:=[0,0,0,   0,0,0,   0,0,0,   0,0,0];

if a = b then

        if a = 1 then

                L0[1]:=2;

        elif a = 2 then

                L0[2]:=2;

        elif a = 3 then

                L0[3]:=2;

        elif a = 4 then

                L0[4]:=2;

 ...

...

else

        if a= 1 or b =1  then

                L0[1]:=1;

        fi;

        if a = 2  or b = 2  then

                L0[2]:=1;

        fi;

        if a = 3  or b = 3  then

                L0[3]:=1;

        fi;

        if a = 4  or b = 4  then

                L0[4]:=1;

        fi;

 ...

...

end if;

if L0[1]+L0[2]+L0[3]+L0[4]+L0[5]+L0[6]+L0[7]+L0[8]+L0[9]+L0[10]+L0[11]+L0[12] <> 2 then

        return"error, not second derivative"

elif L0[1] = 2 or L0[2] = 2 or L0[3] = 2 then

        return T1[1] #The symbolic value of the derivative

elif L0[4] = 2 or L0[5] = 2 or L0[6] = 2 then

        return T1[2]

elif L0[7] = 2 or L0[8] = 2 or L0[9] = 2 then

        return T1[3]

elif L0[10] = 2 or L0[11] = 2 or L0[12] = 2 then

        return T1[4]

elif L0[1]+L0[4] = 2 or L0[2]+L0[5] = 2 or L0[3]+L0[6] = 2 then

        return T2[1]

elif L0[1]+L0[7] = 2 or L0[2]+L0[8] = 2 or L0[3]+L0[9] = 2 then

        return T2[2]

elif L0[1]+L0[10] = 2 or L0[2]+L0[11] = 2 or L0[3]+L0[12] = 2 then

        return T2[3]

elif L0[4]+L0[7] = 2 or L0[5]+L0[8] = 2 or L0[6]+L0[9] = 2 then

        return T2[4]

elif L0[4]+L0[10] = 2 or L0[5]+L0[11] = 2 or L0[6]+L0[12] = 2 then

        return T2[5]

elif L0[7]+L0[10] = 2 or L0[8]+L0[11] = 2 or L0[9]+L0[12] = 2 then

        return T2[6]

else

        return 0

endif:

 

end proc:

 


Clearly there are many more of the same if-then statements that I omitted for brevity, though I can post entire codes if anyone wants to read through them.  There are procedures essentially the same as F5_gen but longer and more complicated for F6_gen, F7_gen, and F8_gen.

Thank you again for all your help.

 

 

I've found a workaround for the problem, though I'm still not entirely certain what causes it.  The problem appears to somehow stem from Maple's method of passing data structures by reference.  

In any event, if--instead of passing variables b1, b2, etc.--we instead define 1=b1, 2=b2, etc., and then have f1 create and pass a list of integers the programs work.

Also, as noted below there is a typo in the original question, f2 calls f3, not itself.

@tomleslie Sorry about that, I made a mistake in how I asked the question above.  Points 1 and 2 answer eachother, the part in which f2 calls itself should have had it calling f3 (that mistake is not present in the code, though).

 

I will try to post the relevant lines of the actual code without adding too much needless complexity to the question--the total code is well over 1,000 lines.

@Stephen Forrest Thanks for your reply.

 

As to (2), there shouldn't be any inplace assignments as you write above.

 

 For (1), f2 currently takes a variable number of arguments, though I could specify the number it takes.

 

As to tracing it, whether I run the program as 

L0:=f1(a1, a2, ..., an):
B:=f2(alpha, beta,...,L0): 

or as 

f2(alpha, beta, ..., [b1, b2, ..., bn])

 

the program trace returns the same results, namely:

{--> enter f2, args = alpha, beta, gamma, delta, [ax, ax, bx, bz, cx]

...

L0:=[ax, ax, bx, bz, cx]

But then when f2 invokes f3, it only returns the correct result in the latter case.

@acer,  @Carl Love  Thanks to both of you for your help.  Both suggestions seem fix the problem.

@Carl Love Yes, I would either call the function using evalhf or evalf, or modify the resulting code to evaluate the result before returning the values.

@Carl Love 

I wouldn't have even thought of changing the digits, which had been set to 15.  However, I played around with values between 14 and 21.  Some values switched the message to "non-fatal error when reading from kernel," but 16 specifically seems to make everything run smoothly.

If you'd still like to see the relevant part of the code, it's:

 

makeSolver := proc (H, d) local dim, d0, vars, dsol, i, j;

dim := LinearAlgebra:-RowDimension(H);

d0 := Array(1 .. dim^2); mat2vec(dim, d, d0);

vars := [seq(seq(dd[i, j], j = 1 .. dim), i = 1 .. dim)];

dsol := dsolve(numeric, number = dim^2, procedure = prop, start = 0, initial = d0, procvars = vars, 'complex' = true, method = rkf45, maxfun = 0);

return eval(dsol);

end:

 

Thanks for the help!

@Carl Love 

I wouldn't have even thought of changing the digits, which had been set to 15.  However, I played around with values between 14 and 21.  Some values switched the message to "non-fatal error when reading from kernel," but 16 specifically seems to make everything run smoothly.

If you'd still like to see the relevant part of the code, it's:

 

makeSolver := proc (H, d) local dim, d0, vars, dsol, i, j;

dim := LinearAlgebra:-RowDimension(H);

d0 := Array(1 .. dim^2); mat2vec(dim, d, d0);

vars := [seq(seq(dd[i, j], j = 1 .. dim), i = 1 .. dim)];

dsol := dsolve(numeric, number = dim^2, procedure = prop, start = 0, initial = d0, procvars = vars, 'complex' = true, method = rkf45, maxfun = 0);

return eval(dsol);

end:

 

Thanks for the help!

Page 1 of 1