Advanced Search

dharr

Dr. David Harrington

8355 Reputation

22 Badges

21 years, 8 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

various small things...

dharr 8355
August 30 2025
2 1

For conditions with <>, just leave them out. They are just assumptions that will prevent the denominator being zero.

Do the equation assumptions as substitutions as you did, but then you do not need them as assumptions. Minor corrections noted on the worksheet.

ode-17.mw

 

Painleve for coupled equations...

dharr 8355
August 28 2025
3 4

Here is how to do coupled equations for the Hirota-Satsuma system. The reference seems the clearest about the algorithm, though I didn't understand how they decided on the a=1/2 condition. For the Schroedinger equation there were equations and not simple integers for the resonant points and I got many, but not all of the results - see schrodinger-test.mw.

Painleve analysis following the notation in doi: 10.2991/jnmp.2006.13.1.8

Hirota-Satsuma system example from that paper.

> 

restart
> 

with(PDEtools); with(LinearAlgebra)
> 

undeclare(prime, quiet)
> 

vars := x, t; declare(u(vars), quiet); declare(v(vars), quiet); declare(g(vars), quiet); declare(U(vars), quiet); declare(V(vars), quiet)
> 

pde1 := a*(6*u(x, t)*(diff(u(x, t), x))+diff(u(x, t), x, x, x))-2*v(x, t)*(diff(v(x, t), x))-(diff(u(x, t), t)); pde2 := -3*u(x, t)*(diff(v(x, t), x))-(diff(v(x, t), `$`(x, 3)))-(diff(v(x, t), t))

a*(6*u(x, t)*(diff(u(x, t), x))+diff(diff(diff(u(x, t), x), x), x))-2*v(x, t)*(diff(v(x, t), x))-(diff(u(x, t), t))

-3*u(x, t)*(diff(v(x, t), x))-(diff(diff(diff(v(x, t), x), x), x))-(diff(v(x, t), t))

> 

depvars := [u(vars), v(vars), g(vars)]

[u(x, t), v(x, t), g(x, t)]

> 

findexp := proc (T) options operator, arrow; simplify(g(vars)*Physics:-diff(T, g(vars))/T) end proc

Step 1

> 

eq31u := u(vars) = chi[1]*g(vars)^alpha[1]; eq31v := v(vars) = chi[2]*g(vars)^alpha[2]

u(x, t) = chi[1]*g(x, t)^alpha[1]

v(x, t) = chi[2]*g(x, t)^alpha[2]

> 

eqG1 := expand(eval(pde1, {eq31u, eq31v})); eqG2 := expand(eval(pde2, {eq31u, eq31v}))

6*a*chi[1]^2*(g(x, t)^alpha[1])^2*alpha[1]*(diff(g(x, t), x))/g(x, t)+a*chi[1]*g(x, t)^alpha[1]*alpha[1]^3*(diff(g(x, t), x))^3/g(x, t)^3+3*a*chi[1]*g(x, t)^alpha[1]*alpha[1]^2*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^2-3*a*chi[1]*g(x, t)^alpha[1]*alpha[1]^2*(diff(g(x, t), x))^3/g(x, t)^3+a*chi[1]*g(x, t)^alpha[1]*alpha[1]*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)-3*a*chi[1]*g(x, t)^alpha[1]*alpha[1]*(diff(diff(g(x, t), x), x))*(diff(g(x, t), x))/g(x, t)^2+2*a*chi[1]*g(x, t)^alpha[1]*alpha[1]*(diff(g(x, t), x))^3/g(x, t)^3-2*chi[2]^2*(g(x, t)^alpha[2])^2*alpha[2]*(diff(g(x, t), x))/g(x, t)-chi[1]*g(x, t)^alpha[1]*alpha[1]*(diff(g(x, t), t))/g(x, t)

-3*chi[1]*g(x, t)^alpha[1]*chi[2]*g(x, t)^alpha[2]*alpha[2]*(diff(g(x, t), x))/g(x, t)-chi[2]*g(x, t)^alpha[2]*alpha[2]^3*(diff(g(x, t), x))^3/g(x, t)^3-3*chi[2]*g(x, t)^alpha[2]*alpha[2]^2*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^2+3*chi[2]*g(x, t)^alpha[2]*alpha[2]^2*(diff(g(x, t), x))^3/g(x, t)^3-chi[2]*g(x, t)^alpha[2]*alpha[2]*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)+3*chi[2]*g(x, t)^alpha[2]*alpha[2]*(diff(diff(g(x, t), x), x))*(diff(g(x, t), x))/g(x, t)^2-2*chi[2]*g(x, t)^alpha[2]*alpha[2]*(diff(g(x, t), x))^3/g(x, t)^3-chi[2]*g(x, t)^alpha[2]*alpha[2]*(diff(g(x, t), t))/g(x, t)

Find all the distinct powers of g for eqG1 and eqG2

> 

exps1 := [(`minus`(map(findexp, {op(eqG1+forceplus)}), {0}))[]]; exps2 := [(`minus`(map(findexp, {op(eqG2+forceplus)}), {0}))[]]

[-2+alpha[1], alpha[1]-3, alpha[1]-1, 2*alpha[1]-1, 2*alpha[2]-1]

[-2+alpha[2], alpha[2]-3, alpha[2]-1, alpha[1]+alpha[2]-1]

Find the slopes wrt alpha[1] and alpha[2] and separate into categories based on the different pairs of slopes. Then find the lowest lines for each case.

> 

mins1 := `~`[min](convert(ListTools:-Classify(proc (q) options operator, arrow; [diff(q, alpha[1]), diff(q, alpha[2])] end proc, exps1), list)); mins2 := `~`[min](convert(ListTools:-Classify(proc (q) options operator, arrow; [diff(q, alpha[1]), diff(q, alpha[2])] end proc, exps2), list))

[2*alpha[2]-1, 2*alpha[1]-1, alpha[1]-3]

[alpha[1]+alpha[2]-1, alpha[2]-3]

Find all intersection points between pairs of these. Each is a list of a set of the alpha values and then the corresponding exponent.

> 

"intersections1:={seq(seq([(ans:=solve(mins1[i]=mins1[j],{alpha[1],alpha[2]})),eval(mins1[i],ans)],i=1..j-1),j=1..nops(mins1))}; intersections2:={seq(seq([(ans:=solve(mins2[i]=mins2[j],{alpha[1],alpha[2]})),eval(mins2[i],ans)],i=1..j-1),j=1..nops(mins2))};"

{[{alpha[1] = -2, alpha[2] = alpha[2]}, -5], [{alpha[1] = alpha[2], alpha[2] = alpha[2]}, 2*alpha[2]-1], [{alpha[1] = 2*alpha[2]+2, alpha[2] = alpha[2]}, 2*alpha[2]-1]}

{[{alpha[1] = -2, alpha[2] = alpha[2]}, alpha[2]-3]}

When an alpha is undetermined, then the algoritm says put successive integers in and accept a possibility (branch) if each equation has at least two leading (most negative) terms.
There is only one solution to the second case, which sets alpha[1] = -2 but doesn't fix alpha[2]. This is also the first solution in the first case, which corresponds to leading order -5.

