dharr

Dr. David Harrington

9067 Reputation

22 Badges

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

I guess my point was that I would not expect eval to work with infinity, and I always use limit with infinity (for the reason that @sand15 mentioned: they are different things). Elsewhere in Maple infinity is handled slightly differently to mathematics, e.g., a*infinity is left as is until simplify(a*infinity) assuming a>0 leads to infinity. But I think this is fairly well thought out and generally doesn't lead to inconsistencies. Another principle that is often found in Maple is that if something is going to be fairly compute intensive it may not be done by default but you can force it with an option or another command. Here I do not think you want to overload eval with special cases since then writing efficient code for long calculations in a loop degrades performance.

@C_R @nm I guess with infinity it is safer to routinely use limit.

A workaround is to use expand

eval(expand((A*x-1)/x),x=infinity);

But it is hard to see any reason to get zero here. Many other places Maple gives sensible results with infinity. For example, infinity/infinity gives undefined, which would be a better result here than zero.

The fractional powers are very hard to deal with. I think you have to transform your ode to get rid of them

simplify(PDEtools:-dchange(V(xi) = U(xi)^3, Fode, [U(xi)])) assuming U(xi) >0

Of course if you really want V(xi) negative there is a mathematical problem about which cube root you want, but I'm guessing the equation was produced with the same disregard for rigor.

I'm struggling to understand the meaning here. It seems you want x=1 at t=0 since the units of the upper and lower limits should be the same and the integral should go to zero at t=0 (from the rhs). But in your approximation for x(t) you have ln(t) which implies you want t=1 to be the initial time?

In your mimimization you minimized t. Probably you meant to mimimize the integral of the squared difference over some range of t? 

Under my guessed interpretation, there is an easy solution:

intsolve.mw

These are not files listed on the File format help page, so you are out of luck for an easy way. I guess if you know the details of the binary file you could write a routine to read them, but I don't think that would count as "easy". 

Just installed 2026.0 in windows 11; run as administrator to install; didn't import previous preferences. All is working fine.

@KIRAN SAJJAN A quick look suggests it is reasonable, assuming the R1(z) and R2(z) functions are correct. I believe you said the BCs depended on z, but I do not see that dependence, though it might be implicit.

@Suryakanth @KIRAN SAJJAN You set up the odes for other conditions, so why not set up for different R or z and solve them? Sorry, I don't understand what you want.

@KIRAN SAJJAN Sorry, the equations you provided do not contain z, so I'm confused. Do you mean it is another independent variable, and you want to solve a PDE? (I don't want to read and figure out the paper.)

@WD0HHU This seemed to be fixed in 2025.2

@salim-barzani Use a set for the 3 plots otherwise it will think it is a parametric plot.

Dr.D-N.mw

You can directly use the procedure for w(eta) returned by dsolve to evaluate w(0.999) etc and the integrals.

parentartery_and_daughter_artery_error.mw

@salim-barzani That explains why odetest didn't work; the rest is OK. The sum of kinetic energy and potential energy should be approximately constant, not each separately.

You undoubtedly know some better numerical values for the parameters. Here the variation in energy is about 20%. You could also just plot the approximate solution and the exact solution (which you know?) together and see how they look.

restart

with(PDEtools)

undeclare(prime, quiet)

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

pde1 := I*(diff(Omega(x, y, t), t))-I*(diff(Omega(x, y, t), x))+alpha[1]*(diff(Omega(x, y, t), `$`(x, 2)))+alpha[2]*(diff(Omega(x, y, t), `$`(y, 2)))+alpha[3]*(diff(Omega(x, y, t), x, y))-alpha[4]*Omega(x, y, t)*U(p*x+q*y-t*v)^2

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

T := xi = p*x+q*y-t*v; T1 := Omega(x, y, t) = U(p*x+q*y-t*v)*exp(I*(-r*t+theta*y+w*x+phi))

P1 := I*(diff(Omega(x, y, t), t))-I*(diff(Omega(x, y, t), x))+alpha[1]*(diff(Omega(x, y, t), `$`(x, 2)))+alpha[2]*(diff(Omega(x, y, t), `$`(y, 2)))+alpha[3]*(diff(Omega(x, y, t), x, y))-alpha[4]*Omega(x, y, t)*U(p*x+q*y-t*v)^2

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

Eq. 37

ode1 := expand(lhs(ode)/exp(I*(-r*t+theta*y+w*x+phi)))

ode11 := evalc(Re(ode1))

ode111 := collect(ode11, {U(xi), diff(diff(U(xi), xi), xi)})

-U(xi)^3*alpha[4]+(-theta^2*alpha[2]-theta*w*alpha[3]-w^2*alpha[1]+r+w)*U(xi)+(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])*(diff(diff(U(xi), xi), xi))

Eq, 5.1

ode1111 := collect(ode111/DEtools:-dcoeffs(ode111, U(xi))[1], U(xi))

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) = ode1111*(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)

From which we can extract K and P.

K := op(1, integrand); P := collect(-op(2, integrand), U(xi))

Eq. 5.6

H := K+P

(1/2)*(diff(U(xi), xi))^2+alpha[4]*U(xi)^4/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*U(xi)^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

