Aragorn

20 Reputation

0 Badges

0 years, 176 days

MaplePrimes Activity


These are questions asked by Aragorn

I have to create a function which records the angle and angular frequency of an aspherical moon orbiting a planet at each periapsis to create a Poincare section plot of theta agaisnt omega. I cannot get eval to work out the solutions to the ODEs to plot the points. Additionally I then have to create another function to plot the FFT of the angular frequency from an array. 

This is what I have tried so far:



Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/Project.mw .
 

Download Project.mw

I have to create a simulation of the orbit of Hyperion around Saturn, which should be reliable for up to thousands of orbits. In practice I don't seem to be able to make the model stay consistent above a few hundred orbits using dsolve for the ODEs, but is there any way to do this so that the orbit stays on the same path?

I have uploaded my code here for your perusal.


 

``

with(plots)

G := 0.667408e-10

MH := 0.55855e19

MSat := 0.56832e27

M := MSat+MH

a := 0.1501e10

ecc := .232

beta := .89

Digits := 10

TH := sqrt(4*3.141592653589793^2*a^3/(G*M))

omegaH := (2*3.141592653589793)/TH

Eqns := diff(xH(t), t) = vxH(t), diff(yH(t), t) = vyH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*MSat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5

ICs := xH(0) = a*(1-ecc), yH(0) = 0., vxH(0) = 0., vyH(0) = sqrt(G*M*(1+ecc)/(a*(1-ecc))), theta(0) = 0., omega(0) = omegaH

soln := dsolve({Eqns, ICs}, numeric, method = rkf45, abserr = 0.1e-9, relerr = 0.1e-9, maxfun = 0)

T := 100*TH

odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000)``

odeplot(soln, [theta(t), omega(t)], 0 .. T, labels = ["theta", "omega"], numpoints = 2000)

``


 

Download Project.mw

The dsolve function is not working to solve a set of differential equations I have written. No error message is showing, but the code does not work.
 

``

with*plots

G := 0.667408e-10

0.667408e-10

(1)

M := 0.56832e27+0.55855e19

0.5683200056e27

(2)

a := 0.1501e10

0.1501e10

(3)

ecc := .232

.232

(4)

orbit := proc (T) local Eqns, ICs, soln; Eqns := diff(xH(t), t) = vxH(t), diff(yH(t), t) = vyH(t), diff(vxH(t), t) = G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(vyH(t), t) = G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1-ecc), yH(0) = 0., vxH(0) = 0., vyH(0) = sqrt(G*M*(1+ecc)/(a*(1-ecc))); soln := dsolve({Eqns, ICs}, numeric, method = classical[rk4]); odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000) end proc

proc (T) local Eqns, ICs, soln; Eqns := diff(xH(t), t) = vxH(t), diff(yH(t), t) = vyH(t), diff(vxH(t), t) = G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(vyH(t), t) = G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1-ecc), yH(0) = 0., vxH(0) = 0., vyH(0) = sqrt(G*M*(1+ecc)/(a*(1-ecc))); soln := dsolve({Eqns, ICs}, numeric, method = classical[rk4]); odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000) end proc

(5)

orbit(20)

odeplot(soln, [0.6662225183e-9*xH(t), 0.6662225183e-9*yH(t)], 0 .. 20, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000)

(6)

``


 

Download Project.mw

I am trying to use Maple to provide numerical solutions to the energy levels of muonic lead by using the Shooting method to work out boundaries for the 1s, 2p, and 3d orbitals. However, I am unable to get Maple to graph or fsolve the roots of the function, since it is "unable to evaluate the function to numeric values in this region" despite the zero being within the bounds.

Here is the code:


 

This worksheet covers the energy levels of a muon bound to a lead atom (Pb-206). All plots have been commented out and inserted as pictures. First define the variables for this atom, and the muonic mass ms in units of the free electron mass:

Z := 82; A := 206; ms := 105.66/(.51100)

206.7710372

(1)

where Z is the atomic number and A is the mass number. The Z and A values can be changed for any other atom or isotope. Modelling the atom as hydrogenic, we have the equation

En := Z^2*ms/(2*n^2)

695164.2270/n^2

(2)

where n is the orbital occupied by the muon and En is the energy in Hartree units. The average radius a of the ground state orbit in Bohr is

a := evalf(1/(Z*ms))

0.5897886916e-4

(3)

This is much smaller than the Bohr radius due to the muonic mass being bigger by 2 orders of magnitude and the higher nuclear charge. The atomic radius of the lead atom is determined by

R := evalf(0.125e-14*A^(1/3)/(0.529177e-10))