If we put it into the second solution in the first case, i.e. alpha[1]=alpha[2]=-2, then the order is 2*alpha[2] - 1 = -5.
If we put it into the third solution in the first case, i.e. alpha[1]=2+2*alpha[2], or -2 = 2+2*alpha[2] or alpha[2] = -2 and the leading order is 2*alpha[2] - 1 = -5.
So alpha[1] = alpha[2] = -2 will work in any case.
Note that alpha[2]>-2 will also work with -5 (we need at least two in each equation having -5)

> 

eval(mins1, {alpha[1] = -2, alpha[2] = -1}); eval(mins2, {alpha[1] = -2, alpha[2] = -1})

[-3, -5, -5]

[-4, -4]

But alpha[2] <-2 won't work - here we have only one of leading order -7 (we have two -5s, but these aren't the most negative).

> 

eval(mins1, {alpha[1] = -2, alpha[2] = -3}); eval(mins2, {alpha[1] = -2, alpha[2] = -3})

[-7, -5, -5]

[-6, -6]


Apparently, alpha[1]=-2, alpha[2]>=0 leads to u[0] or v[0] identically 0, so is excluded. So we have two possibilities (need integers or have to do a special analysis):

> 

b1 := {alpha[1] = -2, alpha[2] = -2}; b2 := {alpha[1] = -2, alpha[2] = -1}

{alpha[1] = -2, alpha[2] = -2}

{alpha[1] = -2, alpha[2] = -1}

Start with the b1 case

> 

eqb := b1

{alpha[1] = -2, alpha[2] = -2}

> 

eq32b := eval({u(vars) = chi[1]*g(vars)^alpha[1], v(vars) = chi[2]*g(vars)^alpha[2]}, eqb)

{u(x, t) = chi[1]/g(x, t)^2, v(x, t) = chi[2]/g(x, t)^2}

 

> 

eq32b := expand(eval([pde1, pde2], eq32b))

[-12*a*chi[1]^2*(diff(g(x, t), x))/g(x, t)^5-24*a*chi[1]*(diff(g(x, t), x))^3/g(x, t)^5+18*a*chi[1]*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4-2*a*chi[1]*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3+4*chi[2]^2*(diff(g(x, t), x))/g(x, t)^5+2*chi[1]*(diff(g(x, t), t))/g(x, t)^3, 6*chi[1]*chi[2]*(diff(g(x, t), x))/g(x, t)^5+24*chi[2]*(diff(g(x, t), x))^3/g(x, t)^5-18*chi[2]*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4+2*chi[2]*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3+2*chi[2]*(diff(g(x, t), t))/g(x, t)^3]

We want only the terms with g^(mindeg), but we cannot collect wrt g, so convert it to a simple name before collecting

> 

eq32c := subs(g[] = G, ToJet(eq32b, depvars, jetnotation = jetvariableswithbrackets))

[-12*a*chi[1]^2*g[x]/G^5-24*a*chi[1]*g[x]^3/G^5+18*a*chi[1]*g[x]*g[x, x]/G^4-2*a*chi[1]*g[x, x, x]/G^3+4*chi[2]^2*g[x]/G^5+2*chi[1]*g[t]/G^3, 6*chi[1]*chi[2]*g[x]/G^5+24*chi[2]*g[x]^3/G^5-18*chi[2]*g[x]*g[x, x]/G^4+2*chi[2]*g[x, x, x]/G^3+2*chi[2]*g[t]/G^3]

> 

`~`[ldegree](eq32c, G); mindeg1, mindeg2 := %[]

[-5, -5]

-5, -5

> 

S1 := FromJet(coeff(collect(eq32c[1], G), G, mindeg1), depvars); S2 := FromJet(coeff(collect(eq32c[2], G), G, mindeg1), depvars)

-24*a*chi[1]*(diff(g(x, t), x))^3-12*a*chi[1]^2*(diff(g(x, t), x))+4*chi[2]^2*(diff(g(x, t), x))

24*chi[2]*(diff(g(x, t), x))^3+6*chi[1]*chi[2]*(diff(g(x, t), x))

> 

solve({S1, S2, chi[1] <> 0, chi[2] <> 0}, {chi[1], chi[2]}, explicit); eqchi := %[1]

{chi[1] = -4*(diff(g(x, t), x))^2, chi[2] = 2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2}, {chi[1] = -4*(diff(g(x, t), x))^2, chi[2] = -2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2}

{chi[1] = -4*(diff(g(x, t), x))^2, chi[2] = 2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2}

Following are u[0]/g^alpha[1] and v[0]/g^alpha[2], which we use in step 2 - Eq 3.19 (1)

> 

eq19ub := rhs(eval(eq31u, `union`(eqchi, eqb))); eq19vb := rhs(eval(eq31v, `union`(eqchi, eqb)))

-4*(diff(g(x, t), x))^2/g(x, t)^2

2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2/g(x, t)^2

So we found (chi) and (alpha) which was our target for first step

Step 2 - finding resonance points. use U as short for u_j

> 

equr := u(vars) = eq19ub+eval(U(vars)*g(vars)^(r+alpha[1]), eqb); eqvr := v(vars) = eq19vb+eval(V(vars)*g(vars)^(r+alpha[2]), eqb)

u(x, t) = -4*(diff(g(x, t), x))^2/g(x, t)^2+U(x, t)*g(x, t)^(r-2)

v(x, t) = 2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2/g(x, t)^2+V(x, t)*g(x, t)^(r-2)

We substitute into the pdes, and seek "the coefficients of the dominant terms",

> 

eqr1 := expand(eval(pde1, {equr, eqvr})); eqr2 := expand(eval(pde2, {equr, eqvr}))

