## 65 Reputation

2 years, 221 days

## @acer Thank you so much!!...

@acer Thank you so much!!

## @acer How to convert it to a larger...

@acer How to convert it to a larger with great structure insight and symmetry??  I was thinking about the  converting sin to exp to see structure with Maple. I had a hard time to make it with Maple. Could you give me advice on this??

## Another question...

Why does subs command not work inside  proj_wigner?

 > restart;
 > with(LinearAlgebra):with(ArrayTools):
 >
 >
 > #Oscillator eigenfunction n stasrts 0 psi := (m1,n,x)->2^(-n/2)*(m1)^(1/4)*HermiteH(n,sqrt(m1)*x)*exp(-m1*x^2/2)/(Pi)^(1/4)/sqrt(n!); #Free particle eigenfunction n stasrts 1 phi:= (L,n,x)->sqrt(2/L)*sin(n*Pi*x/L+n*Pi/2);   #KroneckerDelta delta := (m,n)->piecewise(m=n,1,m<>n,0); TripleProduct := proc(A,B,C) local M:   M := MatrixMatrixMultiply(B,C):   MatrixMatrixMultiply(A,M): end proc:
 (1)
 >
 >
 > #Wignerfucntion integral eq1 := W[m,n](q,p) = simplify(1/Pi*Int(phi(L,m,q+y)*exp(-2*I*p*y)*phi(L,n,q-y),y=-L/2+abs(q)..L/2-abs(q))); eq2 := simplify(convert(eq1,int)) assuming(n,integer,m,integer): wp := simplify(eq2) assuming(q>0); wn := simplify(eq2) assuming(q<0);
 (2)
 > # Wigner Matrix for Projectile proj_wigner := proc(d2,L,z_sol) local w:   st[0] := time[real]():   [seq([seq(subs(m=i,n=j,wp),j=1..d2)],i=1..d2)];     (*f1:= (m,n,L)->wp:   w1:= (m,n)->evalc(f1(m,n,L)):   M1:= map(x->unapply(x,p,q),Matrix(d2,w1));   f2:= (m,n,L)->wn:   w2:= (m,n)->evalc(f2(m,n,L)):   M2:= map(x->unapply(x,p,q),Matrix(d2,w2));   Times := time[real]()-st[0]:   lprint(cat("Calculated Wigner matrix1: ",timeparse(Times)));   st[1] := time[real]():   z_sol_dag := map(HermitianTranspose,z_sol):   z_sol_tranp := map(Transpose,z_sol):   z_sol_conj := map(Transpose,z_sol_dag):   WP := zip((x,y)->Trace(TripleProduct(x,y,M1)),z_sol_tranp,z_sol_conj):   WN := zip((x,y)->Trace(TripleProduct(x,y,M2)),z_sol_tranp,z_sol_conj):*)   Times := time[real]()-st[1]:   lprint(cat("Calculated Wigner matrix2: ",timeparse(Times)));   WP,WN: end proc:
 > M := convert(RandomArray(5,distribution=normal),Matrix):
 > proj_wigner(2,1,M): A:= [seq([seq(subs(m=i,n=j,wp),j=1..2)],i=1..2)]:
 ("Calculated Wigner matrix2: ") || (timeparse(3235.776-st[1]))
 >
 >

## Simplest expression...

@acer Actually I want the simplest form of the eq12 and eq13 or at least I want for them to look similar. W_n,m is my matrix component which would be multiplied by another matrix. W_n,m is also function of p  and q. Because of the long expression, I can not see output of the result of W matrix. My goal is to make two equations as simple as possible.

## How to delete question ?...

I posted the stupid question and code. Is there any way to delete this question and code?

## elif...

@acer How can I use if-statement inside elif? I want to use elif rather than if.

## I intended e^(z^2)....

@acer
My original integration was as like as the below code: Because the calculation took a bit longer time, I tried to calculate intermediate steps of the integration which are repeated in real part and imaginary part using changing variables. Thanks for your help in advance.

 > restart; with(IntegrationTools);
 (1)
 > # too slow L:=1,sigma:=0.01,x0:=0.2,k0:=-100},
 > ICs_cal := proc(n_max) odd_eigen := cos(Pi*n*x/L); even_eigen := sin(Pi*n*x/L); g := exp(-(x-x0)^2/2/sigma^2); q[1] := even_eigen*g*cos(k0*x); #even_ICs_real q[2] := even_eigen*g*sin(k0*x): #even_ICs_imag q[3] := odd_eigen*g*cos(k0*x): #odd_ICs_real q[4] := odd_eigen*g*sin(k0*x): #odd_ICs_imag q1 := convert(q,list): q2 := map(f->unapply(f,x),q1); q3 := map(f->Int(f(x),x=-L/2..L/2),q2): q4 := map(f->Change(f,x=sqrt(2)*sigma*z+x0),q3): q5 := map(f->Expand(f),q4): Indets :=  convert(map(f->op(select(has,indets(f,function),Int)),q5),set); q6 := convert(Indets,list); q7 := map(f->unapply(f,n),q6); (*for j from 1 to n_max do  if(is(j,even)) then q8 := map(f->evalf[8](f(j)),q7);                        q9 := seq(q6[i]=q8[i],i=1..4);                      v[j] := evalf[8](unapply(subs(q9,q5[1]),n)(j))+I*evalf[8](unapply(subs(q9,q5[2]),n)(j));         else q8 := map(f->evalf[8](f(j)),q7);          q9 := seq(q6[i]=q8[i],i=1..4);        v[j] := evalf[8](unapply(subs(q9,q5[3]),n)(j))+I*evalf[8](unapply(subs(q9,q5[4]),n)(j));       end if;  convert(convert(v,list),Vector); end do;*) end proc:
 > ICs_cal(20);
 (2)
 >

