Advanced Search

dharr

Dr. David Harrington

9152 Reputation

22 Badges

21 years, 259 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

near enough...

dharr 9152
7 hours ago
0 2

@salim-barzani The paper seems to have several errors

fp2.mw

(1) volumes and areas...

dharr 9152
May 08 2026
3 6

Here for (1): [edit: updated to use the InRadius.]

> 

restart
> 

with(geom3d); _local(O)

Define o as the origin, and C and O as cube and octahedron. The cube has sides of length 1. The corners of the octahedron will touch the centres of the faces of the cube, so the circumsphere of the octahedron will be equal to 1/2, which can be found as the radius of the in-sphere of the cube.

> 

point(o, 0, 0, 0); cube(C, o, side = 1); inR := InRadius(C); octahedron(O, o, inR)

1/2

> 

draw([C(cutout = .85), O])

And the ratio of the volumes is

> 

vOvC := volume(O)/volume(C)

1/6

For the cube inside the octahedron, its corners touch the centres of faces of the octahedron. So we want the distance of the centre of a face from the origin, which is the InRadius

> 

d := InRadius(O)

(1/6)*3^(1/2)

So the cube inside has circumsphere of radius d

> 

cube(C__2, o, d)
> 

draw([C(cutout = .85), O(cutout = .85), C__2])

The ratio of volumes of successive cubes is

> 

rc := volume(C__2)/volume(C)

1/27

And so the sum of all the volumes of the cubes is a geometric series "1+1/(27)+(1/(27))^(2)+..."

> 

cubesvol := 1/(1-rc)

27/26

The volume of the next octahedron is 1/6 of the volume of C__2, from which we get the ratio of volumes of successive octahedra is also 1/27

> 

vO2 := vOvC*volume(C__2); ro := vO2/volume(O)

1/162

1/27

Summing the geoemetric series, starting at the volume of the first octahedron gives

> 

octsvol := volume(O)/(1-ro)

9/52

And the grand total volume is

> 

allvol := cubesvol+octsvol

63/52

For the surface areas we have for the cubes

> 

ac := area(C__2)/area(C); cubesarea := area(C)/(1-ac)

1/9

27/4

and for the octahedra it is again the same ratio

> 

aO2 := area(O)*area(C__2)/area(C); ao := aO2/area(O)

(1/9)*3^(1/2)

1/9

> 

octsarea := area(O)/(1-ao)

(9/8)*3^(1/2)

And the grand total area is

> 

cubesarea+octsarea

27/4+(9/8)*3^(1/2)

> 

``

Download cubeoct2.mw

Conflict with Statistics...

dharr 9152
May 05 2026
1 1

The error message mentions a problem with Statistics:-Remove, so the problem is that the dataframe Remove and the Statistics Remove commands are being confused. If you remove your with(Statistics) line, you will see that many of your Removes now work. To use it while the Statistics package is loaded you can use

df1:-Remove(df1, A);

or

DataFrame:-Remove(df1, A);

This syntax is because the dataframes are implemented as Maple Objects, but this use of Remove doesn't seem to be well documented. 

ZehfussReplace...

dharr 9152
April 21 2026
2 5

I think this is what you want. Maybe you want to use uneval quotes around Determinant in the last line of the transformer procedure (or use %Determinant)

> 

restart
> 

with(LinearAlgebra)
> 

A := Matrix(2, 2, symbol = a); B := Matrix(3, 3, symbol = b)
> 

ZehfussReplace := proc(expr::uneval)
  uses LinearAlgebra;
  subsindets(expr,
    specfunc(specfunc('KroneckerProduct'),'Determinant'),
    proc(q)
      local AB;
      AB:=op(op(q));
      if not eval([AB])::['Matrix'('square'),'Matrix'('square')] then return eval(q) end if;
      Determinant(AB[1])^(upperbound(AB[2])[1])*Determinant(AB[2])^(upperbound(AB[1])[1])
    end proc
   );
end proc:
> 

ZehfussReplace(Determinant(KroneckerProduct(A, B))+2*Determinant(KroneckerProduct(A, A)))

(a[1, 1]*a[2, 2]-a[1, 2]*a[2, 1])^3*(b[1, 1]*b[2, 2]*b[3, 3]-b[1, 1]*b[2, 3]*b[3, 2]-b[1, 2]*b[2, 1]*b[3, 3]+b[1, 2]*b[2, 3]*b[3, 1]+b[1, 3]*b[2, 1]*b[3, 2]-b[1, 3]*b[2, 2]*b[3, 1])^2+2*(a[1, 1]*a[2, 2]-a[1, 2]*a[2, 1])^4

> 

NULL

Download ZehfussReplace.mw

Convert to...

dharr 9152
April 18 2026
1 4

After inserting the table, select it all with the mouse (all cells now yellow). Then the right-click menu has "Convert To" with option "2-D math input", which will put prompts in each cell ready for 2-D math input. Is this what you want?

Some things...

dharr 9152
April 16 2026
1 18

Here's a corrected version:

second_try.mw

I made the following changes:
1. restart in separate execution group as stated in the help page.

2. I changed d2lambda to dlambda2 (your main error).