8*(diff(g(x, t), x))*(diff(diff(g(x, t), t), x))/g(x, t)^2-8*(diff(g(x, t), x))^2*(diff(g(x, t), t))/g(x, t)^3-15*a*U(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))*r/g(x, t)^4+3*a*U(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))*r^2/g(x, t)^4-4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r/g(x, t)^5-8*V(x, t)*g(x, t)^r*6^(1/2)*a^(1/2)*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4+a*U(x, t)*g(x, t)^r*(diff(diff(diff(g(x, t), x), x), x))*r/g(x, t)^3+3*a*(diff(diff(U(x, t), x), x))*g(x, t)^r*(diff(g(x, t), x))*r/g(x, t)^3+3*a*(diff(U(x, t), x))*g(x, t)^r*(diff(diff(g(x, t), x), x))*r/g(x, t)^3+3*a*(diff(g(x, t), x))^2*(diff(U(x, t), x))*g(x, t)^r*r^2/g(x, t)^4-15*a*(diff(g(x, t), x))^2*(diff(U(x, t), x))*g(x, t)^r*r/g(x, t)^4+a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r^3/g(x, t)^5-9*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r^2/g(x, t)^5+2*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r/g(x, t)^5-30*a*U(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4+6*a*U(x, t)^2*(g(x, t)^r)^2*(diff(g(x, t), x))*r/g(x, t)^5-4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2*(diff(V(x, t), x))*g(x, t)^r/g(x, t)^4+16*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r/g(x, t)^5+2*U(x, t)*g(x, t)^r*(diff(g(x, t), t))/g(x, t)^3-(diff(U(x, t), t))*g(x, t)^r/g(x, t)^2-24*a*(diff(diff(g(x, t), x), x))*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^2+96*a*(diff(diff(g(x, t), x), x))^2*(diff(g(x, t), x))/g(x, t)^3+56*a*(diff(g(x, t), x))^2*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3-8*a*(diff(g(x, t), x))*(diff(diff(diff(diff(g(x, t), x), x), x), x))/g(x, t)^2+a*(diff(diff(diff(U(x, t), x), x), x))*g(x, t)^r/g(x, t)^2-120*a*(diff(g(x, t), x))^3*(diff(diff(g(x, t), x), x))/g(x, t)^4-2*V(x, t)*(g(x, t)^r)^2*(diff(V(x, t), x))/g(x, t)^4+4*V(x, t)^2*(g(x, t)^r)^2*(diff(g(x, t), x))/g(x, t)^5-U(x, t)*g(x, t)^r*(diff(g(x, t), t))*r/g(x, t)^3-2*a*U(x, t)*g(x, t)^r*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3-6*a*(diff(diff(U(x, t), x), x))*g(x, t)^r*(diff(g(x, t), x))/g(x, t)^3-6*a*(diff(U(x, t), x))*g(x, t)^r*(diff(diff(g(x, t), x), x))/g(x, t)^3-6*a*(diff(g(x, t), x))^2*(diff(U(x, t), x))*g(x, t)^r/g(x, t)^4+72*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r/g(x, t)^5+6*a*U(x, t)*(g(x, t)^r)^2*(diff(U(x, t), x))/g(x, t)^4-12*a*U(x, t)^2*(g(x, t)^r)^2*(diff(g(x, t), x))/g(x, t)^5-2*V(x, t)^2*(g(x, t)^r)^2*(diff(g(x, t), x))*r/g(x, t)^5

12*U(x, t)*g(x, t)^r*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3/g(x, t)^5-3*U(x, t)*(g(x, t)^r)^2*V(x, t)*(diff(g(x, t), x))*r/g(x, t)^5-3*V(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))*r^2/g(x, t)^4+15*V(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))*r/g(x, t)^4+6*(diff(diff(V(x, t), x), x))*g(x, t)^r*(diff(g(x, t), x))/g(x, t)^3+6*(diff(V(x, t), x))*g(x, t)^r*(diff(diff(g(x, t), x), x))/g(x, t)^3+2*V(x, t)*g(x, t)^r*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3+2*V(x, t)*g(x, t)^r*(diff(g(x, t), t))/g(x, t)^3-60*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3*(diff(diff(g(x, t), x), x))/g(x, t)^4-12*U(x, t)*g(x, t)^r*6^(1/2)*a^(1/2)*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4-18*V(x, t)*g(x, t)^r*(diff(g(x, t), x))*(diff(diff(g(x, t), x), x))/g(x, t)^4-(diff(diff(diff(V(x, t), x), x), x))*g(x, t)^r/g(x, t)^2-(diff(V(x, t), t))*g(x, t)^r/g(x, t)^2-6*(diff(g(x, t), x))^2*(diff(V(x, t), x))*g(x, t)^r/g(x, t)^4-3*U(x, t)*(g(x, t)^r)^2*(diff(V(x, t), x))/g(x, t)^4-(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r^3/g(x, t)^5-V(x, t)*g(x, t)^r*(diff(g(x, t), t))*r/g(x, t)^3-3*(diff(diff(V(x, t), x), x))*g(x, t)^r*(diff(g(x, t), x))*r/g(x, t)^3-3*(diff(V(x, t), x))*g(x, t)^r*(diff(diff(g(x, t), x), x))*r/g(x, t)^3-V(x, t)*g(x, t)^r*(diff(diff(diff(g(x, t), x), x), x))*r/g(x, t)^3-3*(diff(g(x, t), x))^2*(diff(V(x, t), x))*g(x, t)^r*r^2/g(x, t)^4-14*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r/g(x, t)^5+6*U(x, t)*(g(x, t)^r)^2*V(x, t)*(diff(g(x, t), x))/g(x, t)^5+15*(diff(g(x, t), x))^2*(diff(V(x, t), x))*g(x, t)^r*r/g(x, t)^4+9*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r^2/g(x, t)^5-4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))*(diff(diff(g(x, t), t), x))/g(x, t)^2+4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2*(diff(g(x, t), t))/g(x, t)^3+48*6^(1/2)*a^(1/2)*(diff(diff(g(x, t), x), x))^2*(diff(g(x, t), x))/g(x, t)^3+28*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^3-4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))*(diff(diff(diff(diff(g(x, t), x), x), x), x))/g(x, t)^2-12*6^(1/2)*a^(1/2)*(diff(diff(g(x, t), x), x))*(diff(diff(diff(g(x, t), x), x), x))/g(x, t)^2

Want the dominant (least value).

> 

exps1 := map(findexp, [op(eqr1)]); dom1 := `assuming`([min(select(has, %, r))], [r::posint]); exps2 := map(findexp, [op(eqr2)]); dom2 := `assuming`([min(select(has, %, r))], [r::posint])

[-2, -3, r-4, r-4, r-5, r-4, r-3, r-3, r-3, r-4, r-4, r-5, r-5, r-5, r-4, 2*r-5, r-4, r-5, r-3, r-2, -2, -3, -3, -2, r-2, -4, -4+2*r, 2*r-5, r-3, r-3, r-3, r-3, r-4, r-5, -4+2*r, 2*r-5, 2*r-5]

r-5

[r-5, 2*r-5, r-4, r-4, r-3, r-3, r-3, r-3, -4, r-4, r-4, r-2, r-2, r-4, -4+2*r, r-5, r-3, r-3, r-3, r-3, r-4, r-5, 2*r-5, r-4, r-5, -2, -3, -3, -3, -2, -2]

r-5

So we want all terms with exponent dom; Find when the determinant is zero

> 

eqrr1 := select(proc (q) options operator, arrow; findexp(q) = dom1 end proc, eqr1); eqrr2 := select(proc (q) options operator, arrow; findexp(q) = dom2 end proc, eqr2); M, bb := GenerateMatrix([eqrr1, eqrr2], [U(vars), V(vars)])

-4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r/g(x, t)^5+a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r^3/g(x, t)^5-9*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r^2/g(x, t)^5+2*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r*r/g(x, t)^5+16*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r/g(x, t)^5+72*a*(diff(g(x, t), x))^3*U(x, t)*g(x, t)^r/g(x, t)^5