## @Mariusz Iwaniuk  Thank you! Howe...

Thank you!

However, I would prefer to obtain the result through Maple.

## Thank you so much!...

@Rouben Rostamian  Thank you so much for your detailed and kind explanation in spite of my poor description on the problem due to my poor English. The next step I want to do is to find the density matrices if ths solution Z(t) is obtained.

Any further advice would be very appreciated . Thank you very much again!

## `[Length of output exceeds limit of 1000...

@Rouben Rostamian   Could  you help me one more time?

 > restart;
 > with(IntegrationTools):with(Physics):
 >
 >
 > m1:=1:m2:=1:l:=1:
 > h1 := (m,n)->(n+1/2)*KroneckerDelta[n,m];
 (1)
 > h2 := (mu,nu)->-(nu*Pi/l)^2/(2*m2)*KroneckerDelta[mu,nu];
 (2)
 > v1 := (m,n)->sqrt(min(n,m)!/max(n,m)!)*(2*m1)^(-abs(n-m)/2)*exp(-1/(4*m1))*LaguerreL(min(n,m),abs(n-m),-1/(2*m1));
 (3)
 > v2 := (mu,nu)->4*Pi^2*l*mu*nu*(exp(l/2)-(-1)^(mu+nu)*exp(-l/2))/((Pi*(mu+nu))^2+l^2)/((Pi*(mu-nu))^2+l^2);
 (4)
 > h:=(m,n,mu,nu)->h1(m,n)+h2(mu,nu)+v1(m,n)+v2(mu,nu);
 (5)
 > N:=2:M:=2:
 > Z:= Array(0..N-1,1..M, (i,j)->z[i,j](t));
 (6)
 > the_rhs := Array(0..N-1,1..M,mnu_entry):
 > DE := I*diff(Z,t)= the_rhs;
 (7)
 > IC := eval(Z,t=0) = Array(0..N-1,1..M, (i,j)->`if`(i=j-1,1,0));
 (8)
 > dsol := dsolve({DE,IC});
 (9)

## Thank you a lot!!...

@Rouben Rostamian  Thank you a lot!!

## Z(t) is matrix...

dz:= diff(Z(t),t)= -I* (Multiply(hm1,Z(t))+Multiply(Z(t),hm2)+Multiply(Multiply(vm1,Z(t)),Transpose(vm2)));

## How can I solve this one?...

 > restart;
 > with(IntegrationTools):with(Physics):with(LinearAlgebra):
 >
 >
 > m1:=1:m2:=1:l:=1:
 > h1 := (m,n)->(n+1/2)*KroneckerDelta[n,m];
 (1)
 > h2 := (mu,nu)->-(nu*Pi/l)^2/(2*m2)*KroneckerDelta[mu,nu];
 (2)
 > v1 := (m,n)->sqrt(min(n,m)!/max(n,m)!)*(2*m1)^(-abs(n-m)/2)*exp(-1/(4*m1))*LaguerreL(min(n,m),abs(n-m),-1/(2*m1));
 (3)
 > v2 := (mu,nu)->4*Pi^2*l*mu*nu*(exp(l/2)-(-1)^(mu+nu)*exp(-l/2))/((Pi*(mu+nu))^2+l^2)/((Pi*(mu-nu))^2+l^2);
 (4)
 > #h:=(m,n,mu,nu,m1,m2,l)->evalf(h1(m,n)+h2(mu,nu,m2,l)+v1(m,n,m1)+v2(mu,nu,l));
 >
 > #H:= (m,n,mu,nu)->h(m,n,mu,nu,m1,m2,l);
 > #H(1,1,0,0);
 > #N:=5:M:=5:
 > #eq1:= (m,nu)->diff(z(m,nu,t),t)=-I*Sum(Sum(H(m,n,mu,nu)*z(n,mu,t),n=1..N),mu=0..M);
 > hm1:= Matrix(5,h1);
 (5)
 > hm2:=evalf(Matrix(5,h2));
 (6)
 > vm1:=evalf(Matrix(5,v1));
 (7)
 > vm2:=evalf(Matrix(5,v2));
 (8)
 >
 > dz:= diff(Z(t),t)=Multiply(hm1,Z(t))+Multiply(Z(t),hm2)+Multiply(Multiply(vm1,Z(t)),Transpose(vm2));
 (9)
 > zint:= (m,n)->KroneckerDelta[m,n];
 (10)
 > Z0:=Matrix(5,zint);
 (11)
 > sol:=dsolve({dz,Z(0)=Z0},Z(t));
 >