Eq. 5.7

sol1 := U(xi) = m*cos(delta*xi)

U(xi) = m*cos(delta*xi)

Eq. 5.10

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

Eq. 5.11 - looks correct

H1 := eval(H, sol1)-H0

Eq. 5.12 - looks correct

H2 := collect(algsubs(delta*xi = (1/4)*Pi, H1), delta, proc (x) options operator, arrow; collect(x, m) end proc)

Eq. 5.13

delta = sort([solve(H2, delta)])[2]; eqphi := `assuming`([simplify(%)], [p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2] > 0])

delta = (1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*alpha[4]*m^2+4*theta^2*alpha[2]-4*r)^(1/2)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2)

sol2 := simplify(eval(sol1, eqphi))

U(xi) = m*cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*alpha[4]*m^2+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))

So the kinetic, potential and total energies are

Kapprox := simplify(eval(K, sol2)); Papprox := simplify(eval(P, sol2)); Happrox := simplify(eval(H, sol2))

Which we expect to be approximately constant. Test with some sensible numerical values - her I just choose them all 1

params := `~`[`=`](`minus`(indets(Happrox, name), {xi}), 1)

{m = 1, p = 1, q = 1, r = 1, theta = 1, w = 1, alpha[1] = 1, alpha[2] = 1, alpha[3] = 1, alpha[4] = 1}

Knum := eval(Kapprox, params); Pnum := eval(Papprox, params); Hnum := eval(Happrox, params)

7/48-(7/48)*cos((1/3)*3^(1/2)*7^(1/2)*xi)

(1/12)*cos((1/6)*3^(1/2)*7^(1/2)*xi)^2*(cos((1/6)*3^(1/2)*7^(1/2)*xi)^2+2)

(1/12)*cos((1/6)*3^(1/2)*7^(1/2)*xi)^4-(1/8)*cos((1/6)*3^(1/2)*7^(1/2)*xi)^2+7/24

Total energy varies by about 0.05 compared to K and P p-p amplitudes of 0.25, or 20%.

plot([Knum, Pnum, Hnum], xi = 0 .. 5, view = [default, default], color = [red, blue, black])

NULL

Download u-1.mw

@salim-barzani Most equations look correct but odetest fails - maybe check some assumptions or try some numerical checks, or maybe the theory doesn't work.

restart

with(PDEtools)

undeclare(prime, quiet)

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

pde1 := I*(diff(Omega(x, y, t), t))-I*(diff(Omega(x, y, t), x))+alpha[1]*(diff(Omega(x, y, t), `$`(x, 2)))+alpha[2]*(diff(Omega(x, y, t), `$`(y, 2)))+alpha[3]*(diff(Omega(x, y, t), x, y))-alpha[4]*Omega(x, y, t)*U(p*x+q*y-t*v)^2

I*(diff(Omega(x, y, t), t))-I*(diff(Omega(x, y, t), x))+alpha[1]*(diff(diff(Omega(x, y, t), x), x))+alpha[2]*(diff(diff(Omega(x, y, t), y), y))+alpha[3]*(diff(diff(Omega(x, y, t), x), y))-alpha[4]*Omega(x, y, t)*U(p*x+q*y-t*v)^2

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

T := xi = p*x+q*y-t*v; T1 := Omega(x, y, t) = U(p*x+q*y-t*v)*exp(I*(-r*t+theta*y+w*x+phi))

P1 := I*(diff(Omega(x, y, t), t))-I*(diff(Omega(x, y, t), x))+alpha[1]*(diff(Omega(x, y, t), `$`(x, 2)))+alpha[2]*(diff(Omega(x, y, t), `$`(y, 2)))+alpha[3]*(diff(Omega(x, y, t), x, y))-alpha[4]*Omega(x, y, t)*U(p*x+q*y-t*v)^2

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

-alpha[4]*U(xi)^3*exp(I*(-r*t+theta*y+w*x+phi))+exp(I*(-r*t+theta*y+w*x+phi))*(-alpha[1]*w^2-(theta*alpha[3]-1)*w-theta^2*alpha[2]+r)*U(xi)+exp(I*(-r*t+theta*y+w*x+phi))*(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])*(diff(diff(U(xi), xi), xi))+I*exp(I*(-r*t+theta*y+w*x+phi))*((theta*alpha[3]+2*w*alpha[1]-1)*p+2*q*theta*alpha[2]+q*w*alpha[3]-v)*(diff(U(xi), xi)) = 0

Eq. 37

ode1 := expand(lhs(ode)/exp(I*(-r*t+theta*y+w*x+phi)))

I*(diff(U(xi), xi))*p*theta*alpha[3]+(2*I)*(diff(U(xi), xi))*p*w*alpha[1]+(2*I)*(diff(U(xi), xi))*q*theta*alpha[2]+I*(diff(U(xi), xi))*q*w*alpha[3]-U(xi)^3*alpha[4]-U(xi)*theta^2*alpha[2]-U(xi)*theta*w*alpha[3]-U(xi)*w^2*alpha[1]+(diff(diff(U(xi), xi), xi))*p^2*alpha[1]+(diff(diff(U(xi), xi), xi))*p*q*alpha[3]+(diff(diff(U(xi), xi), xi))*q^2*alpha[2]-I*(diff(U(xi), xi))*p-I*(diff(U(xi), xi))*v+U(xi)*r+U(xi)*w