12*U(x, t)*g(x, t)^r*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3/g(x, t)^5-(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r^3/g(x, t)^5-14*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r/g(x, t)^5+9*(diff(g(x, t), x))^3*V(x, t)*g(x, t)^r*r^2/g(x, t)^5

Matrix(%id = 36893490698808109644), Vector[column](%id = 36893490698808109524)

Find resonant points, but <-1 doesn't count and Painleve doesn't apply; there is not a general solution.

> 

factor(Determinant(M)); [solve(%, r)]; rp := remove(`<`, %, -1); maxrp := max(rp)

-(diff(g(x, t), x))^6*(g(x, t)^r)^2*a*(r-3)*(r-4)*(r-6)*(r-8)*(r+2)*(r+1)/g(x, t)^10

[-2, -1, 3, 4, 6, 8]

[-1, 3, 4, 6, 8]

8

Step 3 -verifying the compatibility condition

> 

equ := u(vars) = eq19ub+eval(add(U[i](vars)*g(vars)^(i+alpha[1]), i = 1 .. maxrp), eqb); eqv := v(vars) = eq19vb+eval(add(V[i](vars)*g(vars)^(i+alpha[2]), i = 1 .. maxrp), eqb)

u(x, t) = -4*(diff(g(x, t), x))^2/g(x, t)^2+U[1](x, t)/g(x, t)+U[2](x, t)+U[3](x, t)*g(x, t)+U[4](x, t)*g(x, t)^2+U[5](x, t)*g(x, t)^3+U[6](x, t)*g(x, t)^4+U[7](x, t)*g(x, t)^5+U[8](x, t)*g(x, t)^6

v(x, t) = 2*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2/g(x, t)^2+V[1](x, t)/g(x, t)+V[2](x, t)+V[3](x, t)*g(x, t)+V[4](x, t)*g(x, t)^2+V[5](x, t)*g(x, t)^3+V[6](x, t)*g(x, t)^4+V[7](x, t)*g(x, t)^5+V[8](x, t)*g(x, t)^6

Substitute back into the pde, and collect wrt g(vars) = G.

> 

coll1 := collect(subs(g[] = G, ToJet(eval(pde1, {equ, eqv}), [g(vars), seq(U[i](vars), i = 1 .. maxrp), seq(V[i](vars), i = 1 .. maxrp)], jetnotation = jetvariableswithbrackets)), G); coll2 := collect(subs(g[] = G, ToJet(eval(pde2, {equ, eqv}), [g(vars), seq(U[i](vars), i = 1 .. maxrp), seq(V[i](vars), i = 1 .. maxrp)], jetnotation = jetvariableswithbrackets)), G); `~`[ldegree]([coll1, coll2], G); mindegree1, mindegree2 := %[]

[-4, -4]

-4, -4

Find the coefficients of the maxrp least-degree terms, since we have maxrp variables to solve for

> 

eqs1 := FromJet([seq(coeff(coll1, G, mindegree1+i), i = 0 .. maxrp-1)], [g(vars), seq(U[i](vars), i = 1 .. maxrp), seq(V[i](vars), i = 1 .. maxrp)]); eqs2 := FromJet([seq(coeff(coll2, G, mindegree2+i), i = 0 .. maxrp-1)], [g(vars), seq(U[i](vars), i = 1 .. maxrp), seq(V[i](vars), i = 1 .. maxrp)])

 

Solving the first equation gives {u[1],v[1]} - Eq 3.24

> 

uv1 := solve({eqs1[1], eqs2[1]}, {U[1](vars), V[1](vars)})

{U[1](x, t) = 4*(diff(diff(g(x, t), x), x)), V[1](x, t) = -2*a^(1/2)*6^(1/2)*(diff(diff(g(x, t), x), x))}

Substitute this in the second equation and solve for U[2] - Eq 3.25

> 

dsubs(uv1, {eqs1[2], eqs2[2]}); uv2 := solve(%, {U[2](vars), V[2](vars)})

{-12*6^(1/2)*a^(1/2)*(diff(diff(g(x, t), x), x))^2*(diff(g(x, t), x))+16*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2*(diff(diff(diff(g(x, t), x), x), x))+12*U[2](x, t)*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3+4*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^2*(diff(g(x, t), t)), 8*V[2](x, t)*6^(1/2)*a^(1/2)*(diff(g(x, t), x))^3+48*(diff(g(x, t), x))^3*U[2](x, t)*a+32*a*(diff(g(x, t), x))^2*(diff(diff(diff(g(x, t), x), x), x))-24*a*(diff(diff(g(x, t), x), x))^2*(diff(g(x, t), x))-8*(diff(g(x, t), x))^2*(diff(g(x, t), t))}

{U[2](x, t) = -(1/3)*(4*(diff(g(x, t), x))*(diff(diff(diff(g(x, t), x), x), x))+(diff(g(x, t), x))*(diff(g(x, t), t))-3*(diff(diff(g(x, t), x), x))^2)/(diff(g(x, t), x))^2, V[2](x, t) = (1/6)*(4*(diff(diff(diff(g(x, t), x), x), x))*(diff(g(x, t), x))*a+2*(diff(g(x, t), t))*(diff(g(x, t), x))*a-3*(diff(diff(g(x, t), x), x))^2*a+(diff(g(x, t), x))*(diff(g(x, t), t)))*6^(1/2)/((diff(g(x, t), x))^2*a^(1/2))}

This one is a resonance

> 

eqs3 := simplify(dsubs(uv1, uv2, {eqs1[3], eqs2[3]})); M3, bb3 := GenerateMatrix(eqs3, [U[3](vars), V[3](vars)])

Matrix(%id = 36893490698808102772), Vector[column](%id = 36893490698808102652)

The determinant = 0 expected for a resonance. Rank of matrix and augmented matrix are the same - if we solve for U[3] and V[3], we find V[3] is arbitrary and U[3] also has a arbitrary function

> 

Determinant(M3); Rank(M3); Rank(`<|>`(M3, bb3)); LinearSolve(M3, bb3)

0

1

1

Vector[column](%id = 36893490698808773556)

This one is also a resonance

> 

eqs4 := simplify(dsubs(uv1, uv2, U[3](vars) = 0, V[3](vars) = 0, {eqs1[4], eqs2[4]})); M4, bb4 := GenerateMatrix(eqs4, [U[4](vars), V[4](vars)])

Matrix(%id = 36893490698894126244), Vector[column](%id = 36893490698894126004)

Here matrix and augmented matrix have different ranks, no solution.

> 

Determinant(M4); Rank(M4); Rank(`<|>`(M4, bb4)); LinearSolve(M4, bb4)

0

1

2

Error, (in LinearAlgebra:-LinearSolve) inconsistent system

We find a solution for U[5] and V[5]

> 

eqs5 := simplify(dsubs(uv1, uv2, U[3](vars) = 0, V[3](vars) = 0, U[4](vars) = 0, V[4](vars) = 0, {eqs1[5], eqs2[5]})); M5, bb5 := GenerateMatrix(eqs5, [U[5](vars), V[5](vars)])

Matrix(%id = 36893490698796769636), Vector[column](%id = 36893490698796769516)

> 

Determinant(M5)

-252*a*(diff(g(x, t), x))^6

> 

uv5 := {seq(Equate(`<,>`(U[5](vars), V[5](vars)), LinearSolve(M5, bb5)))}

This one is a resonance

> 

eqs6 := simplify(dsubs(uv1, uv2, uv5, U[3](vars) = 0, V[3](vars) = 0, U[4](vars) = 0, V[4](vars) = 0, {eqs1[6], eqs2[6]})); M6, bb6 := GenerateMatrix(eqs6, [U[6](vars), V[6](vars)])

Matrix(%id = 36893490698867434788), Vector[column](%id = 36893490698867434668)

Determinant = 0

> 

Determinant(M6); LinearSolve(M6, bb6)

0

Error, (in LinearAlgebra:-LinearSolve) inconsistent system

Solution for U[7], V[7]

> 

eqs7 := simplify(dsubs(uv1, uv2, uv5, U[3](vars) = 0, V[3](vars) = 0, U[4](vars) = 0, V[4](vars) = 0, U[6](vars) = 0, V[6](vars) = 0, {eqs1[7], eqs2[7]})); M7, bb7 := GenerateMatrix(eqs7, [U[7](vars), V[7](vars)])

Matrix(%id = 36893490698867056756), Vector[column](%id = 36893490698867056636)

> 

Determinant(M7); LinearSolve(M7, bb7)

-864*a*(diff(g(x, t), x))^6

Vector[column](%id = 36893490698894124924)

And this is the last resonance

> 

eqs8 := simplify(dsubs(uv1, uv2, uv5, U[3](vars) = 0, V[3](vars) = 0, U[4](vars) = 0, V[4](vars) = 0, U[6](vars) = 0, V[6](vars) = 0, {eqs1[8], eqs2[8]})); M8, bb8 := GenerateMatrix(eqs8, [U[8](vars), V[8](vars)])

Matrix(%id = 36893490698740702252), Vector[column](%id = 36893490698740702132)

> 

Determinant(M8)

0

I don't understand the bit about needing a=1/2 for compatibility

> 

NULL

Download AllStepscKdV.mw

cursor jumping...

dharr 8355
August 26 2025
1 0

Yes this factorial in the denominator jumping drives me crazy. It seems to be a side effect of the new 2-D entry system that abbreviates entering x^2 cursor-right + 2 to  x^2+2. It is still present in Maple 2025.1, so there is currently not a solution.

conversion to trig...

dharr 8355
August 25 2025
2 2

Well, this is probably cheating since I saw the answer, but maybe I would have thought of it anyway.

> 

restart
> 

interface(version)

`Standard Worksheet Interface, Maple 2025.1, Windows 11, June 12 2025 Build ID 1932578`

> 

term := 2*cos(5*arcsin((1/2)*x))

2*cos(5*arcsin((1/2)*x))

> 

term1 := expand(term)

-3*(-x^2+4)^(1/2)*x^2+(-x^2+4)^(1/2)+(-x^2+4)^(1/2)*x^4

Well, combine is usually the opposite of expand, so the first thing I'd try is

> 

combine(term1)

-3*(-x^2+4)^(1/2)*x^2+(-x^2+4)^(1/2)+(-x^2+4)^(1/2)*x^4

I tried a lot of converts without success. Of course not many algebraic things are likely to be expressible as trig functions, so it is just "luck" that term1 came out so simply.

So I'll give Maple a hand, trying to forget that I saw the answer already. Now sqrt(1-x*2) for x = sin/cos(u) gives cos/sin(u), and for the 4's being 2^2 suggests a substitution

> 

tr := x = 2*sin(u); p := simplify(eval(term1, tr))

x = 2*sin(u)

(32*cos(u)^5-40*cos(u)^3+10*cos(u))*csgn(cos(u))

We could play around with assumptions to handle the csgn, but let's be reckless

> 

p2 := simplify(p, symbolic)

2*cos(5*u)

To go back to x we need to inverse transformation

> 

itr := isolate(tr, u)

u = arcsin((1/2)*x)

which gives what we want.

> 

f := eval(p2, itr)

2*cos(5*arcsin((1/2)*x))

If we had used cos instead

> 

tr2 := x = 2*cos(u); q := simplify(eval(term1, tr2), symbolic)

x = 2*cos(u)

2*sin(5*u)

> 

itr2 := isolate(tr2, u); g := eval(q, itr2)

u = arccos((1/2)*x)

2*sin(5*arccos((1/2)*x))

And the plots are on top of each other

> 

plot([f, g], x = -2.5 .. 2.5)

which we can verify exactly

> 

simplify(f-g)

0

> 

NULL

Download test2.mw

extraneous solutions generated...

dharr 8355
August 24 2025
2 0

Of course it is a bug, but notice that solve gives the same result. It's the classic problem of introducing extraneous solutions by squaring both sides of an equation.

> 

restart;
> 

sys := diff~(sqrt(2*a^2 - 8*a + 10) + sqrt(b^2 - 6*b + 10) + sqrt(2*a^2 - 2*a*b + b^2), {a, b});

{(1/2)*(4*a-8)/(2*a^2-8*a+10)^(1/2)+(1/2)*(4*a-2*b)/(2*a^2-2*a*b+b^2)^(1/2), (1/2)*(2*b-6)/(b^2-6*b+10)^(1/2)+(1/2)*(-2*a+2*b)/(2*a^2-2*a*b+b^2)^(1/2)}

> 

solve(sys);

{a = 5/3, b = 5/2}, {a = a, b = 2*a/(a-1)}

If you solve with infolevel[solve]:=5; you see it solved 5 equations, so plenty of opportunity for extraneous solutions.
So I'm guessing it did something like this:

> 

neweqs := {2*a^2-8*a+10 = s1^2, 2*a^2 - 2*a*b + b^2 = s2^2, b^2 - 6*b + 10 = s3^2};
sys2a:=simplify(eval(sys, neweqs)) assuming s1>0,s2>0,s3>0;
sys2:=numer~(%) union (lhs-rhs)~(neweqs);

{2*a^2-8*a+10 = s1^2, b^2-6*b+10 = s3^2, 2*a^2-2*a*b+b^2 = s2^2}

{(2*(a-2)*s2+(2*a-b)*s1)/(s1*s2), ((b-3)*s2-(a-b)*s3)/(s3*s2)}

{2*a^2-s1^2-8*a+10, b^2-s3^2-6*b+10, 2*a^2-2*a*b+b^2-s2^2, 2*a*s1+2*a*s2-b*s1-4*s2, -a*s3+b*s2+b*s3-3*s2}

The last two answers here are the problematic ones, but actually the first six are special cases of a=a, b=2*a/(a-1). The problem is that s1, s2, s2 are not all positive, so the squaring in neweqs with the assumptions led to the problem

> 

ans:=sort([solve(sys2, {a, b, s1, s2, s3},explicit)]);

[{a = 2, b = 4, s1 = 2^(1/2), s2 = 2*2^(1/2), s3 = -2^(1/2)}, {a = 2, b = 4, s1 = -2^(1/2), s2 = -2*2^(1/2), s3 = 2^(1/2)}, {a = 3, b = 3, s1 = -2, s2 = 3, s3 = -1}, {a = 3, b = 3, s1 = -2, s2 = 3, s3 = 1}, {a = 3, b = 3, s1 = 2, s2 = -3, s3 = -1}, {a = 3, b = 3, s1 = 2, s2 = -3, s3 = 1}, {a = 5/3, b = 5/2, s1 = -(2/3)*5^(1/2), s2 = -(5/6)*5^(1/2), s3 = -(1/2)*5^(1/2)}, {a = 5/3, b = 5/2, s1 = (2/3)*5^(1/2), s2 = (5/6)*5^(1/2), s3 = (1/2)*5^(1/2)}, {a = a, b = 2*a/(a-1), s1 = (2*a^2-8*a+10)^(1/2), s2 = -(2*a^2-8*a+10)^(1/2)*a/(a-1), s3 = (2*a^2-8*a+10)^(1/2)/(a-1)}, {a = a, b = 2*a/(a-1), s1 = -(2*a^2-8*a+10)^(1/2), s2 = (2*a^2-8*a+10)^(1/2)*a/(a-1), s3 = -(2*a^2-8*a+10)^(1/2)/(a-1)}]

All these solutions check out for the expanded system

> 

seq(simplify(eval(sys2,ans[i])), i=1..nops(ans));

{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}

Taking numer wasn't a problem

> 

seq(simplify(eval(sys2a,ans[i])), i=1..nops(ans));

{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}

OK in neweqs

> 

seq(simplify(eval(neweqs,ans[i])), i=1..nops(ans));

{2 = 2, 8 = 8}, {2 = 2, 8 = 8}, {1 = 1, 4 = 4, 9 = 9}, {1 = 1, 4 = 4, 9 = 9}, {1 = 1, 4 = 4, 9 = 9}, {1 = 1, 4 = 4, 9 = 9}, {5/4 = 5/4, 20/9 = 20/9, 125/36 = 125/36}, {5/4 = 5/4, 20/9 = 20/9, 125/36 = 125/36}, {(2*a^2-8*a+10)/(a-1)^2 = (2*a^2-8*a+10)/(a-1)^2, 2*a^2*(a^2-4*a+5)/(a-1)^2 = 2*a^2*(a^2-4*a+5)/(a-1)^2, 2*a^2-8*a+10 = 2*a^2-8*a+10}, {(2*a^2-8*a+10)/(a-1)^2 = (2*a^2-8*a+10)/(a-1)^2, 2*a^2*(a^2-4*a+5)/(a-1)^2 = 2*a^2*(a^2-4*a+5)/(a-1)^2, 2*a^2-8*a+10 = 2*a^2-8*a+10}

But not in sys

> 

seq(simplify(eval(sys,ans[i])), i=1..nops(ans));

{0, 2^(1/2)}, {0, 2^(1/2)}, {0, 2}, {0, 2}, {0, 2}, {0, 2}, {0}, {0}, {(a-2)*(2/(2*a^2-8*a+10)^(1/2)+a*2^(1/2)/((a^2*(a^2-4*a+5)/(a-1)^2)^(1/2)*(a-1))), -(1/2)*2^(1/2)*(a-3)*(1/((a^2-4*a+5)/(a-1)^2)^(1/2)+a/(a^2*(a^2-4*a+5)/(a-1)^2)^(1/2))/(a-1)}, {(a-2)*(2/(2*a^2-8*a+10)^(1/2)+a*2^(1/2)/((a^2*(a^2-4*a+5)/(a-1)^2)^(1/2)*(a-1))), -(1/2)*2^(1/2)*(a-3)*(1/((a^2-4*a+5)/(a-1)^2)^(1/2)+a/(a^2*(a^2-4*a+5)/(a-1)^2)^(1/2))/(a-1)}

 

NULL

Download realsolve.mw

Steps 1 and 2 and 3...

dharr 8355
August 23 2025
4 15

[Edit - modified so step 1 works as described in the literature - see more examples with refs below.]

> 

restart
> 

with(PDEtools)
> 

undeclare(prime, quiet)
> 

vars := x, t; declare(u(vars), quiet); declare(Phi(vars), quiet); declare(phi(vars), quiet)

Burgers eqn  - see https://doi.org/10.1063/1.525721

> 

pde1 := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x))-sigma*(diff(u(x, t), x, x))

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x))-sigma*(diff(diff(u(x, t), x), x))

