dharr

Dr. David Harrington

9102 Reputation

22 Badges

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

Social Networks and Content at Maplesoft.com

Maple Application Center
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

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.

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 arc length is 50 m. The 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 minimize

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 := factor(sort([solve(ode4, diff(Y(X), X))])[2]); sol4 := Y(X) = Int(subs(X = x, ans), x = 0 .. X)

-(X^2-2*mu-1)/(-(X-1)*(X+1)*(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)

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-1)*(X+1)*(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)

-(X-1)*(X+1)*(W-4*mu-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)

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, 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-1/2)*EllipticK(1/(4*mu+1)^(1/2))*mu)/(-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-1/2)*EllipticK(1/(4*mu+1)^(1/2))*mu)/(-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-1/2)*EllipticK(1/(4*mu+1)^(1/2))*mu)/(mu^3*EllipticK(1/(4*mu+1)^(1/2))^3)

Now we can maximize this wrt mu.

plot(cost2, mu = -1 .. -.25)

 

Differentiate - we see that mu = -1/4 is a solution but there is a singularity there, so it is not the solution we want.

dcdmu := simplify(numer(normal(diff(cost2, mu))))

mu = -1/4

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

Sadly, after getting an exact solution this far, we cannot get a symbolic solution for the cost

solve(dcdmu/(mu+1/4), mu)

-(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.-4.514161256*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.-4.514161256*I)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.+.6520869849*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.+.6520869849*I)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.+4.514161256*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, -0.+4.514161256*I)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.-4.514161256*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.-4.514161256*I)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.-.6520869849*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.-.6520869849*I)^2, -(1/4)*(RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.+4.514161256*I)^2-1)/RootOf(-5*EllipticK(_Z)^2*_Z^4+2*EllipticK(_Z)*EllipticE(_Z)*_Z^4+6*EllipticK(_Z)^2*_Z^2-12*EllipticK(_Z)*EllipticE(_Z)*_Z^2+5*EllipticE(_Z)^2*_Z^2-EllipticK(_Z)^2+2*EllipticK(_Z)*EllipticE(_Z)-EllipticE(_Z)^2, 0.+4.514161256*I)^2

So we settle for a numeric one

mumax := fsolve(dcdmu/(mu+1/4), mu = -.4 .. -.25)

-.2622683419

And the final cost is

finalcost := simplify(eval(cost2, mu = mumax), zero)

22665.80623

Now we can find x__max

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

33.20348480

Check fence length

eval(eqfence, {mu = mumax, x__max = xmaxnum})

50 = 50.00000002-0.*I

In real units the fence curve is

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

-149.8858841*(-0.9070530060e-3*x^2+1)^(1/2)*(0.9070530060e-3*x^2+0.49073368e-1)^(1/2)*(-.5245366838*EllipticF(0.3011732070e-1*x, 4.514161239*I)+0.49073368e-1*EllipticE(0.3011732070e-1*x, 4.514161239*I))/(-0.8227451557e-6*x^4+0.8625408600e-3*x^2+0.49073368e-1)^(1/2)

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

NULL

NULL

Download meadow2.mw

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.

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

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

For 50% width try

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

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

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

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).

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

Here's one way

plot(t,'title'=InertForm:-Typeset(the_title),size=[300,300]);

You can use SimilarityTransformation directly without using Eq 29. For the second paper, I don't know about prolongations, which seem to be done with Eta_k. I'm guessing by choosing the arbitrary functions appropriately you can get what you want.

Download Lie.mw

