dharr

simplify with side relations...

Answered: dharr 8270

March 11 2025

2 2

>

Warning, The imaginary unit, I, has been renamed _I

>

dS := -P*S*alpha-S*mu+Lambda; dE := alpha*S*P-(-T*eta+1)*beta*E-theta*E-mu*E; dI := (-T*eta+1)*beta*E-delta*I-gamma*I-mu*I; dR := E*theta+I*gamma-R*mu; dP := -P*T*sigma+I*xi-P*tau; dT := r*T*(1-T/K)-phi*T

>

We want the above solution in terms of the following expression.

>

We can use simplify with side relations but we want the above relationship in polynomial form. This can be achieved by multiplying each side by the denominator, which amounts to requiring the following to be zero

>

-(sigma*(r-phi)*K+r*tau)*(delta+gamma+mu)*(R__0-1)*mu/(xi*alpha*beta*(eta*(r-phi)*K-r))

>

Download end.mw

a possibility...

Answered: dharr 8270

March 08 2025

1 3

I don't think there is a simple way to do this, but here is one way to make it convincing. Verifying the actual and your proposed formulas solve the equation is easier and could be a first step.

>

Solve and then replace the variable starting with _Z by k.

>

ans1, ans2 := solve(eqn1, allsolutions); ans1 := eval(ans1, indets(ans1, suffixed(_Z))[] = k); ans2 := eval(ans2, indets(ans2, suffixed(_Z))[] = k)

>

Check each of them individually satisfies the equation (could also apply evalb or is to get true).

>

`assuming`([simplify(eval(eqn1, x = ans1))], [k::integer]); `assuming`([simplify(eval(eqn1, x = ans2))], [k::integer]); `assuming`([simplify(eval(eqn1, x = ans3))], [k::integer])

Of course this doesn't show ans1 and ans2 cover the same values as ans3. We can generate some possibilities to see if this is true, but you have to try different ranges of k to see that they match

>

So this works

>

{-(31/6)*Pi, -(23/6)*Pi, -(19/6)*Pi, -(11/6)*Pi, -(7/6)*Pi, (1/6)*Pi, (5/6)*Pi, (13/6)*Pi, (17/6)*Pi, (25/6)*Pi, (29/6)*Pi} = {-(31/6)*Pi, -(23/6)*Pi, -(19/6)*Pi, -(11/6)*Pi, -(7/6)*Pi, (1/6)*Pi, (5/6)*Pi, (13/6)*Pi, (17/6)*Pi, (25/6)*Pi, (29/6)*Pi}

>

Download sets.mw

convert(...,trigh)...

Answered: dharr 8270

March 07 2025

0 0

The answer to your second question is convert(Bans, trigh). But this can leave complex expressions as arguments, so you might want to also try convert(..., trig).

via products...

Answered: dharr 8270

March 06 2025

4 0

I more-or-less did it by hand here; I'm sure it could be streamlined.

>

restart;

>	s1 := arctan(2/n^2); s2 := expand(convert(s1, ln));

((1/2)*I)*ln(1-(2*I)/n^2)-((1/2)*I)*ln(1+(2*I)/n^2)

Take exp so we can do products instead of sums; find the pieces

>	p:=simplify(exp(s2),exp); p1,e1:=op(op(1,p)); p2,e2:=op(op(2,p));

(1-(2*I)/n^2)^((1/2)*I)*(1+(2*I)/n^2)^(-(1/2)*I)

We can do the products

>	P1:=product(p1, n=1.. infinity); P2:=product(p2, n=1.. infinity);

(-1/2+(1/2)*I)*sin((2+I)*Pi)/Pi

(-1/2-(1/2)*I)*sin((2-I)*Pi)/Pi

>	Take the ln to reverse the exp

>	S:=simplify(e1ln(P1)+e2ln(P2));

>

Download sum.mw

Some tips...

Answered: dharr 8270

March 03 2025

1 0

Here's how I would do it, keeping close to how you did it. When using solve, if you get a warning about assumptions not used, add the useassumptions option.

Levine - Ch. 2.4 - Particle in a Rectangular Well

In this problem, and .

Define

s__1 := sqrt(2*m*(V__0-E))/`&hbar;` = 2^(1/2)*(m*(V__0-E))^(1/2)/`&hbar;`

= -2^(1/2)*(m*(V__0-E))^(1/2)/`&hbar;`