ode11 := evalc(Re(ode1))*o

-U(xi)^3*alpha[4]-U(xi)*theta^2*alpha[2]-U(xi)*theta*w*alpha[3]-U(xi)*w^2*alpha[1]+(diff(diff(U(xi), xi), xi))*p^2*alpha[1]+(diff(diff(U(xi), xi), xi))*p*q*alpha[3]+(diff(diff(U(xi), xi), xi))*q^2*alpha[2]+U(xi)*r+U(xi)*w

ode111 := collect(ode11, {U(xi), diff(diff(U(xi), xi), xi)})

-U(xi)^3*alpha[4]+(-theta^2*alpha[2]-theta*w*alpha[3]-w^2*alpha[1]+r+w)*U(xi)+(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])*(diff(diff(U(xi), xi), xi))

Eq, 5.1

ode1111 := collect(ode111/DEtools:-dcoeffs(ode111, U(xi))[1], U(xi))

-U(xi)^3*alpha[4]/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])+(-theta^2*alpha[2]-theta*w*alpha[3]-w^2*alpha[1]+r+w)*U(xi)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])+diff(diff(U(xi), xi), xi)

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) = ode1111*(diff(U(xi), xi))

diff(f(xi), xi) = (-U(xi)^3*alpha[4]/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])+(-theta^2*alpha[2]-theta*w*alpha[3]-w^2*alpha[1]+r+w)*U(xi)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])+diff(diff(U(xi), xi), xi))*(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)*alpha[4]*U(xi)^4+(1/2)*(-p^2*alpha[1]-p*q*alpha[3]-q^2*alpha[2])*(diff(U(xi), xi))^2+(1/2)*(theta^2*alpha[2]+theta*w*alpha[3]+w^2*alpha[1]-r-w)*U(xi)^2)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])+_C1 = 0

(1/2)*(diff(U(xi), xi))^2-U(xi)^2*(alpha[4]*U(xi)^2+2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*alpha[1]*w^2-2*r-2*w)/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

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

K := op(1, integrand); P := collect(-op(2, integrand), U(xi))

(1/2)*(diff(U(xi), xi))^2

U(xi)^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*U(xi)^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

Eq. 5.6

H := K+P

(1/2)*(diff(U(xi), xi))^2+U(xi)^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*U(xi)^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

Eq. 5.7

sol1 := U(xi) = m*cos(delta*xi)

U(xi) = m*cos(delta*xi)

Eq. 5.10

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

(1/2)*(D(U))(0)^2+U(0)^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*U(0)^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

m^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*m^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

Eq. 5.11 - looks correct

H1 := eval(H, sol1)-H0

(1/2)*m^2*delta^2*sin(delta*xi)^2+m^4*cos(delta*xi)^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])+(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*m^2*cos(delta*xi)^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])-m^4*alpha[4]/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])-(2*theta^2*alpha[2]+2*theta*w*alpha[3]+2*w^2*alpha[1]-2*r-2*w)*m^2/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

Eq. 5.12 - looks correct

H2 := collect(algsubs(delta*xi = (1/4)*Pi, H1), delta, proc (x) options operator, arrow; collect(x, m) end proc)

(1/4)*m^2*delta^2-(3/16)*m^4*alpha[4]/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])-(1/16)*(4*theta^2*alpha[2]+4*theta*w*alpha[3]+4*w^2*alpha[1]-4*r-4*w)*m^2/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])

Eq. 5.13

delta = sort([solve(H2, delta)])[2]; eqphi := `assuming`([simplify(%)], [p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2] > 0])

delta = (1/2)*((p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])*(3*m^2*alpha[4]+4*theta^2*alpha[2]+4*theta*w*alpha[3]+4*w^2*alpha[1]-4*r-4*w))^(1/2)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])

delta = (1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2)

sol2 := simplify(eval(sol1, eqphi))

U(xi) = m*cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))

simplify(odetest(sol2, ode1111), symbolic)

-4*m*cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))*(cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))^2*m^2*alpha[4]+(3/4)*m^2*alpha[4]+2*alpha[1]*w^2+(2*theta*alpha[3]-2)*w+2*theta^2*alpha[2]-2*r)/(4*p^2*alpha[1]+4*p*q*alpha[3]+4*q^2*alpha[2])

simplify(odetest(sol2, ode111), symbolic)

-(cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))^2*m^2*alpha[4]+(3/4)*m^2*alpha[4]+2*alpha[1]*w^2+(2*theta*alpha[3]-2)*w+2*theta^2*alpha[2]-2*r)*m*cos((1/2)*(4*alpha[1]*w^2+(4*theta*alpha[3]-4)*w+3*m^2*alpha[4]+4*theta^2*alpha[2]-4*r)^(1/2)*xi/(p^2*alpha[1]+p*q*alpha[3]+q^2*alpha[2])^(1/2))

 

 

NULL

Download u-1.mw

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