There are 510 paths, which agrees with the formula in your linked page. The approach, which works for any paths, not just Hamiltonian ones, here uses "multiplication" of modified adjacency matrices in a similar fashion to this post. However in this case I track the vertices using the Bits package, with one bit per vertex, and kill a walk whenever it revisits a vertex, and I hard coded the modified matrix multiplication rather than use LinearAlgebra:-Generic. (I'm short on time or I would describe the algorithm in more detail.)

The code can output the vertices visited (coded as an integer) but not their order, which is useless in this case, where all paths visit the same vertices. (For a DAG this would be enough to reconstruct the paths.)

Procedure Paths in startup code region

restart

with(GraphTheory)

G := SpecialGraphs:-GridGraph(11, 3); n := NumberOfVertices(G)

33

DrawGraph(G)

Paths(G, "2,2", "10,2")[n-1]

510

 

Download Paths.mw

Paths:=proc(G0::Graph,V1,V2,{minlength::posint := 1, maxlength::posint := GraphTheory:-NumberOfVertices(G0)-1, output_paths::truefalse := false})
	uses GraphTheory;
	local k, A, Adj, B, n, i, j, SetBit, verts, v1, v2, row, col, mrow, V, s, vec, paths, nrows, G, VerticesG,
		terminals_done, terminals, outdegs, killarcs, new_terminals;
	G := CopyGraph(G0);
	VerticesG := Vertices(G);
	n := NumberOfVertices(G);
	# output table records number of paths (output_paths = false) or Vectors of paths as coded integers (output_paths = true)  
	paths := table('sparse'):
	if _params['V1'] = NULL and _params['V2'] = NULL then
		verts := false;
	elif _params['V1'] = NULL or _params['V2'] = NULL then
	     error "must be zero or two vertices specified"
	else
		if not member(V1,VerticesG,'v1') then error "Vertex %1 not in graph",V1 end if;
		if not member(V2,VerticesG,'v2') then error "Vertex %1 not in graph",V2 end if;
		verts:=true;
	end if;
	# if we are looking for paths between vertices on a digraph, remove extraneous pieces of the graph
	if verts and IsDirected(G) then
		if OutDegree(G, V1) = 0 then error "start vertex has no departures" end if;
		if InDegree(G, V2) = 0 then error "end vertex has no arrivals" end if;
		userinfo(2, procname, "pruning extra vertices from digraph");
		terminals_done:={V1}; # can't delete V1
		do 
			outdegs:=map2(OutDegree, G, VerticesG); # includes isolated vertices
			terminals := {remove(`=`, VerticesG[[ListTools:-SearchAll(0, outdegs)]], V2)[]};
			new_terminals := terminals minus terminals_done;
			if nops(new_terminals) = 0 then break end if;
			killarcs := map(t -> seq([i, t], i in Arrivals(G, t)), new_terminals);
   			DeleteArc(G, killarcs);
			terminals_done := terminals_done union terminals;
  		end do;
  		terminals_done := terminals_done minus {V1};
  		# and delete departures from last vertex so Matrix is zero at the end
  		if OutDegree(G, V2) = 1 then
  			DeleteArc(G, [V2, Departures(G,V2)[]])
  		elif OutDegree(G, V2) > 1 then
  			DeleteArc(G, [seq([V2, i], i in Departures(G,V2))])
  		end if;
  		# update graph, n, v1, v2, VerticesG
		G := DeleteVertex(G, terminals_done);
		userinfo(2, procname, "%1 extra vertices deleted", nops(terminals_done));
		VerticesG := Vertices(G);
		n := NumberOfVertices(G);
		member(V1,VerticesG,'v1');
		member(V2,VerticesG,'v2');
	end if;
	# set up adjacency matrices.
	Adj:=AdjacencyMatrix(G);
	SetBit := p -> integermul2exp(1, p - 1);
	A := Matrix(n, n, (i, j) -> if Adj[i, j] = 1 then SetBit(j); else 0; end if);
	if verts then # for B we only need the row corresponding to the first vertex
		nrows:=1;
		B :=Matrix(1, n, 'storage' = 'sparse', order = 'C_order');
		for j to n do
    			if Adj[v1, j] = 1 then B[1, j] := Vector([Bits:-Or(SetBit(v1), SetBit(j))]); end if;
    		end do;
    		if not output_paths then
    			paths[1]:= Adj[v1, v2];		
    		elif minlength = 1 then
    			paths[1] := B[1, v2] # 0 or a vector
    		end if  			
    	else
    		nrows := n;
		B := Matrix(nrows, n, 'storage' = 'sparse', order = 'C_order');
		for i to n do
			for j to n do
    				if Adj[i, j] = 1 then B[i, j] := Vector([Bits:-Or(SetBit(i), SetBit(j))]); end if;
    			end do
    		end do;
    		if not output_paths then
    			paths[1] := add(Adj);
    		elif minlength = 1 then
    			paths[1] := Vector([entries(B,'nolist')]);
    		end if
    	end if;
	# do the matrix multiplications and collect the paths of each length
	for k from 2 to n-1 do   
		for row to nrows do # rows of B diving down ...
			mrow := B[row]; # (save row so can do in-place)
			for col to n do # ... cols of A
				V := Vector():
				s := 0;
				for i, vec in mrow do
					if A[i, col] <> 0 then
						for j to numelems(vec) do
							if Bits:-And(vec[j], A[i,col]) = 0 then # not a repeat vertex
								s := s+1;
								V(s) := Bits:-Or(vec[j], A[i,col]); # add the vertex
							end if
						end do
					end if
				end do:
				if numelems(V) = 0 then
					B[row,col] := 0;
				else	 
					B[row,col] := V;
				end if;
			end do
		end do;
		if andmap(`=`,B,0) then break end if;
		if verts then
			if not output_paths then
				paths[k] := ifelse(B[1, v2] = 0, 0, numelems(B[1, v2]));
			elif k >= minlength and k <= maxlength then
				paths[k] := B[1,v2];
			end if;
		else
			if not output_paths then
				s := 0;
				for i, vec in B do
					s + numelems(vec)
				end do;
				paths[k] := s;
			elif k >= minlength and k <= maxlength then
				paths[k] := Vector([entries(B,'nolist')]);
			end if;	
		end if;
	end do;
	eval(paths)
end proc;

Here's one way.

results2 := Matrix(map2(eval, [N, tbo, two], [seq(results)]));

PrimesQuestion.mw

The following URL finds my Orbitals package

https://www.maplesoft.com/Applications/Detail.aspx?id=4865

but substituting your ID doesn't find anything, suggesting that the ID is no longer valid or the application is not there.

It is much easier to use a code edit region. But you can do it your intended way with EvalString instead of Run. GetVariable returns the variable to Maple. Here are both ways:

python.mw

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