and are wavefunctions.

=

As can be seen above, there are four unknown constants, and .

We can obtain values for these unknowns by imposing boundary conditions.

Some of the boundary conditions involve the first derivatives of the wavefunctions.

=

The boundary conditions are

`ψ__1`(0) = `ψ__2`(0)*`ψ__2`(l) and `ψ__2`(0)*`ψ__2`(l) = `ψ__3`(l)*(D(`ψ__1`))(0) and `ψ__3`(l)*(D(`ψ__1`))(0) = (D(`ψ__2`))(0)*(D(`ψ__2`))(l) and (D(`ψ__2`))(0)*(D(`ψ__2`))(l) = (D(`ψ__3`))(l)

My question is how to find the constants in Maple.

Solving the first boundary condition is easy.

$expr1 := solve(`ψ__1`(0) = `ψ__2`(0), {C}, useassumptions)$ =

Next, I solve the third boundary condition and sub in the result from solving the first boundary condition. We get in terms of .

$expr3 := simplify(eval(solve(`ψ__1,d`(0) = `ψ__2,d`(0), {B}), expr1))$ = ${B = A*(V__0-E)^(1/2)/E^(1/2)}$

Already here it is not clear to me why Maple does not cancel the 's - use simplify.

Next, consider the second boundary condition

$expr2 := solve(`ψ__2`(l) = `ψ__3`(l), {G}, useassumptions)$ = ${G = (A*cos(s*l)+B*sin(s*l))/exp(s__2*l)}$

We can obtain in terms of just by subbing in previously found relationships, as follows

${G = exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)*A*(cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2)+(V__0-E)^(1/2)*sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`))/E^(1/2)}$

Finally, consider the fourth boundary condition.

$solve(`ψ__2,d`(l) = `ψ__3,d`(l), {G}, useassumptions)$

${G = (m*E)^(1/2)*(A*sin(2^(1/2)*(m*E)^(1/2)*l/`&hbar;`)-B*cos(2^(1/2)*(m*E)^(1/2)*l/`&hbar;`))/((m*(V__0-E))^(1/2)*exp(-2^(1/2)*(m*(V__0-E))^(1/2)*l/`&hbar;`))}$

$expr4 := simplify(subs(expr3, solve(`ψ__2,d`(l) = `ψ__3,d`(l), {G}, useassumptions)))$

${G = (-(V__0-E)^(1/2)*cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)+sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2))*A*exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)/(V__0-E)^(1/2)}$

${exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)*A*(cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2)+(V__0-E)^(1/2)*sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`))/E^(1/2) = (-(V__0-E)^(1/2)*cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)+sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2))*A*exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)/(V__0-E)^(1/2)}$

Subtract rhs-lhs, and divide through or multiply by the factors to cancel

result2 := simplify((rhs-lhs)(result[])*sqrt(V__0-E)*sqrt(E)/(exp(sqrt(2)*sqrt(m)*sqrt(V__0-E)*l/`&hbar;`)*A))

For example, I would like to replace with . But s already has a value, so you need to use a different symbol or unassign it

Since this is essentially replacing all the ms

s = sqrt(2*m*E)/`&hbar;`; isolate(%, m); simplify(eval(result2, %))

(2*E-V__0)*sin((s^2/E)^(1/2)*E^(1/2)*l)-2*cos((s^2/E)^(1/2)*E^(1/2)*l)*E^(1/2)*(V__0-E)^(1/2)

This is a first task (that I have asked about in a separate question).

Then, I would want Maple to do the cancellations for me, but I have already used "simplify".

Finally, I think I would like to collect terms involving and .

Here is what happens when I try to get Maple to solve the equations for me in one go

$solve({`ψ__1`(0) = `ψ__2`(0), `ψ__1,d`(0) = `ψ__2,d`(0), `ψ__2`(l) = `ψ__3`(l), `ψ__2,d`(l) = `ψ__3,d`(l)}, {A, B, C, G}, useassumptions)$

Download solve_constants.mw

Error fixed; LambertW...

Answered: dharr 8270

February 28 2025

2 1

