dharr

Dr. David Harrington

9202 Reputation

22 Badges

21 years, 300 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 replies submitted by dharr

@salim-barzani The case here in the paper with e[1]=0 does not need the full analysis since it reduces to a quadratic in w^2, so it is simple to go through all the cases. I used here e[i] but in the worksheets with the pdes I used e__i.

restart

For a complex-conjugate pair of roots w1+`&+-`(I*w2) we can write this as the form on the rhs, which is used in the table

(w-w1-I*w2)*(w-w1+I*w2) = (w-w1)^2+w2^2; simplify((rhs-lhs)(%))

(w-w1-I*w2)*(w-w1+I*w2) = (w-w1)^2+w2^2

0

So the cases in the table are

(I) a complex-conjugate pair that occurs twice

(II) 1 single real root of multiplicity 4 (the root must be zero)

(III) a real root of multiplicity 2 and a different real root of multiplicity 2

(IV) a real root of multiplicity 3 and a real root of multiplicity 1

(V) a real root of multiplicity 2 and a complex-conjugate pair

(VI) four real roots all different

(VII) two different real roots and a complex conjugate pair

(VIII) two different complex conjugate pairs

(IX) a real root of multiplicity 2 and two other different real roots.

p := w^4+w^2*e[2]+w*e[1]+e[0]

w^4+w^2*e[2]+w*e[1]+e[0]

The discriminant is zero iff multiplicity > 1. This is four times theD__4 in the paper

d := discrim(p, w)

16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3

But for e[1] = 0, our polynomial and the discriminant simplifies to

p1 := eval(p, e[1] = 0); factor(eval(d, e[1] = 0))

w^4+w^2*e[2]+e[0]

16*e[0]*(-e[2]^2+4*e[0])^2

So if e[0]<>0 and e[0] <> e[2]^2/4 then all roots are distinct. This is cases (VI), (VII) and (VIII). As in the table D4<>0 in these cases and D4 = 0 for all the other cases.
If we have 4*e[0] = e[2]^2 then we have multiple roots. But in this case the polynomial is a perfect square and we have cases (I), (II), or (III). If e[2] is positive then the roots are imaginary and we have case (I) with w1=0. So anywhere in the paper for case (I) that w1 appears it should be set equal to zero. If e[2] is positive then we have case (III). If e[2] = 0 (and therefore e[0] = 0) then we have case (II).

eval(p1, e[0] = (1/4)*e[2]^2); factor(%); solve(%, w)

w^4+w^2*e[2]+(1/4)*e[2]^2

(1/4)*(2*w^2+e[2])^2

(1/2)*(-2*e[2])^(1/2), -(1/2)*(-2*e[2])^(1/2), (1/2)*(-2*e[2])^(1/2), -(1/2)*(-2*e[2])^(1/2)

The other way for the discriminant to be zero is e[0] = 0. So the other cases (IV), (V), and (IX) with multiple roots but not perfect squares must arise from e[0] = 0 and e[0] <> e[2]^2/4.

Now if e[0]=0 then two of the roots are zero anyway and the other two are +/- sqrt(-e[2])

factor(eval(p1, e[0] = 0)); solve(%, w)

w^2*(w^2+e[2])

0, 0, (-e[2])^(1/2), -(-e[2])^(1/2)

So if e[2] = 0 we have the case (II) we saw already. If e[2] > 0 we have two roots zero and two imaginary which is case (V) with w1 = w2 = 0.

If e[2] < 0 then we have case (IX) with w1 =0, w2 = -w3 = sqrt(-e[2]).

We have exhausted all cases; case (IV) cannot arise iin this paper where e[1]=0.

 

Note that D3 given in the paper simplifies for e[1] = 0 to

D2 := factor(eval(-2*e[2]^3+8*e[0]*e[2]-9*e[1]^2, e[1] = 0))

2*e[2]*(-e[2]^2+4*e[0])

so if e[0] = e[2]^2/4 then both D3 and D4 are zero.
Note also that D2 = -e[2] so we can check that with the above.

NULL

Download DiscrimChain.mw

I see now you flipped the sign of c in your substitution, to make every thing consistently x+y-c*t.So the authors consistently used x+y+c*t. To get the plot orientation right you need x+y+c*t or use your convention and a negative c value in producing the plot.

Here it is altogether - you can change the line with xtc to change the convention.

derivation.mw

@salim-barzani Excellent. You used xi = x+y-c*t. I figured out that for the plots you need xi = x+y+c*t to get the correct orientation, but then pde2 is not satisfied. So the authors must have used xi = x+y-c*t in deriving the substitution in 3.4 but xi = x+y+c*t elsewhere. It is in any case just a matter of changing the sign of c.