Routine to find symbolic exponent of phi

> 

findexp := proc (T) options operator, arrow; simplify(phi(vars)*Physics:-diff(T, phi(vars))/T) end proc

Step 1

> 

equ1 := u(vars) = chi*phi(vars)^a

u(x, t) = chi*phi(x, t)^a

> 

eqphi := expand(eval(pde1, equ1))

chi*phi(x, t)^a*a*(diff(phi(x, t), t))/phi(x, t)+chi^2*(phi(x, t)^a)^2*a*(diff(phi(x, t), x))/phi(x, t)-sigma*chi*phi(x, t)^a*a^2*(diff(phi(x, t), x))^2/phi(x, t)^2-sigma*chi*phi(x, t)^a*a*(diff(diff(phi(x, t), x), x))/phi(x, t)+sigma*chi*phi(x, t)^a*a*(diff(phi(x, t), x))^2/phi(x, t)^2

Find all the distinct powers of phi

> 

exps := [(`minus`(map(findexp, {op(eqphi+forceplus)}), {0}))[]]

[a-2, a-1, 2*a-1]

We are supposed to compare all pairs of lowest exponents. So first find the lowest line for each slope.

> 

convert(ListTools:-Classify(diff, exps, a), list); expmins := `~`[min](%)

[{a-2, a-1}, {2*a-1}]