You did not change u(x,y,t) to u(x,y,z,t) in the code that extracts the linear part. After fixing that, there is no solution.
(The answer to the question about removing LambertW is that you can't).

remove2.mw

using indets...

Answered: dharr 8270

February 26 2025

2 2

Basically, using indets gives you the (more complicated) subexpressions that you can directly substitute. I did it twice so there are no exps or square roots left, but you can go further if you wish.

Download expFT_compact.mw

P.S. If the matrix were diagonalizable, I could probably do something more intelligent via the eigenvectors. Is the last row really intended to be zero?

convert(..,radical)...

Answered: dharr 8270

February 25 2025

2 0

The normal way to deal with a RootOf if you want an explicit solution would be to convert it to radical form. In your case, it is the root(s) of a quadratic, so you can use convert(E[u], radical) to get an explicit solution in terms of the quadratic formula. If you want both roots, you can use allvalues(E[u]). Or better yet, you probably got that result from solve, so you could add the option explicit to the solve command and not generate the RootOf in the first place. (remove_RootOf is just giving you an equivalent equation to solve, so you are not getting ahead by that method).

In future please upload your worksheet using the green up-arrow in the Mapleprimes editor, select your file, click upload, and then insert link or insert contents.

One lump answer...

Answered: dharr 8270

February 21 2025

1 51

Here the derivation for the one lump case (sections in green). The two-lump case can be done similarly, I assume.

>

pde := diff(u(x, y, t), t)-(diff(diff(u(x, y, t), `$`(x, 4))+5*u(x, y, t)*(diff(u(x, y, t), `$`(x, 2)))+(5/3)*u(x, y, t)^3+5*(diff(u(x, y, t), x, y)), x))-5*u(x, y, t)*(diff(u(x, y, t), y))+5*(int(diff(u(x, y, t), `$`(y, 2)), x))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

>

diff(u(x, y, t), t)-(diff(diff(diff(diff(diff(u(x, y, t), x), x), x), x), x))-5*(diff(diff(diff(u(x, y, t), x), x), y))+5*(int(diff(diff(u(x, y, t), y), y), x)), -5*(diff(u(x, y, t), x))*(diff(diff(u(x, y, t), x), x))-5*u(x, y, t)*(diff(diff(diff(u(x, y, t), x), x), x))-5*u(x, y, t)^2*(diff(u(x, y, t), x))-5*u(x, y, t)*(diff(u(x, y, t), y))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

>

Bij := proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

>

f(x, y, t) = (-5*t*p[1]^2+y*p[1]+x)*(-5*t*p[2]^2+y*p[2]+x)+(-6*p[1]-6*p[2])/(p[1]-p[2])^2

>

Eqs 16,17

>	$eqp := {x = xp-(5(a^2+b^2))t, y = 10at+yp}; eqfc2 := f(x, y, t) = simplify(eval(rhs(eqfcomplex), eqp))$

f(x, y, t) = (b^4*yp^2+(a*yp+xp)^2*b^2+3*a)/b^2

>

eq17 := u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2

u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2

Eq 18

>

u(x, y, t) = -12*(-b^4*yp^2+(a*yp+xp)^2*b^2-3*a)*b^2/(b^4*yp^2+(a*yp+xp)^2*b^2+3*a)^2

Alternatively,

>

eqfp := dchange(eqp, eqfcomplex, [xp, yp], params = [a, b], simplify); eq17p := dchange(eqp, eq17, [xp, yp], params = [a, b], simplify); eqt := simplify(eval(eq17p, eqfp))

f(xp, yp, t) = (b^4*yp^2+(a*yp+xp)^2*b^2+3*a)/b^2

u(xp, yp, t) = (6*(diff(diff(f(xp, yp, t), xp), xp))*f(xp, yp, t)-6*(diff(f(xp, yp, t), xp))^2)/f(xp, yp, t)^2

u(xp, yp, t) = -12*(-b^4*yp^2+(a*yp+xp)^2*b^2-3*a)*b^2/(b^4*yp^2+(a*yp+xp)^2*b^2+3*a)^2

Eq 17 is the parametric equation for a line.
velocities are derivatives wrt t. (after eq 19)

Intercept is at the origin (x=0,y=0 at t=0), and slope is dy/dx (Eq 19).

>

vx, vy := diff(eval(x, eqp), t), diff(eval(y, eqp), t); dydx := simplify(vy/vx)

>

Download Y-pde-line.mw

Maple can do this one...

Answered: dharr 8270

February 21 2025

3 0

Download limit.mw

answer to the followup question...

Answered: dharr 8270

February 21 2025

2 1

Here is an answer to a followup question, now deleted, in which the evalc and pdetest were much slower because of a more complicated function f. I'm guessing the integrals in the original pde are problematic, so suggest working with the pde in f. For the evalc, if you do it on the simpler f and not u it is also manageable.

Question file: stop-working_question.mw

>

pde := diff(u(x, y, t), t)-(diff(diff(u(x, y, t), `$`(x, 4))+5*u(x, y, t)*(diff(u(x, y, t), `$`(x, 2)))+(5/3)*u(x, y, t)^3+5*(diff(u(x, y, t), x, y)), x))-5*u(x, y, t)*(diff(u(x, y, t), y))+5*(int(diff(u(x, y, t), `$`(y, 2)), x))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

>

eq17 := u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2

The integrals in pde are problematic for fast evaluation, so work in the pde for f

>

Bij := proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

>

memory used=0.63GiB, alloc change=483.86MiB, cpu time=98.80s, real time=18.07s, gc time=296.88ms

evalc is slow on the complicated expression, so do it on the expression for f first

>

And now substitute to get the expression for u

>

Download stop-working.mw

identical...

Answered: dharr 8270

February 19 2025

1 5

Download identical.mw

_options and _rest...

Answered: dharr 8270

February 19 2025

0 0

Maple has the _options mechanism for this purpose, see the help page. Passing _options to your embedded procedure passes all the keyword arguments. If you just want to pass some of them, such as color, then use the syntax _options['color'].

If you don't need to override the default values used by plot, then the simplest is to use _rest. In this case, use

ecdf := proc(S) ... end proc;

and pass _rest to the plot routines; the rest of the arguments (aside from S) that were passed are passed on to the plot routines. This of course doesn't do any argument checking, so has some downsides. Or you can use some combined strategy.

>

restart

>

ecdf := proc(S, {color::string:="blue", thickness::posint:=1})
  local N   := numelems(S):
  local M   := sort(S):
  local i;
  plots:-display(
    plot([ seq(op([[M[i], (i-1)/N], [M[i], i/N], [M[i+1], i/N]]), i=1..N-1) ], _options)
    , plot([ [M[N], (N-1)/N], [M[N], 1] ], _options)
    , plot([ [M[N], 1], [(max+(max-min)*0.2)(S), 1] ], _options, 'linestyle'=3)
    , plot([ [M[1], 0], [(min-(max-min)*0.2)(S), 0] ], _options, 'linestyle'=3)
  )
end proc:

>	S := Statistics:-Sample(Uniform(0, 1), 10): ecdf(S)

>	ecdf(S, color="red", thickness=3)

>

Download ecdf.mw

ilog10...

Answered: dharr 8270

February 16 2025

0 1

If I understand correctly, you want 99.9999999 to return 1, not 2. Then ilog10 will work.

Some comments....

Answered: dharr 8270

February 16 2025

0 0

In general, if there is a symbolic solution, it is preferred, since the sorts of numerical issues you get here don't occur. Here my workarounds, and comments in italics.

Note that ConditionNumber(Snum) returns a very high value, which is probably why LinearSolve (and fsolve) have difficulty finding both solutions; this will generally be the case when you are solving (A-lambda*I)x =0, (where lambda is an eigenvalue).

Below is an eigenvector equation for for the eigenvalue -3.4376...

S := AA-evs[1]*IdentityMatrix(2); Rank(S); Snum := AA-evsnum[1]*IdentityMatrix(2); Rank(Snum)

When I solve the system below, I expect to get the vector

, just like the result of the Eigenvectors command used above.

Instead I get the zero vector. You can use NullSpace to get around this.

If you do it symbolically, you find there is a degree of freedom with LinearSolve. Different values of the variable starting with an underline can give both the zero vector, and a multiple of the eigenvector (whose length is arbitrary). Since the zero vector is not in the nullspace, you only get the eigenvector with NullSpace.

In addition, when I compute the reduced row echelon form, I expect to have one of the rows be zero since matrix S has rank 1.

We get the expected results with a symbolic solution

Therefore, I guess my question is maybe about calculations and not about Maple. That being said, I did the calculations manually, then I asked chatgpt, which agreed with me.

Why are the calculations not coinciding with Maple?

I can do the calculations "manually" with Maple, For numeric calculations use fsolve