0.1395076832e-3

(4)

This is over twice the above average orbital radius, meaning that the ground state orbit lies within the nucleus. This means we cannot accurately treat this atom as hydrogenic, since the muon will only be attracted towards the portion of the nuclear charge closer to the centre. Instead, we must use a model of nuclear structure that accounts for the muon being inside the nucleus, treating the nucleus as a spherical charge distribution. Solving the radial Schrodinger equation for the sphere using a uniform spherical charge distribution, then applying Gauss's law for the charge enclosed by the nuclear surface gives a potential distribution which can be modelled using a piecewise function for the different potentials inside and outside the nucleus:

V := proc (r) options operator, arrow; piecewise(r < R, (-1)*.5*(Z/R)*(3-r^2/R^2), -Z/r) end proc

proc (r) options operator, arrow; piecewise(r < R, (-1)*.5*(Z/R)*(3-r^2/R^2), -Z/r) end proc

(5)

with(plots)

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(6)

plot([V(r), -Z/r], r = 0 .. 6*R, V = 1.5*V(0.) .. 0, scaling = unconstrained, axes = boxed, color = [red, black], legend = ["sphere", "point"], title = "Spherical and point potentials", labels = ["r [a.u.]", "V [a.u.]"])

``

While the point charge model unrealistically goes to an infinite potential as the muon gets closer to the centre of the nucleus, the spherical model instead shows it tailing off inside the nucleus as more of the nuclear charge is further from the centre than the muon itself. To model the bounds of the energy of the ground state 1s, we use the Shooting method for boundary value problems.

plotsol := proc (E::float, l::nonnegint, r0::numeric, S::numeric, c::symbol) local Eqns, ICs, fnl, gnl, r, soln; Eqns := diff(fnl(r), r) = gnl(r), diff(gnl(r), r) = l*(l+1)*fnl(r)/r^2-2*ms*(E-V(r))*fnl(r); ICs := fnl(r0) = r0, gnl(r0) = l+1; soln := dsolve({Eqns, ICs}, numeric); plots[odeplot](soln, [r, fnl(r)], r0 .. S, color = c) end proc

NULL

We are looking for a solution with l=0, and the offset r0 to prvent division by 0 must be much smaller than the Bohr radius given the scale of the orbit. We need to estimate bounds for the energy and set an upper bound for the effective range of the potential. This upper bound must be low enough for Maple to calculate the function, which unfortunately means only a few Bohr radii out.

plotsol(-38.5, 0, 0.1e-9, 7.0, red)NULL

 

plotsol(-38.6, 0, 0.1e-9, 7.0, blue)

 

These bounds show that the ground state binding energy is between -39 and -38 Hartrees, since fnl tends to negative infinity at one bound and positive infinity at the other. To home in on the correct value, we plot the function between the bounds and find the zero.

NULL

fnl_at_r := rhs(soln[2])

Error, invalid input: rhs received soln, which is not valid for its 1st argument, expr

 

fnl_at_S := proc (E, r0, l, S) global fnl_at_r; fnl_at_r(parameters = [E, r0, l]); return fnl_at_r(S) end proc

plot(proc (x) options operator, arrow; fnl_at_S(x, 0.1e-9, 0, 4.0) end proc, -38.6 .. -38.5)

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Elevel := fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.10e-19, 0, 7.0) end proc, -10.00 .. 0.)

fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.10e-19, 0, 7.0) end proc, -10.00 .. 0.)

(7)

evalf(Elevel)

fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.10e-19, 0, 7.0) end proc, -10.00 .. 0.)

(8)

The plotted function tends to a straight line with no nodes, so is the solution to 1s. This gives a binding energy of -11.94 Hartrees for the ground 1s state. The 2p state is the lowest with l=1, so the same method using l=1 will find the energy of the 2p state. First find upper and lower bounds:

plotsol(-26.0, 1, 0.1e-9, 8.0, red)

plotsol(-27.0, 1, 0.1e-9, 8.0, blue)

Then solve the differential equation:

plot(proc (x) options operator, arrow; fnl_at_S(x, 0.1e-9, 1, 8.0) end proc, -12.0 .. -10.0)

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

 

 

Elevel := fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.1e-9, 0, 8.0) end proc, -27.0 .. -26.0)

fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.1e-9, 0, 8.0) end proc, -27.0 .. -26.0)

(9)

evalf(Elevel)

fsolve(proc (E) options operator, arrow; fnl_at_S(E, 0.1e-9, 0, 8.0) end proc, -27.0 .. -26.0)

(10)

NULL


 

Download Exercise_3.mw

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

 

Page 1 of 1