[a-2, 2*a-1]

Now compare all pairs of these to find the intersection with least a

> 

eqa := a = min(seq(seq(solve(expmins[i] = expmins[j], a), i = 1 .. j-1), j = 1 .. nops(expmins)))

a = -1

> 

eq6 := u(vars) = eval(chi*phi(vars)^a, eqa)

u(x, t) = chi/phi(x, t)

The following has all the terms, not just those from the highest derivative and higher nonlinearity

> 

eq7 := expand(eval(pde1, eq6))

-chi*(diff(phi(x, t), t))/phi(x, t)^2-chi^2*(diff(phi(x, t), x))/phi(x, t)^3-2*sigma*chi*(diff(phi(x, t), x))^2/phi(x, t)^3+sigma*chi*(diff(diff(phi(x, t), x), x))/phi(x, t)^2

We want only the terms with phi^(mindeg), but we cannot collect wrt phi, so convert it to a simple name before collecting

> 

eq7c := subs(phi[] = F, ToJet(eq7, phi(vars), jetnotation = jetvariableswithbrackets))

-chi*phi[t]/F^2-chi^2*phi[x]/F^3-2*sigma*chi*phi[x]^2/F^3+sigma*chi*phi[x, x]/F^2

> 

mindeg := ldegree(eq7c, F)

-3

> 

S := FromJet(coeff(collect(eq7c, F), F, mindeg), phi(vars))

-2*sigma*chi*(diff(phi(x, t), x))^2-chi^2*(diff(phi(x, t), x))

> 

eqchi := chi = solve(S, chi)[2]

chi = -2*(diff(phi(x, t), x))*sigma

Following is u[0]/phi^a, which we use in step 2

> 

u0phia := rhs(eval(eq6, eqchi))

-2*(diff(phi(x, t), x))*sigma/phi(x, t)

So we found (chi) and (a) which was our target for first step

Step 2 - finding resonance points. use Phi as short for phi_j

> 

eq11 := u(vars) = u0phia+eval(Phi(vars)*phi(vars)^(j+a), eqa)

u(x, t) = -2*(diff(phi(x, t), x))*sigma/phi(x, t)+Phi(x, t)*phi(x, t)^(j-1)

We substitute into the pde, and seek "the coefficients of the dominant terms",

> 

eq11b := expand(eval(pde1, eq11))

Following routine finds the exponent of phi in a term T

> 

findexp := proc (T) options operator, arrow; simplify(phi(vars)*Physics:-diff(T, phi(vars))/T) end proc

