Thanks. Now im struggling with another problem. I want to create a fortran subroutine that gives me a series of pure numerical matrix that i compute  in matrix. i would like to get something like

subroutine EXAMPLE (name1, name2,name3,name4,name5)
                doubleprecision x
                doubleprecision name1(10,10)
                doubleprecision name1(10,10)
                doubleprecision name1(10,10)
                doubleprecision name1(10,10)
                doubleprecision name1(10,10)
                name1(1,1)  = 1.D0
                name1(2,1)  = 1.D0

                .....
                name2(1,1)  = 1.D0

                ......

                ......

              end

I need this since i have to create 10 and maybe more matrices in maple that are different any time i change some parameters. So i want to create only one subroutine (corresponding to a specific set of parameters) and i want  to past and copy  in fortran each of the subroutines (maybe i will have twenty subroutine). (note every element of the matrix comes form an analytical integration of polynomials so i can not just put the parameters as argument o the subroutine but i have to compute analitically each matrix in maple and create a fortran subroutine). You understand that one thing is to copy and paste 20 subroutines, another is to copy and paste 20x10 or more matrices. and if i do something wrong i have to do it again.

restart;
 with(LinearAlgebra):
 with(codegen);
 with(CodeGeneration);
 Kt11inv := Array([[41/4, -7/4, -7/4, -3], [-7/4, 41/4, -7/4, -3], [-7/4, -7/4, 41/4, -3], [3, 3, 3, 0]]);
 Kt11invMASS11 := Array([[-17/40, 16/45, -7/120, 16/45, -2/45, -7/120], [-7/120, 16/45, -17/40, -2/45, 16/45, -7/120], [-7/120, -2/45, -7/120, 16/45, 16/45, -17/40], [-23/90, 13/45, -23/90, 13/45, 13/45, -23/90]]);
 Kt11invKt10 := Array([[-1, 0, 0, 0, 0, 0], [0, 0, -1, 0, 0, 0], [0, 0, 0, 0, 0, -1], [1/9, -4/9, 1/9, -4/9, -4/9, 1/9]]);

f:=makeproc(Kt11invKt10):
g:=makeproc(Kt11invMASS11);
 h:=makeproc(Kt11inv);
 joinprocs([f,g],result_type=array);
codegen[prep2trans](%);
FORT:= codegen[split](%);
Fortran(FORT,optimize,deducetypes=false,defaulttype=numeric,limitvariablelength=false);

where Kt11inv,Kt11invMASS11, Kt11invKt10 are three of my matrices. Like that it work partially, two main problems:

1)it creates a lot of extra stuff in the subroutine

2)if i use more then two matrix it does not work. i tried with joinprocs([f,g,h],result_type=array); and it gives me an error.

any suggestions?

thanks so much

alberto

no i mean I want basex (or another  name that i can provide) instead of cgret. In your example you still have cgret! i.e.:

subroutine basex (x, basex)
doubleprecision x
doubleprecision basex(10)
doubleprecision t14
doubleprecision t2
doubleprecision t4
doubleprecision t6
t2 = x ** 2
t4 = t2 * x
t6 = t2 ** 2
t14 = t6 ** 2
basex(1) = 0
basex(2) = x - 0.1D1
basex(3) = t2 - 0.1D1
basex(4) = t4 - 0.1D1
basex(5) = t6 - 0.1D1
basex(6) = t6 * x - 0.1D1
basex(7) = t6 * t2 - 0.1D1
basex(8) = t6 * t4 - 0.1D1
basex(9) = t14 - 0.1D1
basex(10) = t14 * x - 0.1D1
end

(now i did it pasting and copying them one by one...)

 I have to do kind of one hudreds of this souboutine for different polynomials having different names i dont want to  do replace every time that i need it, i want to have the final soubroutine with the name that i want.

thanks

alberto

You know what's the problem with using  basex:=unapply(<seq(x^i-1,i=0..9)>,x) as you suggested ?

I'm not able to prescribe the name of the vector of elements in the fortran output. i uses codegen[makeproc] since i want the declaration of all the variables t1,t2,t3 that are created after optimization like

    doubleprecision t1
      doubleprecision t2
      doubleprecision t3

infact sometimes they are tens of variables so its long to copy and past all of them! but now im not able to provide a different name for the vector, Fortran  provides cgret by default!!

do you know how to change it? this is the code i wrote that provide cgret by default

restart;
 with(LinearAlgebra):
 basex:=unapply(<seq(x^i-1,i=0..9)>,x);
 aaa:=basex(x):
 codegen[makeproc](aaa);codegen[optimize](%); 
 codegen[prep2trans](%);
 FORT:= codegen[split](%);
 CodeGeneration:-Fortran(FORT, optimize,defaulttype=numeric,deducetypes=false,resultname=THETA);

thanks

alberto
 

ooo OK MAPLE 12 stuff...im still using 11

thanks

alberto

thanks.

1)sorry what do you mean for: "rtables, including Matrices and Vectors, use () for INDEXES now" ?why indexes? if i write bases (2), it computes bases using x=2, it is not the index I don't understand.

2) and if I want for example(it happens often to me) to use if statements to exclude any i, for example i=3, can i do it using seq in a compact way?or should i do something like this:

restart;
 basex:=Vector(10);
 for i from 0 to 9 do
 if i<>3 then
 basex[i+1]:=x^i-1;
 else
 basex[i+1]:=x^i-9;
 fi:
 od:
 basex;basex:=unapply(basex,x);
 CodeGeneration:-Fortran(basex(x), optimize,
     deducetypes=false,defaulttype=numeric,resultname=THETA);

thank you very much

alberto

Thanks. You are right about set and list. It's a small but important difference,
And why do you think linalg is deprecated?what the problem with using vector?I know that know there is LinearAlgebra, but for the seek of knowledge why is linalg deprecated?