@salim-barzani You have u and v in pde1 and pde2 in a way that I didn't understand. The paper says xi = x+y+c*t but you have x+y-c*t (which seems better to me). Either way pde1 is satisfied but pde2 is not. Perhaps something is not right about the derivation in Eqs 3.1-3.5?

P5-pde.mw

Here's further back to psi and fig 1. I had v=u^3/3 instead of v=u^2/3 earlier; now everything looks correct.

Subcases (1) and (3)

restart

with(PDEtools)

undeclare(prime, quiet)

declare(u(xi), quiet); declare(w(`&xi;__1`), quiet)

Eq 3.5

ode := (diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+g[0]

(diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+g[0]

Eq 3.6 with epsilon = 1

itr := {`&xi;__1` = delta^(1/4)*xi, w(`&xi;__1`) = delta^(1/4)*u(xi)}

{xi__1 = delta^(1/4)*xi, w(xi__1) = delta^(1/4)*u(xi)}

tr := solve(itr, {xi, u(xi)})

{xi = xi__1/delta^(1/4), u(xi) = w(xi__1)/delta^(1/4)}

Eq. 3.7

ode2 := dchange(tr, ode, [`&xi;__1`, w(`&xi;__1`)], params = {delta})

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+(beta*p^3+p*q-r)*w(xi__1)^2/(delta^(1/2)*(3*beta*p+1))+g[0]

Eq 3.7 in terms of e2 and e0

ode3 := subs(1/sqrt(delta) = e__2*(3*beta*p+1)/(beta*p^3+p*q-r), g[0] = e__0, ode2)

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+e__2*w(xi__1)^2+e__0

F := subs(w(`&xi;__1`) = W, rhs(ode3))

W^4+W^2*e__2+e__0

F2 := collect(((W-w__1)^2+w__2^2)^2, W)

W^4-4*w__1*W^3+(6*w__1^2+2*w__2^2)*W^2-4*(w__1^2+w__2^2)*w__1*W+(w__1^2+w__2^2)^2

 

For the coefficients of w^3 to be the same we need w1 = 0. Then for the coeffs of w^2 and 1 we need to solve

F21 := eval(F2, w__1 = 0); eqs := {coeff(F21, W^2) = coeff(F, W^2), coeff(F21, W, 0) = coeff(F, W, 0)}

W^4+2*W^2*w__2^2+w__2^4

{w__2^4 = e__0, 2*w__2^2 = e__2}

There is no solution for w2 unless we have some special relationship between e[2] and e[0]

eliminate(eqs, w__2); eq := %[1][2][]; eq0 := e__0 = solve(eq, e__0)

[{w__2 = -(1/2)*2^(1/2)*e__2^(1/2)}, {-e__2^2+4*e__0}], [{w__2 = (1/2)*2^(1/2)*e__2^(1/2)}, {-e__2^2+4*e__0}]

-e__2^2+4*e__0

e__0 = (1/4)*e__2^2

So with this substitution we can work out the integral. The result has a slightly different form since w1 =0

factor(eval(F, eq0)); `assuming`([Int(1/sqrt(%), W)], [2*W^2+e__2 > 0]); eqxi1 := `&xi;__1` = value(%)+`&xi;__10`

(1/4)*(2*W^2+e__2)^2

Int(1/(W^2+(1/2)*e__2), W)

xi__1 = 2^(1/2)*arctan(W*2^(1/2)/e__2^(1/2))/e__2^(1/2)+xi__10

So the solution to ode3 for this is

sol3 := w(`&xi;__1`) = solve(eqxi1, W)

w(xi__1) = (1/2)*tan((1/2)*e__2^(1/2)*(xi__1-xi__10)*2^(1/2))*e__2^(1/2)*2^(1/2)

odetest(sol3, eval(ode3, eq0))

0

Convert back to earlier constants. We have e0 = g[0] =e2^2/4 so change

eqe2 := e__2 = (beta*p^3+p*q-r)/(sqrt(delta)*(3*beta*p+1)); eqg0 := g[0] = rhs(eval(eq0, eqe2)); ode2b := eval(ode2, eqg0); sol2 := simplify(eval(sol3, eqe2))

e__2 = (beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1))

g[0] = (1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

(diff(w(xi__1), xi__1))^2 = w(xi__1)^4+(beta*p^3+p*q-r)*w(xi__1)^2/(delta^(1/2)*(3*beta*p+1))+(1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

w(xi__1) = (1/2)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(xi__1-xi__10)*2^(1/2))*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*2^(1/2)

odetest(sol2, ode2b)

0

So now transform back to u(xi)

itr2 := {`&xi;__10` = delta^(1/4)*`&xi;__0`}

solu := dchange(`union`(itr, itr2), sol2/delta^(1/4), [u(xi), xi, `&xi;__0`], params = {delta}); odeb := eval(ode, eqg0)

u(xi) = (1/2)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*xi-delta^(1/4)*xi__0)*2^(1/2))*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*2^(1/2)/delta^(1/4)