3. Digits := 10. Although Digits = 15 means hardware precision, some algorithms require a few guard digits. At Digits=10 (default) you get hardware precision anyway but some reserve in error tolerance. (If you try Digits = 13, you will get an error about not achieving the required error tolerance - if you really need another digit or so you could play with specifying error tolerances for the integrals, but it's probably not worth it.)

4. I used evalc() instead of Re() in your integrands. This reduces the time from about 40 seconds to 3 seconds. Your final result is real (as I think you expected). Doing all those complex calculations is tough. Each integral apparently has imaginary parts, but if you know they simplify to zero, you can do evalc(Re()) instead, which gives the same result slightly faster.

cost = 22666 - no it's 15468...

dharr 9152
April 14 2026
3 9

After Rouben's comment, I've corrected my answer. See meadow4.mw below for how to get the solution in terms of elliptic functions directly from dsolve.

https://www.mapleprimes.com/questions/242539-How-Much-Will-The-Meadow-Cost

> 

restart
> 

with(VariationalCalculus); with(DEtools); with(IntegrationTools); with(PDEtools); with(plots)

Work with the upper part only. So a rectangle at x with base dx and height y has area dA = y*dx and costs dc = k*x*y*dx with k = 1(currency unit) per cubic metre. Suppose the fence crosses the x axis at x = x__max.

We need to maximize the cost

> 

c := Int(k*x*y(x), x = 0 .. x__max)

Int(k*x*y(x), x = 0 .. x__max)

subject to the constraint that the (upper part) arc length is 50 m. That arc length is

> 

s := Int(sqrt(1+(diff(y(x), x))^2), x = 0 .. x__max)

Int((1+(diff(y(x), x))^2)^(1/2), x = 0 .. x__max)

Using a Lagrange multiplier `λ__1` with units (currency unit) per metre we want to maximize

> 

s*`λ__1`+c

(Int((1+(diff(y(x), x))^2)^(1/2), x = 0 .. x__max))*lambda__1+Int(k*x*y(x), x = 0 .. x__max)

We can define a new lambda = `λ__1`/k (units: square metre) and maximize

> 

expand(c/k)+lambda*s; L := combine(%)

Int(x*y(x), x = 0 .. x__max)+lambda*(Int((1+(diff(y(x), x))^2)^(1/2), x = 0 .. x__max))

Int(lambda*(1+(diff(y(x), x))^2)^(1/2)+x*y(x), x = 0 .. x__max)

EulerLagrange takes the integrand and gives the differential equation to solve (=0 is implied)

> 

ode := EulerLagrange(GetIntegrand(L), x, y(x))[]

x+lambda*(diff(y(x), x))^2*(diff(diff(y(x), x), x))/(1+(diff(y(x), x))^2)^(3/2)-lambda*(diff(diff(y(x), x), x))/(1+(diff(y(x), x))^2)^(1/2)

Let's change variables to X = x/x__max, Y(X) = y(x)/x__max, mu = lambda/x__max^2 (all dimensionless) to find a new ode, a dimensionless arc length S and a dimensionless cost C.

> 

itr := {X = x/x__max, mu = lambda/x__max^2, Y(X) = y(x)/x__max}; tr := solve(itr, {lambda, x, y(x)}); ode2 := dchange(tr, ode/x__max, [X, Y(X), mu], params = [x__max], simplify); S := dchange(tr, s/x__max, [X, Y(X), mu], params = [x__max], simplify); C := dchange(tr, c/(k*x__max^3), [X, Y(X), mu], params = [x__max], simplify)

{X = x/x__max, mu = lambda/x__max^2, Y(X) = y(x)/x__max}

{lambda = mu*x__max^2, x = X*x__max, y(x) = Y(X)*x__max}

((1+(diff(Y(X), X))^2)^(3/2)*X-(diff(diff(Y(X), X), X))*mu)/(1+(diff(Y(X), X))^2)^(3/2)

Int((1+(diff(Y(X), X))^2)^(1/2), X = 0 .. 1)

Int(X*Y(X), X = 0 .. 1)

Solve numerically with boundary conditions Y(0) = 0, Y(1) = 0 for some mu values.

> 

muvals := [-10, -5, -2, -.5]; colors := [red, blue, black, magenta]; display(seq(odeplot(dsolve({eval(ode2, mu = muvals[i]), Y(0) = 0, Y(1) = 0}, Y(X), numeric), color = colors[i]), i = 1 .. nops(muvals)))

[-10, -5, -2, -.5]

[red, blue, black, magenta]

We learn from this that mu is negative.
However, the slope needs to be -infinity at X = 1, otherwise the above-axis and below-axis parts of the solution will have a discontinuity in slopes.

So we want the boundary condition at X = 1 to be that the slope is -∞ there.

> 

sol2 := `assuming`([factor(sort([dsolve({ode2, Y(0) = 0, (D(Y))(1) = -infinity}, Y(X))])[1])], [mu < 0])

Y(X) = int((_z1^2+2*mu-1)/(-(_z1-1)*(_z1+1)*(_z1^2+4*mu-1))^(1/2), _z1 = 0 .. X)

I'm a bit shocked that thing with the infinity worked. Here's an alternative method. We find the first integral

> 

ode3 := firint(ode2)

(diff(Y(X), X))/(1+(diff(Y(X), X))^2)^(1/2)-(1/2)*X^2/mu+c__1 = 0

As proc (X) options operator, arrow; 1 end proc we want "diff(Y(X),X)->-infinity". Use this condition to find c__1.

> 

eval(limit(eval(ode3, diff(Y(X), X) = dydx), dydx = -infinity), X = 1); _C1 = solve(%, _C1); ode4 := eval(ode3, %)

(1/2)*(2*c__1*mu-2*mu-1)/mu = 0

c__1 = (1/2)*(2*mu+1)/mu

(diff(Y(X), X))/(1+(diff(Y(X), X))^2)^(1/2)-(1/2)*X^2/mu+(1/2)*(2*mu+1)/mu = 0

Solve ode4 for diff(Y(X), X) and then integrate to get Y(X). The solution is similar to the solution obtained above, but the sign of mu is changed.

> 

ans := sort([solve(ode4, diff(Y(X), X))])[2]; sol4 := Y(X) = factor(Int(subs(X = x, ans), x = 0 .. X))

-(X^2-2*mu-1)/(-X^4+4*X^2*mu+2*X^2-4*mu-1)^(1/2)

Y(X) = Int((-x^2+2*mu+1)/((x-1)*(x+1)*(-x^2+4*mu+1))^(1/2), x = 0 .. X)

Which is correct? Try plotting them. sol4 looks good.

> 

plot(eval(rhs(sol4), mu = -.3), X = 0 .. 1)

But no plot for sol2.

> 

sol2b := convert(eval(rhs(sol2), mu = -.3), Int); plot(sol2b, X = 0 .. 1)

Warning, unable to determine if -1 is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if 1 is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if (-4*mu+1)^(1/2) is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if -(-4*mu+1)^(1/2) is between 0 and X; try to use assumptions or use the AllSolutions option

Int((_z1^2-1.6)/(-(_z1-1)*(_z1+1)*(_z1^2-2.2))^(1/2), _z1 = 0 .. X)

Verify that the solution is complex

> 

evalf(eval(rhs(sol2), {X = .5, mu = -.3}))

Warning, unable to determine if -1 is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if 1 is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if (-4*mu+1)^(1/2) is between 0 and X; try to use assumptions or use the AllSolutions option

Warning, unable to determine if -(-4*mu+1)^(1/2) is between 0 and X; try to use assumptions or use the AllSolutions option

.5448754549*I

So we continue with sol4. (sol2 works with positive mu, so it looks as though the assumption about mu did not work?)

Is the solution real as proc (mu) options operator, arrow; 0 end proc? We need the part under the square root to be positive

> 

z := denom(ans)^2

-X^4+4*X^2*mu+2*X^2-4*mu-1

We have only even powers so we can consider it as a quadratic in X^2 = W. Since the coefficient of W^2 is negative, it is an inverted parabola, which is positive between its roots.

> 

algsubs(X^2 = W, z); solve(%, W)

-W^2+4*W*mu+2*W-4*mu-1

1, 4*mu+1

So if we have mu <= -1/4, then it is positive over the whole range 0 <= W and W <= 1, which is also 0 <= X and X <= 1. But for mu > -1/4 the integrand will be somewhere complex and we do not have a real solution.

> 

sol5 := `assuming`([Y(X) = value(rhs(sol4))], [mu < -1/4])

Y(X) = -(-X^2+1)^(1/2)*(X^2-4*mu-1)^(1/2)*(2*mu*EllipticF(X, I/(-4*mu-1)^(1/2))-4*mu*EllipticE(X, I/(-4*mu-1)^(1/2))-EllipticE(X, I/(-4*mu-1)^(1/2)))/((-4*mu-1)^(1/2)*(-X^4+4*X^2*mu+2*X^2-4*mu-1)^(1/2))

Still good

> 

plot(eval(rhs(sol5), mu = -.3), X = 0 .. 1)

What is the limit as proc (X) options operator, arrow; 1 end proc?

> 

Y1 := limit(rhs(sol5), X = 1, left)

-(2*mu*EllipticK(I/(-4*mu-1)^(1/2))-4*mu*EllipticE(I/(-4*mu-1)^(1/2))-EllipticE(I/(-4*mu-1)^(1/2)))/(-4*mu-1)^(1/2)

So there is only one mu value that works

> 

plot(Y1, mu = -5 .. .25)

No symbolic solution

> 

solve(Y1, mu)

-(1/4)*(RootOf(_Z^2*EllipticK(I*_Z)-2*EllipticE(I*_Z)+EllipticK(I*_Z), -2.179660432-0.*I)^2+1)/RootOf(_Z^2*EllipticK(I*_Z)-2*EllipticE(I*_Z)+EllipticK(I*_Z), -2.179660432-0.*I)^2, -(1/4)*(RootOf(_Z^2*EllipticK(I*_Z)-2*EllipticE(I*_Z)+EllipticK(I*_Z), 2.179660432-0.*I)^2+1)/RootOf(_Z^2*EllipticK(I*_Z)-2*EllipticE(I*_Z)+EllipticK(I*_Z), 2.179660432-0.*I)^2

> 

eqmu := mu = fsolve(Y1, mu = -.5 .. -.25)

mu = -.3026213915

> 

plot(eval(rhs(sol5), eqmu), X = 0 .. 1)

The dimensionless arc length is

> 

factor(eval(S, sol4)); S2 := `assuming`([simplify(value(%))], [mu < -1/4])

Int((-4*mu^2/((X-1)*(X+1)*(X^2-4*mu-1)))^(1/2), X = 0 .. 1)

-2*mu*EllipticK(1/(4*mu+1)^(1/2))/(-4*mu-1)^(1/2)

which means the real arc length, which we need to be 50 m (for the top half only), must satisfy

> 

eqfence := 50 = x__max*S2

50 = -2*x__max*mu*EllipticK(1/(4*mu+1)^(1/2))/(-4*mu-1)^(1/2)

The dimensionless cost is

> 

C2 := `assuming`([simplify(value(eval(C, sol5)))], [mu < -1/4])

-(4/3)*((mu^2-(3/4)*mu-1/4)*EllipticE(1/(4*mu+1)^(1/2))-mu*(mu-1/2)*EllipticK(1/(4*mu+1)^(1/2)))/(-4*mu-1)^(1/2)

which means the real cost (taking k = 1) and including the factor of 2 to account also for the bottom half is

> 

cost := 2*x__max^3*C2

-(8/3)*x__max^3*((mu^2-(3/4)*mu-1/4)*EllipticE(1/(4*mu+1)^(1/2))-mu*(mu-1/2)*EllipticK(1/(4*mu+1)^(1/2)))/(-4*mu-1)^(1/2)

We can solve the fence equation for x__max and use that to find the cost in terms of mu.

> 

cost2 := eval(cost, x__max = solve(eqfence, x__max))

(125000/3)*(-4*mu-1)*((mu^2-(3/4)*mu-1/4)*EllipticE(1/(4*mu+1)^(1/2))-mu*(mu-1/2)*EllipticK(1/(4*mu+1)^(1/2)))/(mu^3*EllipticK(1/(4*mu+1)^(1/2))^3)

And for the numerical value of mu the final cost is

> 

finalcost := simplify(eval(cost2, eqmu), zero)

15468.55912

Now we can find x__max

> 

xmaxnum := simplify(solve(eval(eqfence, eqmu), x__max), zero)

39.15937451

Check fence length (top half)

> 

eval(eqfence, {eqmu, x__max = xmaxnum})

50 = 49.99999996-0.*I

In real units the fence curve is

> 

ynum := eval(x__max*(eval(rhs(sol5), {eqmu, X = x/x__max})), x__max = xmaxnum)

-85.35413912*(-0.6521214830e-3*x^2+1)^(1/2)*(0.6521214830e-3*x^2+.210485566)^(1/2)*(-.6052427830*EllipticF(0.2553666938e-1*x, 2.179660431*I)+.210485566*EllipticE(0.2553666938e-1*x, 2.179660431*I))/(-0.4252624286e-6*x^4+0.5148593235e-3*x^2+.210485566)^(1/2)

> 

plot([ynum, -ynum], x = 0 .. xmaxnum, numpoints = 500, labels = [x, y], scaling = constrained)

> 

NULL

meadow3.mw

Here's an update that uses more assumptions to get dsolve to directly give a solution in terms of elliptic functions, and simplifies the expressions somewhat:

meadow4.mw

Old plot:

> 

NULL

Old worksheet:

Download meadow2.mw

uneval...

dharr 9152
April 11 2026
1 1

Use of showstat(ExpressionTools:-Compare) gives

ExpressionTools:-Compare := proc(XX::uneval, VV::uneval,...

so your normal expectations about full evaluation do not apply. Your last example gives the workaround.

Mostly...

dharr 9152
March 28 2026
2 10

Please use the notation in the paper to make it easier to follow. A doi link would also be nice.
Edit: Now correctly runs through and is validated by pdetest.

> 

restart
> 

with(PDEtools)
> 

undeclare(prime, quiet)
> 

declare(u(x, y, t), quiet); declare(v(x, y, t), quiet)
> 

pde1 := diff(u(x, y, t), t, y)-(diff(u(x, y, t), `$`(x, 2), y))+2*(diff(u(x, y, t)*(diff(u(x, y, t), x)), y))+2*(diff(x, y, t, `$`(x, 2)))

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

> 

pde2 := diff(v(x, y, t), t)+diff(v(x, y, t), `$`(x, 2))+2*(diff(u(x, y, t)*v(x, y, t), x))

diff(v(x, y, t), t)+diff(diff(v(x, y, t), x), x)+2*(diff(u(x, y, t), x))*v(x, y, t)+2*u(x, y, t)*(diff(v(x, y, t), x))

> 

Hexpr := h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4

h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4

> 

eq3 := diff(phi(xi), xi) = sqrt(Hexpr)

diff(phi(xi), xi) = (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2)

We have balance in pde1 is n=1 and in pde2 n=2 giving Eq 24

> 

Seriesu := a[0](y, t)+sum(a[j](y, t)*phi(xi)^j+b[j](y, t)*phi(xi)^(j-1)*sqrt(Hexpr), j = 1 .. 1)

a[0](y, t)+a[1](y, t)*phi(xi)+b[1](y, t)*(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2)

> 

Seriesv := A[0](y, t)+sum(A[j](y, t)*phi(xi)^j+B[j](y, t)*phi(xi)^(j-1)*sqrt(Hexpr), j = 1 .. 2)

A[0](y, t)+A[1](y, t)*phi(xi)+B[1](y, t)*(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2)+A[2](y, t)*phi(xi)^2+B[2](y, t)*phi(xi)*(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2)

And we want

> 

eqxi := xi = P(y, t)*x+Q(y, t)

xi = P(y, t)*x+Q(y, t)

> 

substs := eval({u(x, y, t) = Seriesu, v(x, y, t) = Seriesv}, eqxi)

{u(x, y, t) = a[0](y, t)+a[1](y, t)*phi(P(y, t)*x+Q(y, t))+b[1](y, t)*(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(1/2), v(x, y, t) = A[0](y, t)+A[1](y, t)*phi(P(y, t)*x+Q(y, t))+B[1](y, t)*(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(1/2)+A[2](y, t)*phi(P(y, t)*x+Q(y, t))^2+B[2](y, t)*phi(P(y, t)*x+Q(y, t))*(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(1/2)}

Substitute into the pdes - notice we will need up to the third derivatives of phi

> 

eqsa := eval([pde1, pde2], substs); indets(%)

{t, x, y, h[0], h[1], h[2], h[3], h[4], 1/(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(5/2), 1/(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(3/2), 1/(h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(1/2), (h[0]+h[1]*phi(P(y, t)*x+Q(y, t))+h[2]*phi(P(y, t)*x+Q(y, t))^2+h[3]*phi(P(y, t)*x+Q(y, t))^3+h[4]*phi(P(y, t)*x+Q(y, t))^4)^(1/2), P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(P(y, t)*x+Q(y, t)), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t), (D(phi))(P(y, t)*x+Q(y, t)), ((D@@2)(phi))(P(y, t)*x+Q(y, t)), ((D@@3)(phi))(P(y, t)*x+Q(y, t))}

We can go back to xi now. Do it the hard way to maintain the third derivative

> 

eqsb := subs(((D@@3)(phi))(P(y, t)*x+Q(y, t)) = diff(phi(xi), `$`(xi, 3)), ((D@@2)(phi))(P(y, t)*x+Q(y, t)) = diff(phi(xi), `$`(xi, 2)), (D(phi))(P(y, t)*x+Q(y, t)) = diff(phi(xi), xi), phi(P(y, t)*x+Q(y, t)) = phi(xi), eqsa); indets(%)

{t, x, xi, y, h[0], h[1], h[2], h[3], h[4], 1/(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(5/2), 1/(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(3/2), 1/(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2), P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(diff(phi(xi), xi), xi), xi), diff(diff(phi(xi), xi), xi), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(phi(xi), xi), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(xi), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

Use eq 3 and dsubs to deal with all the derivatives of phi

> 

eqsc := dsubs(eq3, eqsb); indets(%)

{t, x, xi, y, h[0], h[1], h[2], h[3], h[4], 1/(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(3/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(5/2), P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(xi), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

Normalize and take the numerators

> 

eqsd := map(`@`(numer, normal), eqsc); indets(%)

{t, x, xi, y, h[0], h[1], h[2], h[3], h[4], (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(3/2), (h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(5/2), P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(xi), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

> 

eqse := simplify(eqsd, {Hexpr = H}); indets(%)

{H, t, x, xi, y, h[1], h[2], h[3], h[4], H^(1/2), H^(3/2), H^(5/2), P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(xi), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

Need to get rid of the half-integer powers. rH = sqrt(H)

> 

eqsf := `assuming`([simplify(eval(eqse, H = rH^2))], [rH > 0]); indets(%)

{rH, t, x, xi, y, h[1], h[2], h[3], h[4], P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), phi(xi), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

Finally, we can collect powers of rH and phi and find their coefficients.

> 

cfs1 := coeffs(collect(eqsf[1], {rH, phi(xi)}, distributed), {rH, phi(xi)}, 'monomials1'); monomials1; cfs2 := coeffs(collect(eqsf[2], {rH, phi(xi)}, distributed), {rH, phi(xi)}, 'monomials2'); monomials2

rH^4, rH^6, rH^5, rH^6*phi(xi)^2, rH^6*phi(xi), rH^5*phi(xi)^3, rH^5*phi(xi)^2, rH^5*phi(xi), rH^4*phi(xi)^6, rH^4*phi(xi)^5, rH^4*phi(xi)^4, rH^4*phi(xi)^3, rH^4*phi(xi)^2, rH^4*phi(xi)

rH^2, rH, phi(xi)*rH^2, rH^3*phi(xi), rH*phi(xi)^5, rH*phi(xi)^4, rH*phi(xi)^3, rH*phi(xi)^2, rH*phi(xi), rH^2*phi(xi)^2, rH^2*phi(xi)^4, rH^2*phi(xi)^3, rH^3, rH^4

The coefficients are an overdetermined set of pdes

> 

pdes := {cfs1, cfs2}; indets(%)

{t, x, y, h[1], h[2], h[3], h[4], P(y, t), Q(y, t), diff(P(y, t), t), diff(P(y, t), y), diff(Q(y, t), t), diff(Q(y, t), y), diff(diff(P(y, t), t), y), diff(diff(Q(y, t), t), y), diff(diff(a[0](y, t), t), y), diff(diff(a[1](y, t), t), y), diff(diff(b[1](y, t), t), y), diff(A[0](y, t), t), diff(A[1](y, t), t), diff(A[2](y, t), t), diff(B[1](y, t), t), diff(B[2](y, t), t), diff(a[0](y, t), t), diff(a[0](y, t), y), diff(a[1](y, t), t), diff(a[1](y, t), y), diff(b[1](y, t), t), diff(b[1](y, t), y), A[0](y, t), A[1](y, t), A[2](y, t), B[1](y, t), B[2](y, t), a[0](y, t), a[1](y, t), b[1](y, t)}

Solve them

> 

sol := [pdsolve(pdes)]; nops(%)

15

There are different solutions, varying in signs, with the following as one of the desired ones.

> 

sol1 := sol[9]

{P(y, t) = c__1, Q(y, t) = f__1(t), A[0](y, t) = (1/8)*f__3(y)*(4*h[2]*h[4]-h[3]^2)/h[4]^(3/2), A[1](y, t) = (1/2)*h[3]*f__3(y)/h[4]^(1/2), A[2](y, t) = f__3(y)*h[4]^(1/2), B[1](y, t) = f__3(y), B[2](y, t) = 0, a[0](y, t) = -(1/4)*(c__1^2*h[3]+2*h[4]^(1/2)*(diff(f__1(t), t)))/(c__1*h[4]^(1/2)), a[1](y, t) = -h[4]^(1/2)*c__1, b[1](y, t) = 0}

We choose c__1,  f__3 f__1 as follows, which gives agreement with Eq 25 except for Q, which is not quite the same.
I interpret K__2[t](t) in Eq 25 as diff(K__2(t),t).

> 

sol2 := eval(sol1, {c__1 = P, f__1(t) = K__2(t), f__3(y) = -(1/2)*P*sqrt(h[4])*K__1(y)})

{P(y, t) = P, Q(y, t) = K__2(t), A[0](y, t) = -(1/16)*P*K__1(y)*(4*h[2]*h[4]-h[3]^2)/h[4], A[1](y, t) = -(1/4)*h[3]*P*K__1(y), A[2](y, t) = -(1/2)*P*h[4]*K__1(y), B[1](y, t) = -(1/2)*P*h[4]^(1/2)*K__1(y), B[2](y, t) = 0, a[0](y, t) = -(1/4)*(P^2*h[3]+2*h[4]^(1/2)*(diff(K__2(t), t)))/(P*h[4]^(1/2)), a[1](y, t) = -h[4]^(1/2)*P, b[1](y, t) = 0}

Which means xi is

> 

eqxi2 := eval(eqxi, sol2)

xi = P*x+K__2(t)

Solutions as a function of phi(xi)

> 

equxi := eval(Seriesu, sol2); eqvxi := eval(Seriesv, sol2)

-(1/4)*(P^2*h[3]+2*h[4]^(1/2)*(diff(K__2(t), t)))/(P*h[4]^(1/2))-h[4]^(1/2)*P*phi(xi)

-(1/16)*P*K__1(y)*(4*h[2]*h[4]-h[3]^2)/h[4]-(1/4)*h[3]*P*K__1(y)*phi(xi)-(1/2)*P*h[4]^(1/2)*K__1(y)*(h[0]+h[1]*phi(xi)+h[2]*phi(xi)^2+h[3]*phi(xi)^3+h[4]*phi(xi)^4)^(1/2)-(1/2)*P*h[4]*K__1(y)*phi(xi)^2

Top p 194

> 

case1 := {h[0] = r^2, h[1] = 2*r*p, h[2] = p^2+2*q*r, h[3] = 2*p*q, h[4] = q^2}

{h[0] = r^2, h[1] = 2*r*p, h[2] = p^2+2*q*r, h[3] = 2*p*q, h[4] = q^2}

> 

ode1 := `assuming`([simplify(eval(eq3, case1))], [q*phi(xi)^2+p*phi(xi)+r > 0])

diff(phi(xi), xi) = q*phi(xi)^2+p*phi(xi)+r

> 

S1 := phi(xi) = -(p+sqrt(p^2-4*q*r)*tanh((1/2)*sqrt(p^2-4*q*r)*xi))/(2*q)

phi(xi) = -(1/2)*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*xi))/q

> 

`assuming`([simplify(odetest(S1, ode1))], [p^2-4*q*r > 0])

0

> 

final := {u(x, y, t) = eval(eval(equxi, `union`({S1}, case1)), eqxi2), v(x, y, t) = eval(eval(eqvxi, `union`({S1}, case1)), eqxi2)}

{u(x, y, t) = -(1/4)*(2*P^2*p*q+2*(q^2)^(1/2)*(diff(K__2(t), t)))/(P*(q^2)^(1/2))+(1/2)*(q^2)^(1/2)*P*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))/q, v(x, y, t) = -(1/16)*P*K__1(y)*(4*q^2*(p^2+2*q*r)-4*p^2*q^2)/q^2+(1/4)*p*P*K__1(y)*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))-(1/2)*P*(q^2)^(1/2)*K__1(y)*(r^2-r*p*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))/q+(1/4)*(p^2+2*q*r)*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))^2/q^2-(1/4)*p*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))^3/q^2+(1/16)*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))^4/q^2)^(1/2)-(1/8)*P*K__1(y)*(p+(p^2-4*q*r)^(1/2)*tanh((1/2)*(p^2-4*q*r)^(1/2)*(P*x+K__2(t))))^2}

> 

`assuming`([simplify(pdetest(final, [pde1, pde2]))], [q > 0, p^2-4*q*r > 0])

[0, 0]

Try the xi in the paper

> 

final2 := subs(K__2(t) = K__2(t)+Int(K__1(y), y), final)

Doesn't work - following is not zero, though perhaps we are missing some assumptions.

> 

`assuming`([simplify(pdetest(final2, [pde1, pde2]), symbolic)], [q > 0, p^2-4*q*r > 0])
> 

NULL

Download F4.mw

One way...

dharr 9152
March 26 2026
3 1

Here's one way.

> 

restart;
> 

S:=RealRange(-2*Pi,Open(-7/4*Pi)), RealRange(Open(-5/4*Pi),Open(-3/4*Pi)),
RealRange(Open(-1/4*Pi),Open(1/4*Pi)), RealRange(Open(3/4*Pi),Open(5/4*Pi)),
RealRange(Open(7/4*Pi),2*Pi);

RealRange(-2*Pi, Open(-(7/4)*Pi)), RealRange(Open(-(5/4)*Pi), Open(-(3/4)*Pi)), RealRange(Open(-(1/4)*Pi), Open((1/4)*Pi)), RealRange(Open((3/4)*Pi), Open((5/4)*Pi)), RealRange(Open((7/4)*Pi), 2*Pi)

> 

W:=[eval(S,{Open=(()->_passed),RealRange=`[]`})];

[[-2*Pi, -(7/4)*Pi], [-(5/4)*Pi, -(3/4)*Pi], [-(1/4)*Pi, (1/4)*Pi], [(3/4)*Pi, (5/4)*Pi], [(7/4)*Pi, 2*Pi]]

> 

lprint(W);

[[-2*Pi, -7/4*Pi], [-5/4*Pi, -3/4*Pi], [-1/4*Pi, 1/4*Pi], [3/4*Pi, 5/4*Pi], [7/
4*Pi, 2*Pi]]

> 

NULL
> 

NULL

Download RealRange.mw

DocumentTools:-Tabulate...

dharr 9152
March 20 2026
3 2

For 50% width try

DocumentTools:-Tabulate([p1,p2],width=50):

convert to differential equation...

dharr 9152
March 18 2026
1 0

The simplest method of solution if the indefinite integral is not known is to convert to a differential equation and solve that. As for approximating by the form of equation you want, that is difficult because of the tendency of NLPSolve to encounter non-numeric or complex values. What about series solutions, Pade approximants, minimax etc; would those be of interest?

https://www.mapleprimes.com/questions/242328-Strange-Integral-Equation

We want to find x(t) that solves int((3*x(t)+5)*(diff(x(t), t)), t = 1 .. x(t)) = 2*t. Let's be careful about using a different symbol for the dummy integration variable (X) than the limits (x). "Cancelling dt's" (changing variables from t to x")" and assuming the initial condition is x=1 at t=0

x

> 

restart
> 

q := Int(3*X+5, X = 1 .. x(t)); value(%) = 2*t; solve(%, x(t)); x := MakeFunction(%[1], t)

Int(3*X+5, X = 1 .. x(t))

(3/2)*x(t)^2-13/2+5*x(t) = 2*t

-5/3+(2/3)*(16+3*t)^(1/2), -5/3-(2/3)*(16+3*t)^(1/2)

proc (t) options operator, arrow; -5/3+(2/3)*(16+3*t)^(1/2) end proc

> 

simplify(x(0))

1

> 

restart

Now suppose we do not know the indefinite integral - we can differentiate and solve the differential equation (numerically if there is no analytical solution)

> 

inteq := Int(3*X+5, X = 1 .. x(t)) = 2*t; de := diff(%, t); exact := rhs(dsolve({%, x(0) = 1}, x(t)))

Int(3*X+5, X = 1 .. x(t)) = 2*t

(diff(x(t), t))*(3*x(t)+5) = 2

-5/3+(2/3)*(16+3*t)^(1/2)

Suppose we seek an approximate solution in the following form that is good over the range t=0..tmax

> 

x_approx := a*t^(1/5)+b*t^(1/3)+c*ln(t)+d*t

a*t^(1/5)+b*t^(1/3)+c*ln(t)+d*t

We want to minimize the following (take tmax=5)

> 

tmax := 5; Int((eval(lhs(inteq), x(t) = x_approx)-rhs(inteq))^2, t = 0 .. tmax); resid := MakeFunction(%, [a, b, c, d])

5

Int((Int(3*X+5, X = 1 .. a*t^(1/5)+b*t^(1/3)+c*ln(t)+d*t)-2*t)^2, t = 0 .. 5)

proc (a, b, c, d) options operator, arrow; Int((Int(3*X+5, X = 1 .. a*t^(1/5)+b*t^(1/3)+c*ln(t)+d*t)-2*t)^2, t = 0 .. 5) end proc

> 

Optimization:-NLPSolve(resid, initialpoint = [.1, .1, .5, 1], method = nonlinearsimplex)

Error, (in Optimization:-NLPSolve) complex value encountered

> 

NULL

Download intsolve.mw

attempt...

dharr 9152
March 10 2026
2 9

Eq, 48 looks wrong to me but maybe it is a sign issue earlier. Need a better paper.

[Edit -varpi^2 missing from the second Eq. 49 added, but still wrong.]

> 

restart
> 

with(PDEtools)
> 

undeclare(prime, quiet)
> 

declare(u(x, t), quiet); declare(U(xi), quiet); declare(V(xi), quiet); declare(Omega(x, t), quiet); declare(Q(x, t), quiet); declare(R(x, t), quiet); declare(CHI(x, t), quiet); declare(Z(x, t), quiet)
> 

pde1 := diff(Omega(x, t), `$`(t, 2))-(diff(Omega(x, t), `$`(x, 2)))-alpha^2*Omega(x, t)+beta^2*Omega(x, t)^3
> 

G1 := U(-epsilon*t+x) = U(xi); G2 := (D(U))(-epsilon*t+x) = diff(U(xi), xi); G3 := ((D@@2)(U))(-epsilon*t+x) = diff(U(xi), `$`(xi, 2)); G4 := ((D@@3)(U))(-epsilon*t+x) = diff(U(xi), `$`(xi, 3)); G5 := ((D@@4)(U))(-epsilon*t+x) = diff(U(xi), `$`(xi, 4)); G6 := ((D@@5)(U))(-epsilon*t+x) = diff(U(xi), `$`(xi, 5))
> 

T := xi = -epsilon*t+x; T1 := Omega(x, t) = U(-epsilon*t+x)
> 

P1 := diff(Omega(x, t), `$`(t, 2))-(diff(Omega(x, t), `$`(x, 2)))-alpha^2*Omega(x, t)+beta^2*Omega(x, t)^3
> 

P11 := eval(P1, {T, T1})
> 

P111 := subs({G1, G2, G3, G4, G5, G6}, P11)
> 

pde1 := P111 = 0
> 

D1 := numer(lhs(pde1))*denom(rhs(pde1)) = numer(rhs(pde1))*denom(lhs(pde1)); simplify(%, 'symbolic'); ode := collect(%, {U(xi), diff(diff(U(xi), xi), xi)})

beta^2*U(xi)^3-alpha^2*U(xi)+(epsilon^2-1)*(diff(diff(U(xi), xi), xi)) = 0

Eq. 37

> 

ode1 := expand(lhs(ode)/(alpha^2*beta^2))

(diff(diff(U(xi), xi), xi))*epsilon^2/(alpha^2*beta^2)-(diff(diff(U(xi), xi), xi))/(alpha^2*beta^2)-U(xi)/beta^2+U(xi)^3/alpha^2

Note that each term in the Hamiltonian integrand when differentiated wrt U(xi) gives a term in ode. (I have no idea why we do this.) So differentiating wrt xi would give a term in the ode multiplied by diff(U(xi),xi) (using the chain rule). So we can set up a differential equation for the integrand f(xi)

> 

def := diff(f(xi), xi) = ode1*(diff(U(xi), xi))

diff(f(xi), xi) = ((diff(diff(U(xi), xi), xi))*epsilon^2/(alpha^2*beta^2)-(diff(diff(U(xi), xi), xi))/(alpha^2*beta^2)-U(xi)/beta^2+U(xi)^3/alpha^2)*(diff(U(xi), xi))

We want the first integral of this

> 

DEtools:-firint(def, f(xi)); integrand := collect(solve(eval(%, {_C1 = 0, c__1 = 0}), f(xi)), diff, simplify)

f(xi)-((1/4)*beta^2*U(xi)^4+(1/2)*(epsilon^2-1)*(diff(U(xi), xi))^2-(1/2)*alpha^2*U(xi)^2)/(alpha^2*beta^2)+_C1 = 0

(1/2)*(epsilon^2-1)*(diff(U(xi), xi))^2/(alpha^2*beta^2)+(1/4)*U(xi)^2*(beta^2*U(xi)^2-2*alpha^2)/(alpha^2*beta^2)

From which we can extract chi and delta, Eq 41.

> 

chi := op(1, integrand); delta := expand(-op(2, integrand))

(1/2)*(epsilon^2-1)*(diff(U(xi), xi))^2/(alpha^2*beta^2)

-(1/4)*U(xi)^4/alpha^2+(1/2)*U(xi)^2/beta^2

Eq. 42.

> 

H := chi+delta

(1/2)*(epsilon^2-1)*(diff(U(xi), xi))^2/(alpha^2*beta^2)-(1/4)*U(xi)^4/alpha^2+(1/2)*U(xi)^2/beta^2

Eq. 43

> 

sol1 := U(xi) = varpi*cos(phi*xi)

U(xi) = varpi*cos(phi*xi)

Eq. 45 and 46

> 

eval(convert(H, D), xi = 0); H0 := eval(%, {U(0) = varpi, (D(U))(0) = 0})

(1/2)*(epsilon^2-1)*(D(U))(0)^2/(alpha^2*beta^2)-(1/4)*U(0)^4/alpha^2+(1/2)*U(0)^2/beta^2

-(1/4)*varpi^4/alpha^2+(1/2)*varpi^2/beta^2

Eq. 47.

> 

H1 := eval(H, sol1)-H0

(1/2)*(epsilon^2-1)*varpi^2*phi^2*sin(phi*xi)^2/(alpha^2*beta^2)-(1/4)*varpi^4*cos(phi*xi)^4/alpha^2+(1/2)*varpi^2*cos(phi*xi)^2/beta^2+(1/4)*varpi^4/alpha^2-(1/2)*varpi^2/beta^2

Eq. 48. Something wrong with Eq. 48

> 

H2 := collect(expand(algsubs(phi*xi = (1/4)*Pi, H1)), phi)

((1/4)*varpi^2*epsilon^2/(alpha^2*beta^2)-(1/4)*varpi^2/(alpha^2*beta^2))*phi^2+(3/16)*varpi^4/alpha^2-(1/4)*varpi^2/beta^2

Eq. 49 also is not right?

> 

eqphi := phi = sort([solve(H2, phi)])[2]

phi = (1/2)*((epsilon^2-1)*(-3*beta^2*varpi^2+4*alpha^2))^(1/2)/(epsilon^2-1)

> 

sol2 := eval(sol1, eqphi)

U(xi) = varpi*cos((1/2)*((epsilon^2-1)*(-3*beta^2*varpi^2+4*alpha^2))^(1/2)*xi/(epsilon^2-1))

> 

simplify(odetest(sol2, ode))

(1/4)*varpi*cos(((epsilon^2-1)*(-3*beta^2*varpi^2+4*alpha^2))^(1/2)*xi/(2*epsilon^2-2))*(4*beta^2*varpi^2*cos(((epsilon^2-1)*(-3*beta^2*varpi^2+4*alpha^2))^(1/2)*xi/(2*epsilon^2-2))^2+3*varpi^2*beta^2-8*alpha^2)

What if we adjust to agree with Eq. 49.

> 

eqphi2 := phi = sqrt((-2*alpha^2+beta^2)*`&varpi;`^2)/(sqrt(2)*sqrt(`&epsilon;`^2-1))

phi = (1/2)*((-2*alpha^2+beta^2)*varpi^2)^(1/2)*2^(1/2)/(epsilon^2-1)^(1/2)

> 

sol3 := eval(sol1, eqphi2)

U(xi) = varpi*cos((1/2)*((-2*alpha^2+beta^2)*varpi^2)^(1/2)*2^(1/2)*xi/(epsilon^2-1)^(1/2))

Still not right.

> 

simplify(odetest(sol3, ode))

cos((1/2)*((-2*alpha^2+beta^2)*varpi^2)^(1/2)*2^(1/2)*xi/(epsilon^2-1)^(1/2))*((alpha^2+beta^2*(cos((1/2)*((-2*alpha^2+beta^2)*varpi^2)^(1/2)*2^(1/2)*xi/(epsilon^2-1)^(1/2))^2-1/2))*varpi^2-alpha^2)*varpi

> 

NULL

NULL

Download f-s.mw

complex analysis...

dharr 9152
March 05 2026
3 1

Maple's default behavior is to work in the complex domain, so working with series, polynomials etc is done without any additional commands. For contour integrals, you can usually just use the residue theorem. Of course for a simple contour like a circle around the origin, you can manually convert it to a normal integral in r,theta. Standard Maple does not have any packages for complex analysis, though there might be something in the application center. Useful commands are

singular

residue

series

numapprox:-laurent (just a special case of series).

using pds:-value...

dharr 9152
March 02 2026
3 3

Use

uu := pds:-value(u(x,t));

This returns a procedure you can give x and t values to. For example,

uu(0.0001,0.1);

gives

value.mw

1 2 3 4 5 6 7 Last Page 1 of 87