MaplePrimes Questions

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

 

Let X be a random such that P[X=x]=x/10, x =1,2,3,4

And let S be the sum 

I want to calculate probabilities like P[S=2] where S is the sum of 2 independent random variable that are distributed like X. The solution should be P[S=2]=P[X=1]*P[X=1]=1/10*1/10=1/100. Or P[S=3]=P[X=1]*P[X=2]+P[X=2]*P[X=1]=4/100

I need this not only for 2 summands but also for 3 and 4 summands. 

How can I code this? I struggle to write a code so that X the random variable has the desired probabilities. Thanks in advance!

 

Thanks in advance

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

First 357 358 359 360 361 362 363 Last Page 359 of 2308