## 30 Reputation

10 years, 283 days

## Nature of b1, b2, ..., bn...

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

## Can be solved by passing integers instea...

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.

## Good catch on error in question...

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

## f2 takes a variable number of arguments...

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.

## Thanks!...

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

## Yes...

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

## Thanks!...

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!

## Thanks!...

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
﻿