Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

I have to compare the probabilty density of a data set of the normalised gaps between eigenfunction energy levels of random matrices to the Wigner surmise Pi/2*n*exp((-Pi)/4*n^2). I have collected the data set, but cannot find a function to give a numerical value for the probability density function to compare. I have tried using KernelDensity but it doesn't seem to work.

Here is the code:


 

restart; with(RandomTools); with(LinearAlgebra); with(Statistics); with(plots)

NULL

First create a procedure for the Gaussian orthogonal ensemble, GOE, with mean 0 and variance 1. This procedure extracts the normalised eigenvalues for a symmetric matrix added to its transpose and divided by 2:

GOE := proc (N::integer) local A, B, C; A := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); B := Transpose(A); C := (1/2)*A+(1/2)*B; return Re(Eigenvalues(C)/sqrt(2*N)) end proc

data := [seq(GOE(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

NULL

The histogram shows a normal distribution centred around 0, as expected. Increasing the value of N makes the distribution trend more and more towards the ideal Gaussian distribution as N tends to infinity, since a bigger sample size tends to the statistical mean. This procedure can be repeated for uniformly distributed random numbers in the range -1 to +1:

GOE2 := proc (N::integer) local A, B, C; A := Matrix(N, N, Generate(distribution(Uniform(-1, 1)), makeproc = true)); B := Transpose(A); C := (1/2)*A+(1/2)*B; return Re(Eigenvalues(C)/sqrt(2*N)) end proc

data := [seq(GOE2(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

This creates a curved distribution steeper and narrower than the above normal distribution, centred around 0 and tailing off below -0.5 and above +0.5. As N tends to infinity, this curve becomes smoother and more rounded, as with the Gaussian distribution.

 

The same principle also applies to complex Hermitian matrices in the Gaussian unitary emsemble (GUE), though complex numbers cannot be normally graphed on histograms. By taking the real and imaginary parts of the eigenvalues separately, it is possible to plot the distribution on a histogram, since both parts have the same distribution.

GUE := proc (N::integer) local A, B, C, D, E; A := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); B := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); C := A+I*B; D := Transpose(C); E := (1/2)*C+(1/2)*D; return Re(Eigenvalues(E)/sqrt(4*N)); return Im(Eigenvalues(E)/sqrt(4*N)) end proc

data := [seq(GUE(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

The result is a similar normal distribution to the GOE, but narrower due to the increased normalisation constant.

 

Finally, a procedure can return the spacings between the nth energy levels of the eigenvalues, taking the nth level spacing as "lambda[n+1]-lambda[n]."This procedure can work for the GOE or GUE distributions, and must check for these as valid arguments.

spacing := proc (n::integer, N::integer, ensemble) local A, B; if ensemble = GOE then A := GOE(N); return A[n+1]-A[n] elif ensemble = GUE then B := GUE(N); return B[n+1]-B[n] else print("Invalid argument. Only GOE and GUE are accepted.") end if end proc

The Wigner surmise of probability density for the GOE and GUE predicts that the probablity density function follows specific formulae for each spacing. This can be compared by taking many results for a given spacing and dividing by the mean, the using the kernel density function with a Gaussian kernel to estimate the probability density of the data set.

data := [seq(spacing(2, 4, GOE), i = 1 .. 500)]

normaliseddata := data/Mean(data)

F := KernelDensity(normaliseddata, bandwidth = 1/4, kernel = 'gaussian', left = -1000, right = 1000, method = exact)

evalf(F(2.0))

0.400007351200000e-304+0.*I

(1)

WignerGOE := proc (n::integer) return (1/2)*Pi*n*exp(-(1/4)*Pi*n^2) end proc

WignerGOE(2)

Pi*exp(-Pi)

(2)

NULL


 

Download matrices.mw
 

restart; with(RandomTools); with(LinearAlgebra); with(Statistics); with(plots)

NULL

First create a procedure for the Gaussian orthogonal ensemble, GOE, with mean 0 and variance 1. This procedure extracts the normalised eigenvalues for a symmetric matrix added to its transpose and divided by 2:

GOE := proc (N::integer) local A, B, C; A := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); B := Transpose(A); C := (1/2)*A+(1/2)*B; return Re(Eigenvalues(C)/sqrt(2*N)) end proc

data := [seq(GOE(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

NULL

The histogram shows a normal distribution centred around 0, as expected. Increasing the value of N makes the distribution trend more and more towards the ideal Gaussian distribution as N tends to infinity, since a bigger sample size tends to the statistical mean. This procedure can be repeated for uniformly distributed random numbers in the range -1 to +1:

GOE2 := proc (N::integer) local A, B, C; A := Matrix(N, N, Generate(distribution(Uniform(-1, 1)), makeproc = true)); B := Transpose(A); C := (1/2)*A+(1/2)*B; return Re(Eigenvalues(C)/sqrt(2*N)) end proc

data := [seq(GOE2(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

This creates a curved distribution steeper and narrower than the above normal distribution, centred around 0 and tailing off below -0.5 and above +0.5. As N tends to infinity, this curve becomes smoother and more rounded, as with the Gaussian distribution.

 

The same principle also applies to complex Hermitian matrices in the Gaussian unitary emsemble (GUE), though complex numbers cannot be normally graphed on histograms. By taking the real and imaginary parts of the eigenvalues separately, it is possible to plot the distribution on a histogram, since both parts have the same distribution.

GUE := proc (N::integer) local A, B, C, D, E; A := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); B := Matrix(N, N, Generate(distribution(Normal(0, 1)), makeproc = true)); C := A+I*B; D := Transpose(C); E := (1/2)*C+(1/2)*D; return Re(Eigenvalues(E)/sqrt(4*N)); return Im(Eigenvalues(E)/sqrt(4*N)) end proc

data := [seq(GUE(3), i = 1 .. 5000)]

Histogram(data, binwidth = 0.5e-1, range = -2 .. 2)

The result is a similar normal distribution to the GOE, but narrower due to the increased normalisation constant.

 

Finally, a procedure can return the spacings between the nth energy levels of the eigenvalues, taking the nth level spacing as "lambda[n+1]-lambda[n]."This procedure can work for the GOE or GUE distributions, and must check for these as valid arguments.

spacing := proc (n::integer, N::integer, ensemble) local A, B; if ensemble = GOE then A := GOE(N); return A[n+1]-A[n] elif ensemble = GUE then B := GUE(N); return B[n+1]-B[n] else print("Invalid argument. Only GOE and GUE are accepted.") end if end proc

The Wigner surmise of probability density for the GOE and GUE predicts that the probablity density function follows specific formulae for each spacing. This can be compared by taking many results for a given spacing and dividing by the mean, the using the kernel density function with a Gaussian kernel to estimate the probability density of the data set.

data := [seq(spacing(2, 4, GOE), i = 1 .. 500)]

normaliseddata := data/Mean(data)

F := KernelDensity(normaliseddata, bandwidth = 1/4, kernel = 'gaussian', left = -1000, right = 1000, method = exact)

evalf(F(2.0))

0.400007351200000e-304+0.*I

(1)

WignerGOE := proc (n::integer) return (1/2)*Pi*n*exp(-(1/4)*Pi*n^2) end proc

WignerGOE(2)

Pi*exp(-Pi)

(2)

NULL


 

Download matrices.mw

 

We have just released updates to Maple and MapleSim.

Maple 2020.2 includes corrections and improvements to printing and export to PDF, support for macOS 11.0, more MATLAB connectivity, resolves issues with the installation of the Maplesoft Physics Updates, and more.  We recommend that all Maple 2020 users install these updates.

This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2020.2 download page, where you can also find more details.

If you are also a MapleSim user, this Maple update will be installed automatically when you update your MapleSim installation to the newly released MapleSim 2020.2. The MapleSim update also includes many updates to MapleSim, the MapleSim CAD Toolbox, specialized MapleSim libraries, MapleSim connectivity tools, and MapleSim Insight.  You will find details about new features and improvements, as well as instructions on obtaining the update, on the MapleSim 2020.2 page.

I have an assignment for Q that after subsequent other assignments and substitutions  results in 

-XC1*R1^2/((R1^2 + XC1^2)*(R1*XC1^2/(R1^2 + XC1^2) + 12960.54302))

 

when I type Q.

 

I would like to solve this for XC1, for values of Q that make  XC1 is real. 

How do I do this?  Can I rearrange this assignment?   

I guess I could do something like this:

eq1:= -XC1*R1^2/((R1^2 + XC1^2)*(R1*XC1^2/(R1^2 + XC1^2) + 12960.54302))

solve(eq1=Q,XC1)

but Q as a function of XC1 is a derived from other relationships.

 

The worksheet probably makes what I'm asking more clear.    I was able to get the result, but I'm sure there is a better, more elegant  way to do what I needed to do...

 

Thanks
 

``

restart

Series-Parallel Conversion Equations as a function of Q

• 

Q = Xs/Rs = Rp/Xp;

• 

Rp := Rs*(Q^2 + 1);

``

``

Impedance Transpormation Equations

Rp := proc (Rs, Xs) options operator, arrow; (Rs^2+Xs^2)/Rs end proc

Xp := proc (Rs, Xs) options operator, arrow; (Rs^2+Xs^2)/Xs end proc

Rs := proc (Rp, Xp) options operator, arrow; Rp*Xp^2/(Rp^2+Xp^2) end proc

Xs := proc (Rp, Xp) options operator, arrow; Xp*Rp^2/(Rp^2+Xp^2) end proc

``

zL := 1.343+I*131.925

1.343+131.925*I

(1)

``

XL0p := Xp(Re(zL), Im(zL))

131.9386717

(2)

RLp := Rp(Re(zL), Im(zL))

12960.54302

(3)

QLp := RLp/XLp

12960.54302/XLp

(4)

XC2 := -XL0p

-131.9386717

(5)

Q := XL1/(R1s+RLp)

XL1/(R1s+12960.54302)

(6)

R1s := Rs(R1, XC1)

R1*XC1^2/(R1^2+XC1^2)

(7)

XC1s := Xs(R1, XC1)

XC1*R1^2/(R1^2+XC1^2)

(8)

XL1 := -XC1s

-XC1*R1^2/(R1^2+XC1^2)

(9)

Q

-XC1*R1^2/((R1^2+XC1^2)*(R1*XC1^2/(R1^2+XC1^2)+12960.54302))

(10)

``

tmp1 := solve(-XC1*R1^2/((R1^2+XC1^2)*(R1*XC1^2/(R1^2+XC1^2)+12960.54302)) = Tmp, XC1)

(-25000.*R1+(-0.3240135755e14*R1*Tmp^2+625000000.*R1^2-0.4199391884e18*Tmp^2)^(1/2))*R1/(Tmp*(50000.*R1+648027151.)), -1.*(25000.*R1+(-0.3240135755e14*R1*Tmp^2+625000000.*R1^2-0.4199391884e18*Tmp^2)^(1/2))*R1/(Tmp*(50000.*R1+648027151.))

(11)

"solve(indets(numer(tmp1[1]))"

{R1, Tmp, (-0.3240135755e14*R1*Tmp^2+625000000.*R1^2-0.4199391884e18*Tmp^2)^(1/2)}

(12)

simplify(solve(op(3, indets(numer(tmp1[1])))^2 > 0, Tmp, parametric))

piecewise(R1 <= -8398783768000/648027151, [[Tmp = Tmp]], R1 < 0, [[50*5^(1/2)*R1/(648027151*R1+8398783768000)^(1/2) < Tmp, Tmp < -50*5^(1/2)*R1/(648027151*R1+8398783768000)^(1/2)]], R1 = 0, [], 0 < R1, [[-50*5^(1/2)*R1/(648027151*R1+8398783768000)^(1/2) < Tmp, Tmp < 50*5^(1/2)*R1/(648027151*R1+8398783768000)^(1/2)]])

(13)

``


 

Download pi-filter_anal_copy.mw

 

 

EDIT TIME: 14:30 CET

where P(1) and P(2) are NxN matrix functions.

 

My trial code for STEP 0

 alpha times integral w.r.t. x

 Int_x__alpha:=proc(term,alpha): 
 return
select(has,term,x).P(alpha)^T.remove(has,term,x)
end proc: 

 alpha times integral w.r.t. t

Int_t__alpha:=proc(term,alpha):   #alpha times integral w.r.t. t
 return
remove(has,term,t).P(alpha).select(has,term,t)
end proc:

when I run the last procedure for the testing

 Int_t__alpha(Psi(x)^T.C. Psi(t),2);

I get

But it must be 

Because the multiplication is not commutative in Matrices. So, the last procedure must be corrected. 

DETAILS for the procedures:

-------------------------------------------------------------------------------------------------------------------------------------

MAIN QUESTION:

Suppose that we have a PDE as follows 

                                             ...(3)

 

subject to appropriate Initial and Boundary conditions.

-------------------------------------------------------------------------------------------------------------------------------------

STEP 1

  • Find the highest derivative w.r.t. x and w.r.t. t. Then, Let the trial function be the summation derivative of these highest derivatives. I mean

Trial Function:                                           ...(4)

where Psi(x), Psi(t) are Nx1 vectors and C is a NxN matrix.

I can't write a maple code for selecting the trial function.  May be you can.
trial_function:=diff(u(x,t),x,t)=Psi(x)^T.C. Psi(t); 

# I deliberately used ^T instead of ^+ for Transpose.
# If I use ^+, the transpose sign doesn't appear in 2d output. May be you have an other idea.

 

-------------------------------------------------------------------------------------------------------------------------------------

STEP 2

  • Integrate the Eq.4  w.r.t. t from 0 to t, we have
STEP2:=int(  lhs(trial_function) ,t=0..t)=Int_t__alpha(rhs(trial_function),1);


The code must be improved. Firstly, substitute t=s in lhs(trial_function) and then integrate s=0..t 

-------------------------------------------------------------------------------------------------------------------------------------

STEP 3

 

int(lhs(STEP2),x=0..x)= Int_x__alpha(rhs(STEP2),1);

The code must be improved.

-------------------------------------------------------------------------------------------------------------------------------------
STEP 4

  • Integrate Eq.4 w.r.t x

-------------------------------------------------------------------------------------------------------------------------------------
STEP 5

  • Substitute Eq. (5), Eq. (6), Eq. (7) to Eq. (3),
  • I mean substituting  u_x(x,t), u_t(x,t), u(x,t) to PDE.

-------------------------------------------------------------------------------------------------------------------------------------

STEP 6

DOWNLOAD ALL MAPLE CODE: all_code.mw

Any ideas why the attached package won't install?

We've been using Maple Cloud previously, but that doesn't work anymore apparently (same problem as stated here, https://www.mapleprimes.com/questions/226154-Installing-Packages-From-The-Cloud-Hangs, support case filed previously).

But I can't even get it installed manually either.

NODELibrary.zip

how I can sort differential equation [q] in terms of u(x,y,z,t) and its derivatives.

In other words, we should find three relations including L11(u) ,L12 (v), and L13(w). 

For instance in L13{w) only w(x,y,z,t) and its derivatives would appear...

For example;

sort.mw


 

q := h*(A11*(diff(u(x, y, z, t), x, x)-(diff(w(x, y, z, t), x))/R+(diff(w(x, y, z, t), x))*(diff(w(x, y, z, t), x, x)))+A12*(diff(v(x, y, z, t), x, y)-(diff(w(x, y, z, t), x))/a+(diff(w(x, y, z, t), y))*(diff(w(x, y, z, t), x, y))))+h^2*(-B11*(diff(w(x, y, z, t), x, x, x))+B12*(-(diff(w(x, y, z, t), x, y, y))-(diff(v(x, y, z, t), x, y))/a))+B66*h^2*(-2*(diff(w(x, y, z, t), x, y, y))-(diff(v(x, y, z, t), x, y))/a)+A66*h*(diff(u(x, y, z, t), y, y)+diff(v(x, y, z, t), x, y)+(diff(w(x, y, z, t), y, y))*(diff(w(x, y, z, t), x))+(diff(w(x, y, z, t), y))*(diff(w(x, y, z, t), x, y))) = rho*(diff(u(x, y, z, t), t, t))-e0^2*a^2*(rho*(diff(u(x, y, z, t), t, t, x, x))+rho*(diff(u(x, y, z, t), t, t, y, y)))

h*(A11*(diff(diff(u(x, y, z, t), x), x)-(diff(w(x, y, z, t), x))/R+(diff(w(x, y, z, t), x))*(diff(diff(w(x, y, z, t), x), x)))+A12*(diff(diff(v(x, y, z, t), x), y)-(diff(w(x, y, z, t), x))/a+(diff(w(x, y, z, t), y))*(diff(diff(w(x, y, z, t), x), y))))+h^2*(-B11*(diff(diff(diff(w(x, y, z, t), x), x), x))+B12*(-(diff(diff(diff(w(x, y, z, t), x), y), y))-(diff(diff(v(x, y, z, t), x), y))/a))+B66*h^2*(-2*(diff(diff(diff(w(x, y, z, t), x), y), y))-(diff(diff(v(x, y, z, t), x), y))/a)+A66*h*(diff(diff(u(x, y, z, t), y), y)+diff(diff(v(x, y, z, t), x), y)+(diff(diff(w(x, y, z, t), y), y))*(diff(w(x, y, z, t), x))+(diff(w(x, y, z, t), y))*(diff(diff(w(x, y, z, t), x), y))) = rho*(diff(diff(u(x, y, z, t), t), t))-e0^2*a^2*(rho*(diff(diff(diff(diff(u(x, y, z, t), t), t), x), x))+rho*(diff(diff(diff(diff(u(x, y, z, t), t), t), y), y)))

(1)

``


 

Download sort.mw

 

This is the context. I am doing reduction of order on an ODE. This sometimes converts the ode to form  where common factor containing only  x can be moved outside, making the ode looks like  f(x)*(y''(x)+....etc...) = 0 then now I can eliminate f(x)  and only solve (y''(x)+....etc...) = 0 which is simpler.

I found if I can call factor on the ode, it works. It does remove any common terms. The problem I am having is how to cleanly determine the factors obtained. In the above example, it will be f(x) and (y''(x)+....etc...) 

For an example, the ODE   x^2 y'' + x y' = 0 can be written (using factor) as  x ( x y'' + y') =0 and now canceling x<>0, the ode becomes simpler x y''+ y' = 0.

I am now trying to find the two factors, using using op() on the result of factor.

But it does not work for all cases. There should be a way more robust way to obtain the factors.  I give 2 examples to better explain.

Example 1

ode:=x^2*diff(y(x),x$2)-(x-3/16)*y(x)=0:
ode:=expand(algsubs( y(x)=v(x)*x^(1/4)*exp(2*sqrt(x)),ode));
ode:=factor(ode);

After factoring the original ODE, there is common term found, which is the one shown UP. I need to find this to cancel it and keep the rest.

Currently I do this, which does not work for call cases

#check it is factored OK. If the type is `*` then Maple
#found common factor.
if type(lhs(ode),`*`) then 
   #add code to extract the two parts
   LHS:=op([1..-2],lhs(ode));
   RHS:=op(-1,lhs(ode));
fi;

Even though this worked here. I can now check it is the LHS which needs to be canceled. But all what I have are the operands. I do not know how to reconstruct the LHS from the operands. They could be + or *.  If it was `*` between the oprands for LHS, then I can do LHS:=mul(LHS) and this gives 

But I got lucky here. I do not know if it will be `*` all the time for LHS operands.

example 2

ode:=x^2*diff(y(x),x$2)-x*(x+2)*diff(y(x),x)+(x+2)*y(x)=0;
ode:=algsubs( y(x)=v(x)*x,ode):
ode:=factor(ode);

For this, the code I have works, but this is only because the LHS was simple.

if type(lhs(ode),`*`) then #factored OK
   LHS:=op([1..-2],lhs(ode));
   RHS:=op(-1,lhs(ode));
fi;

My question is: I expect factor, if there is common factor, to generate 2 expressions with `*` between them. I am looking for a good way to find what these two factors are. Once I do, it is easy for me to find which is the ODE and which is not and cancel the one which is not out.

I tried collect, but this does not work in general. Since I do not know beforehand, what is the common term, if any, present.

If this needs more clarification please feel free to ask. I have many more examples. This is all done by coding. Non-interactive. So solution based on looking at the output then do something, will not work for me.

 

 

I have been working on some code today and I keep getting errors in my code trying to solve it. I am trying to solve this problem:

y' = (-2xy) + 1 , y(0) = 0

Solve this problem numerically using step size of h=0.1 (101 values) and h = (0.05) (201 values). Create a plot of the two approximate solutions and the analytic solution. Show all work!

I think I got the plot code correct, I must be entering the beginning code incorect as It says:

"Error, (in DEtools/DEplot/CheckDE) extra unknowns found: xy"

Here is my code. Somebody please help me or lead me in the right direction! Thank you!

 

TEST_CODE_1.mw

 

 

 

Dear all

I compute the solution of a system of ordinary differential equaiton using Picard iteration. I have the following error:

Error, invalid input: f uses a 2nd argument, y, which is missing

My code 

Picard_iteration.mw

Thank you for any help 

 

Mac OS 11 Big Sur ist not listed as a compatible system with Maple 2020. Are there known issues? Could someone who already upgraded give feedback?

 

How can I plot the following function?

delta := t -> piecewise(t = 0, 1, 0)
Suspect that the point  at Y=1 automatically scales when zooming in.

 

Thanks

 

 

Currently to parse initial conditions of ODE entered bu user in the form y(x0)=y0 in order to determine x0 and y0 I use patmatch on lhs and rhs separatly which is not a big deal. I wanted to see if I can use an equation inside patmatch.

I could not figure how to use patmatch on both sides of equation on one call.  Help says expression to patmatch is algebraic expression. Does this means an equation can't be used as whole? Here is an example to make it more clear

restart;
ic:=y(0) = 1;
patmatch(ic,y(x0::anything)=y0::anything,'la')

Which gives false. i.e. no match found.

Is there a trick to make the above work if possible? This is just a basic example. I might want to apply this on more complicated equation where more things are on lhs and rhs of =  but patmatch does not like `=` in the middle there. 

 

 

 

Why won't Maple solve any of these inequalities for Q?

At first I tried solving the system of equations, but then I tried solving the inequalities individually for Q, and those too could not be solved by Maple.  What am I doing wrong?
 

restart

 

assume(R1, real, R2, real, RL, real, XC1, real, XC2, real, XL1, real)

additionally(R1 > 0, R2 > 0, RL > 0, XC1 > 0, XC2 > 0, XL1 > 0)

 

 

eq3 := R1*(4*Q^2+1)-RL > 0

0 < R1*(4*Q^2+1)-RL

(1)

eq4 := 4*Q^2*R1*RL-(R1-RL)^2 >= 0

0 <= 4*Q^2*R1*RL-(R1-RL)^2

(2)

eq5 := Rs^2*(RL-R1)/Q^2+R1^2*RL > 0

0 < Rs^2*(RL-R1)/Q^2+RL*R1^2

(3)

 

``

sys1 := {eq3, eq4, eq5}

`assuming`([solve(sys1, {Q})], [R1 > RL])

`assuming`([solve(sys1, {Q})], [R1 < RL])

solve(sys1, {Q})

Warning, solutions may have been lost

 

solve(eq3, Q)

Warning, solutions may have been lost

 

solve(eq4, Q)

Warning, solutions may have been lost

 

solve(eq5, Q)

Warning, solutions may have been lost

 

``


 

Download pi-filter_anal_copy.mw

howca i determine the range of gains if g(x)=1/(x+x^2+x^3+2)

the answer is like closed loop

I need to use the fact that the expression under the radical in a larger expression is >0.

How do I programmaticaly "grab" the radical term to include it in an equation such as for example:

eq1:= b^2-4 a c>0

where lhs term would have come from a previous derivation or solution ?

Thank you...

First 471 472 473 474 475 476 477 Last Page 473 of 2217