(diff(u(xi), xi))^2 = delta*u(xi)^4+(beta*p^3+p*q-r)*u(xi)^2/(3*beta*p+1)+(1/4)*(beta*p^3+p*q-r)^2/(delta*(3*beta*p+1)^2)

odetest(solu, odeb)

0

So v(xi) is, from Eq. 3.2

solv := v(xi) = eval(-(1/3)*u(xi)^2, solu)

v(xi) = -(1/6)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*xi-delta^(1/4)*xi__0)*2^(1/2))^2*(beta*p^3+p*q-r)/(delta*(3*beta*p+1))

And so psi is perhaps? (There is something about taking the real part I didn't understand.)

psi := MakeFunction(eval(rhs(solv), xi = c*t+x+y), x, y, t)

proc (x, y, t) options operator, arrow; -(1/6)*tan((1/2)*((beta*p^3+p*q-r)/(delta^(1/2)*(3*beta*p+1)))^(1/2)*(delta^(1/4)*(c*t+x+y)-delta^(1/4)*xi__0)*2^(1/2))^2*(beta*p^3+p*q-r)/(delta*(3*beta*p+1)) end proc

Figure 1 (psi may be scaled differently). Not sure which is x and which is t

psi1 := eval(psi(x, 0, t), {`&xi;__0` = 0, beta = 1, c = 1/4, delta = 1/6, p = -1/2, q = 2*sqrt(6)*(1/3), r = -1/8})

-(2/3)*tan((1/12)*4^(1/2)*6^(3/4)*((1/4)*t+x)*2^(1/2))^2*6^(1/2)

plot3d(psi1, t = -10 .. 10, x = -10 .. 10, view = [default, default, -40 .. 2], orientation = [55, 55])

plot(eval(psi1, t = 0), x = -10 .. 10, -40 .. 2)

Subcase (3) is the same function, just written in a different way. Fig. 3:

psi4 := eval(psi(x, 0, t), {`&xi;__0` = 0, beta = 1, c = 1/4, delta = 1/6, p = -1/2, q = -(1/2)*sqrt(6), r = -1/8})

(1/2)*tan((1/12)*(-3)^(1/2)*6^(3/4)*((1/4)*t+x)*2^(1/2))^2*6^(1/2)

plot3d(psi4, t = -5 .. 5, x = -5 .. 5, view = [default, default, -2 .. 0], orientation = [55, 55])

NULL

Download fp3.mw

@salim-barzani Here's how to transform back to Eq 3.5, which was the original starting point of this question

fp3.mw

I think the subcases in the paper are not all different when e1=0, so it is simpler to go through them without the classification scheme. Knowing when the quartic is a perfect square is useful to simplify the integral, but it is not clear to me that knowing when the various roots are real or not has much to do with evaluating the integral. Probably the solutions can be interconverted. In the end, Maple can just solve the ode directly, so I'm not sure I can see the point.

roots.mw

@Alfred_F That's a nice observation. e0 is lost, but there are now two integration constants. But the square root is gone, which I think is more important.

Here's a simpler way to find two of the solutions, and of course the general solution is easy to find.

roots2.mw

@dharr I found a slightly different tan solution which only works for a special relationship between e0 and e2, so probably that is what is meant. 

fp3.mw

@salim-barzani OK, lets take epsilon=1 then I agree with Eq 3.7. I needed W because you can't integrate wrt a function in Maple, i.e. int(..., w(xi__1)) doesn't work.

But as I showed there is no way to transform F to F2, so the integral that gives the arctan result is invalid, and so Eq 3.13 is incorrect, so it will not satisfy the ode. You can perfom the integral with e2 and e0 in (just change the Int to int). The problem is then that you get xi__1 = f(w) but here and in general f(w) has multiple w in, so you cannot easily invert it to get w as a function of xi__1.

The solution to this is to solve the ode in Eq 3.7. dsolve will do this directly, or you could take the square root and integrate. Since you get a general solution it will work for arbitrary e2 and e0. If you want to then find cases that are real for various choices of e2 and e0 you can make assumptions that make the things under the square roots positive.

@salim-barzani Almost there.

@salim-barzani Well if it were changed I might look at it. I have several times pointed out that it makes it harder to help you, and I have ignored several of your recent questions for this reason. In the present case using w instead of r will be especially difficult since the paper uses w later for a different purpose.

I notice that as usual you changed the notation in your worksheet to make it harder to follow the paper. Perhaps you wish to correct that.

@sand15 Sorry; I must have google searched it instead of clicking the link

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