proc (T) options operator, arrow; simplify(phi(vars)*Physics:-diff(T, phi(vars))/T) end proc

Want the dominant (least value).

> 

exps := map(findexp, [op(eq11b)]); dom := `assuming`([min(select(has, %, j))], [j::posint])

[-1, j-1, -1, j-2, j-3, j-2, j-3, j-2, -2, 2*j-2, -3+2*j, j-1, -2, j-2, j-3, j-2, -3+2*j]

j-3

So we want all terms with exponent dom; Solve the equation with these terms to find the resonant points.

> 

eq12 := select(proc (q) options operator, arrow; findexp(q) = dom end proc, eq11b); rp := [solve(eq12, j)]; maxrp := max(rp)

(diff(phi(x, t), x))^2*sigma*Phi(x, t)*phi(x, t)^j*j/phi(x, t)^3-sigma*Phi(x, t)*phi(x, t)^j*(diff(phi(x, t), x))^2*j^2/phi(x, t)^3+2*(diff(phi(x, t), x))^2*sigma*Phi(x, t)*phi(x, t)^j/phi(x, t)^3

[2, -1]

2

Step 3 -verifying the compatibility condition

> 

eq14 := u(vars) = u0phia+eval(add(Phi[j](vars)*phi(vars)^(j+a), j = 1 .. maxrp), eqa)

u(x, t) = -2*(diff(phi(x, t), x))*sigma/phi(x, t)+Phi[1](x, t)+Phi[2](x, t)*phi(x, t)

Substitute back into the pde, and collect wrt phi(vars) = F.

> 

depvars := [phi(vars), seq(Phi[i](vars), i = 1 .. maxrp)]; coll := collect(subs(phi[] = F, ToJet(eval(pde1, eq14), depvars, jetnotation = jetvariableswithbrackets)), F); mindegree := ldegree(coll, F)

Phi[2][]*Phi[2][x]*F^2+(Phi[2][t]+Phi[1][]*Phi[2][x]+Phi[2][]*(phi[x]*Phi[2][]+Phi[1][x])-sigma*Phi[2][x, x])*F+Phi[1][t]+Phi[2][]*phi[t]-2*phi[x]*sigma*Phi[2][x]+Phi[1][]*(phi[x]*Phi[2][]+Phi[1][x])-2*Phi[2][]*phi[x, x]*sigma-sigma*(2*phi[x]*Phi[2][x]+phi[x, x]*Phi[2][]+Phi[1][x, x])+(-2*phi[x, t]*sigma-2*phi[x]*sigma*(phi[x]*Phi[2][]+Phi[1][x])-2*Phi[1][]*phi[x, x]*sigma+2*Phi[2][]*phi[x]^2*sigma+2*sigma^2*phi[x, x, x])/F+(-2*sigma^2*phi[x]*phi[x, x]+2*sigma*phi[x]^2*Phi[1][]+2*sigma*phi[t]*phi[x])/F^2

-2

Find the coefficients of the maxrp least-degree terms, since we have maxrp variables to solve for

> 

eqs := FromJet([seq(coeff(coll, F, mindegree+i), i = 0 .. maxrp-1)], depvars)

[-2*(diff(phi(x, t), x))*sigma^2*(diff(diff(phi(x, t), x), x))+2*Phi[1](x, t)*(diff(phi(x, t), x))^2*sigma+2*(diff(phi(x, t), x))*sigma*(diff(phi(x, t), t)), -2*(diff(diff(phi(x, t), t), x))*sigma-2*(diff(phi(x, t), x))*sigma*(Phi[2](x, t)*(diff(phi(x, t), x))+diff(Phi[1](x, t), x))-2*Phi[1](x, t)*(diff(diff(phi(x, t), x), x))*sigma+2*Phi[2](x, t)*(diff(phi(x, t), x))^2*sigma+2*sigma^2*(diff(diff(diff(phi(x, t), x), x), x))]

For the Burger's equation, solving the first equation gives u[1] - Eq. 3.7

> 

u1 := Phi[1](vars) = solve(eqs[1], Phi[1](vars))

Phi[1](x, t) = ((diff(diff(phi(x, t), x), x))*sigma-(diff(phi(x, t), t)))/(diff(phi(x, t), x))

If we substitute this in the second equation, it simplifies to zero, which is why we cannot just solve the system.
This is because j=2 is a resonance point.
We have found Eq. 3.8; the factor in the first line below is the derivative wrt x of the eqn for u1.

> 

simplify(eqs[2]); simplify(dsubs(u1, eqs[2]))

-2*sigma*((diff(phi(x, t), x))*(diff(Phi[1](x, t), x))+Phi[1](x, t)*(diff(diff(phi(x, t), x), x))-sigma*(diff(diff(diff(phi(x, t), x), x), x))+diff(diff(phi(x, t), t), x))

0

> 

NULL

Download AllstepsBurgerNew.mw

Example 2: OP's Hirota Bilinear eqn

AllstepsBurgerNew.mw

Example 3: KdV equation

AllStepsKdVnew2.mw

 

shape = symmetric...

dharr 8355
August 13 2025
3 15

If you create the Matrix with option shape = symmetric, then you will get only real eigenvalues and eigenvectors, and Maple will use a more efficient method specialized for symmetric matrices.

(If you just want to get rid of the +0*I after the fact, you can use simplify(..., zero)).

vector graphics...

dharr 8355
August 06 2025
1 0

You didn't say how you want to export your plot, but if I assume you want vector graphics then there are two possibilities, encapsulated postscript (EPS) or SVG. (possibly also PDF but I've never tried it)

If I make a plot with no background, say plot(sin(x),axes=none): and export it as EPS with the right-click menu, then look at it with CorelDraw I find it has a white background (which I usually delete and then process further), so this does not work for you.

However, exporting as SVG gives a file that does not have a background, just the curve, so this is what you want.

If you actually want to export as a bitmap image (with loss of resolution), then I don't have much experience with that. Which formats would you want then?

Improvements...

dharr 8355
August 04 2025
2 1

An extract from the ?updates,Maple2017,Performance help page: "Maple 2017 includes a new compiled C implementation of the FGLM algorithm for computing a lexicographic Groebner basis from a total degree basis when there are a finite number of solutions.  This routine is used automatically by Groebner:-Basis when applicable and by the solve command when solving polynomial systems.  The example below runs about 200 times faster in Maple 2017."

The ?updates,Maple2018,index help page and the current help page for ?Groebner,Basis say "The Groebner[Basis] command was updated in Maple 2018", which I assume means the command itself has changed but not necessarily the algorithms.