thanks a lot

alberto

thanks. a lot.

1) you are right. thanks

2)i tried with your suggestion

sol:= solve({seq(eq[i], i = 1..np)}, {seq(a[i], i = 1..np)}):

and it works. i guess there is a bug with solve and convert/set since when i use :

Thanks i opted to use the  first hack (using eval) (NOTE: it work only writting b=Q not b=Q(x,y,t) ). It slows down a little bit the code(because he computes the derivatives of signum(1,Q) and all the stuff that moltiplies it,and the code was already slow) but that's ok. You suggested me also to normalize, how do you suggest me to do?using "normal"?because my aim is to have the minor computational cost when  I create the fortran subroutine and i used simplify(%,size). Do you think in general is better to use also normal?before or after simplify(%,size)? I thought simplify(%,size) was the best and included the normalization).

Thanks

alberto

I attach you the file with the correction you suggest me if you need it.

Download 8216_cauchyKOVALEVSKIshallowFONDOmobPETER-MULLERfriction_ord3.mws
View file details

thanks I should use 3/2 you are right.It seems working better but now i have another problem with assume: I have declared Q real and <>0 because I want that signum(1,Q) is zero (i remember you: my code just derives a few times the differential equations to get the high order derivatives),otherwise when it derives it keeps all the term that multiply signum(1,Q) and it become very complicated. But I would like also to declare (as your suggestion) all the other variables like:
 

assume(Q::real,Q<>0,H::real,B::real,H-B>0,S::real,S>0,AA::real,Ks::real,Ks>0,Diam::real,Diam>0,DELTA::real,DELTA>0,TScr::real,TScr>0,qs::real);
 

, but if i declare them I don't know why there is a conflict and signum(1,Q) is not zero,it is zero only if I declare only assume(Q::real,Q<>0). Maybe because Q,H,S,B are function of x,y,t and the assumption goes on x,y,t, but it's  a weird behaviour. Also if i use the piecewise function i still have the same problem.How can i declare all the variables and get signum(1,Q)equal to zero?

thanks
 

ps i attached the file, you can check if signum(1,Q) is zero or not when it computes :
 

eq[3,1]:= diff(S,t) = - diff(S*Q/(H-B),x);

just at the beginning of the sheet. Thanks.
 

Download 8216_cauchyKOVALEVSKIshallowFONDOmobPETER-MULLERfriction_ord3.mws
View file details
 

For Pagan: of course I meant  Tcr>0 sorry.

For jakubi: First of all thanks. for the qs expression it can be what i want   but it is not general! I have to use it inside a for cycle where i don't know a priori the results  and the powers change every time for different inputs and they are not always 3/2 and 71/5 so I can't use this method. is there no way to collect products of powers having equal basis in a general way?

thanks a lot

Alberto

So there is no way to collect products of powers having equal basis? I mean is it possible to collect "If" in the first equation:

eq2:=8*(-(-If*(q*ks^4/If^(1/2))^(3/5)+Tcr*Delta*Ds*ks^3)/Delta/Ds/ks^3)^(3/2) = qs/(Delta*GRAVIT*Ds^3)^(1/2);

and then solve?

because now i have  a similar problem, i got this expression after a few calculation and I have to derive it in H and Q and B but I want before to collect the term 1/(H-B)^(7/3)*abs(q)^2/ks^2/Delta/Ds:

qs := c1*(1/(H-B)^(7/3)*abs(q)^2/ks^2/Delta/Ds)^(3/2)*(1/(H-B)^(7/3)*abs(q)^2/ks^2/Delta/Ds/tr)^(71/5)*signum(q)*(Delta*GRAVIT*Ds^3)^(1/2)

if I don't collect the term obviusly  the derivative becomes too complicate.

thanks

alberto

Thanks a lot.

In this case:

> restart;
> eq2:=8*(-(-If*(q*ks^4/If^(1/2))^(3/5)+Tcr*Delta*Ds*ks^3)/Delta/Ds/ks^3)^(3/2) = qs/(Delta*GRAVIT*Ds^3)^(1/2);
> assume(If>0,ks>0,q>0,Tcr>o,Delta>0,Ds>0,qs>0,GRAVIT>0);;
> use RealDomain in If:=solve(eq2,If) end use;
 

Do you know how can i collect the powers?I tried with combine(%,power); and simplify(%,size); but it doesnt work. I would like to get an output like this(gained hand-collecting the variable "If" before solving).

>  restart:

>  eq2:=(If^(7/10)*(q*ks^4)^(3/5)) = ((qs/8/sqrt(Delta*GRAVIT*Ds^3))^(2/3)+Tcr)*Delta*Ds*ks^3;

> If:= solve(eq2,If);

thanks a lot

alberto

I had already  tried in this way but its not correct because all the derivative of Q are zero now!! (for example diff(Q,x); gives zero!) . why?

Thanks:

How can i use assume together with declare to get the compact form in the output?i wrote this:

restart:
 printlevel := 1;with(PDEtools):
with(codegen,fortran):
with(CodeGeneration,Fortran):
alias(Q = Q(x, y, t), H = H(x, y, t), S = S(x,y,t) , B = B(x, y, t)):
declare(Q(x, y, t), H(x, y, t), S(x,y,t), B(x, y,t)): 
assume (Q::real,Q<>0);
 interface(showassumed=2);
qs:=aa/(1-p)*abs(Q)^mm/(H-B)^mm*abs(1,Q);
diff(H, t) = -(diff(Q, x))-diff(qs,x);

i get the simplifications that i wanted but not the compact form in the output

thanks

alberto
 

I have found the way for a fast and immediate substistution. Let my know anyway if the option will be available int the future.