(If you are using the bases for solving systems of polynomial equations, then there are new methods that don't use Groebner bases that may be more efficient.)

ctrl-shift J/K...

dharr 8355
July 31 2025
2 1

If I understand correctly what you mean, use ctrl-shift-K to insert above and ctrl-shift-J to insert below. (command-shift-K and command-shift-J on Mac). See ?worksheet,reference,hotmac or ..hotwin or .. hotunix.

setup(computederivatives = false);...

dharr 8355
July 31 2025
3 0

The default behavior was changed in Maple 2020, see ?updates,Maple2020,IntegralTransforms.

To get the behavior you want, use setup(computederivatives = false);

Download laplace.mw

contourplot...

dharr 8355
July 30 2025
1 0

This is an answer to the followup question, now deleted, about plotting the solution in cartesian coordinates. Here is an example:

PDE_upd_5.mw

Note that VectorCalculus:-Laplacian(u(xi, eta), 'bipolar'[xi, eta]) assumes a=1, so I substituted your earlier laplacian.

parametric plot...

dharr 8355
July 30 2025
0 0

Since you have already the functions x(xi,eta) and y(xi,eta), you can use the parametric version of contourplot.

PDE_upd_5.mw

I would probably play around with custom contour values because your large spike dominates and so the contours near zero do not have good resolution, but the I think you cannot then use the colorbar.

NOTE: I replaced VectorCalculus:-Laplacian(u(xi, eta), 'bipolar'[xi, eta]) with your earlier hand-calculated Laplacian, since the built-in one uses a = 1, which you had before but is not applicable here. VectorCalculus is supposed to have a mechanism for changing the a parameter (SetCoordinateParameters) but I couldn't get it to work with Laplacian. 

extra x...

dharr 8355
July 27 2025
1 3

Yes, your substitutions missed an x. You can capture this with 

eqs := {coeffs(collect(numer(normal(eq)), {X, Y, Z, x}, 'distributed'), {X, Y, Z, x})}

but I think that you cannot constrain tanh(k*x + p*y - t*w) together as one variable, since x, y, and t are independent. If the ln was not there you could expand(tanh(k*x + p*y - t*w)) and substitute the individual sinh, cosh etc, or just convert to exp, but in both cases the ln is a problem.

Verified...

dharr 8355
July 26 2025
1 1

I would almost always use simplify for this sort of thing, not is.

> 

restart
> 

d := proc (i, j) options operator, arrow; (x[i]-x[j])^2+(y[i]-y[j])^2 end proc

proc (i, j) options operator, arrow; (x[i]-x[j])^2+(y[i]-y[j])^2 end proc

> 

q1:=d(1,2)*d(1,2)*d(3,4)+d(1,3)*d(1,3)*d(2,4)+d(1,4)*d(1,4)*d(2,3)+d(2,3)*d(2,3)*d(1,4)
+d(2,4)*d(2,4)*d(1,3)+d(3,4)*d(3,4)*d(1,2)+d(1,2)*d(2,3)*d(3,1)+d(1,2)*d(2,4)*d(4,1)
+d(1,3)*d(3,4)*d(4,1)+d(2,3)*d(3,4)*d(4,2);

((x[1]-x[2])^2+(y[1]-y[2])^2)^2*((x[3]-x[4])^2+(y[3]-y[4])^2)+((x[1]-x[3])^2+(y[1]-y[3])^2)^2*((x[2]-x[4])^2+(y[2]-y[4])^2)+((x[1]-x[4])^2+(y[1]-y[4])^2)^2*((x[2]-x[3])^2+(y[2]-y[3])^2)+((x[2]-x[3])^2+(y[2]-y[3])^2)^2*((x[1]-x[4])^2+(y[1]-y[4])^2)+((x[2]-x[4])^2+(y[2]-y[4])^2)^2*((x[1]-x[3])^2+(y[1]-y[3])^2)+((x[3]-x[4])^2+(y[3]-y[4])^2)^2*((x[1]-x[2])^2+(y[1]-y[2])^2)+((x[1]-x[2])^2+(y[1]-y[2])^2)*((x[2]-x[3])^2+(y[2]-y[3])^2)*((x[3]-x[1])^2+(y[3]-y[1])^2)+((x[1]-x[2])^2+(y[1]-y[2])^2)*((x[2]-x[4])^2+(y[2]-y[4])^2)*((x[4]-x[1])^2+(y[4]-y[1])^2)+((x[1]-x[3])^2+(y[1]-y[3])^2)*((x[3]-x[4])^2+(y[3]-y[4])^2)*((x[4]-x[1])^2+(y[4]-y[1])^2)+((x[2]-x[3])^2+(y[2]-y[3])^2)*((x[3]-x[4])^2+(y[3]-y[4])^2)*((x[4]-x[2])^2+(y[4]-y[2])^2)

> 

q2:=d(1,2)*d(2,3)*d(3,4)+d(1,3)*d(3,2)*d(2,4)+d(1,2)*d(2,4)*d(4,3)+d(1,4)*d(4,2)*d(2,3)
+d(1,3)*d(3,4)*d(4,2)+d(1,4)*d(4,3)*d(3,2)+d(2,3)*d(3,1)*d(1,4)+d(2,1)*d(1,3)*d(3,4)
+d(2,4)*d(4,1)*d(1,3)+d(2,1)*d(1,4)*d(4,3)+d(3,1)*d(1,2)*d(2,4)+d(3,2)*d(2,1)*d(1,4);

((x[1]-x[2])^2+(y[1]-y[2])^2)*((x[2]-x[3])^2+(y[2]-y[3])^2)*((x[3]-x[4])^2+(y[3]-y[4])^2)+((x[1]-x[3])^2+(y[1]-y[3])^2)*((x[3]-x[2])^2+(y[3]-y[2])^2)*((x[2]-x[4])^2+(y[2]-y[4])^2)+((x[1]-x[2])^2+(y[1]-y[2])^2)*((x[2]-x[4])^2+(y[2]-y[4])^2)*((x[4]-x[3])^2+(y[4]-y[3])^2)+((x[1]-x[4])^2+(y[1]-y[4])^2)*((x[4]-x[2])^2+(y[4]-y[2])^2)*((x[2]-x[3])^2+(y[2]-y[3])^2)+((x[1]-x[3])^2+(y[1]-y[3])^2)*((x[3]-x[4])^2+(y[3]-y[4])^2)*((x[4]-x[2])^2+(y[4]-y[2])^2)+((x[1]-x[4])^2+(y[1]-y[4])^2)*((x[4]-x[3])^2+(y[4]-y[3])^2)*((x[3]-x[2])^2+(y[3]-y[2])^2)+((x[2]-x[3])^2+(y[2]-y[3])^2)*((x[3]-x[1])^2+(y[3]-y[1])^2)*((x[1]-x[4])^2+(y[1]-y[4])^2)+((x[2]-x[1])^2+(y[2]-y[1])^2)*((x[1]-x[3])^2+(y[1]-y[3])^2)*((x[3]-x[4])^2+(y[3]-y[4])^2)+((x[2]-x[4])^2+(y[2]-y[4])^2)*((x[4]-x[1])^2+(y[4]-y[1])^2)*((x[1]-x[3])^2+(y[1]-y[3])^2)+((x[2]-x[1])^2+(y[2]-y[1])^2)*((x[1]-x[4])^2+(y[1]-y[4])^2)*((x[4]-x[3])^2+(y[4]-y[3])^2)+((x[3]-x[1])^2+(y[3]-y[1])^2)*((x[1]-x[2])^2+(y[1]-y[2])^2)*((x[2]-x[4])^2+(y[2]-y[4])^2)+((x[3]-x[2])^2+(y[3]-y[2])^2)*((x[2]-x[1])^2+(y[2]-y[1])^2)*((x[1]-x[4])^2+(y[1]-y[4])^2)

> 

simplify(q1-q2)

0

> 

``

Download distances.mw

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