dharr

Dr. David Harrington

9484 Reputation

22 Badges

22 years, 39 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

Someone else found a similar problem here and I reported this bug to Maple. The same workaround might work here - put simplify(...) around the expression with the square root.

You can use the output=Array option of dsolve to do this.

> 

restart

Specify required x values

> 

xpts := Array([seq(0 .. 3.0, numelems = 20)])

Array(%id = 36893490082029227172)

> 

ode := {diff(y(x), x) = -a*y(x), y(0) = 1}

{diff(y(x), x) = -a*y(x), y(0) = 1}

> 

NULL

> 

ans20 := dsolve(eval(ode, a = 2.0), numeric, output = xpts); ans21 := dsolve(eval(ode, a = 2.1), numeric, output = xpts)

extract vectors of x and y values

> 

xx := ans20[2, 1][() .. (), 1]; y20 := ans20[2, 1][() .. (), 2]; y21 := ans21[2, 1][() .. (), 2]

Plot the difference

> 

plot(xx, y21-y20)

Plot the ratio (enter \/~)

> 

plot(xx, `~`[`/`](y21, y20))

Same thing using the new (Maple 2024) elementwise

> 

plot(xx, elementwise(y21/y20))

> 

NULL

Download dsolveArray.mw

If you expand Maple's result you get (1/4)*x^4 + x^2 + 1, differing from your expected value by a constant. Since you and Maple are both omitting arbitrary constants, both answers are correct.

@KIRAN SAJJAN I would write the 6th bc as (D@@2)(f)(1) = 0. But the equations have a third order derivative of f and a second order derivative of theta, so Maple is right in only expecting 5 boundary conditions. So the paper has an error, perhaps in the odes or perhaps in the bcs (redundant one perhaps?, but I didn't check.)

I recommend putting the parameters in a set so you can look at the odes without the numbers in to check them. Then it is easier to find typos. I found two. The second lambda_val just before Bstar in the second equation should be omitted, and the first equation has an incorrect term. Corrected here:

thin_film_base_paper_comparision.mw

[Updated] If I understand you correctly, then this classic method should work.

> 

restart;

> 

with(StringTools):

> 

s:="from -\\infty to \\infty and then from \\infty to -\\infty";

"from -\infty to \infty and then from \infty to -\infty"

> 

SubstituteAll(s, "-\\infty", "#@#"):
SubstituteAll(%, "\\infty", "+\\infty"):
s2:=SubstituteAll(%, "#@#","-\\infty");

"from -\infty to +\infty and then from +\infty to -\infty"

> 

NULL

Download Substitute.mw

The Heunc help page has some detailed instructions for finding polynomial solutions.

> 

restart;

> 

local gamma:

> 

Hc := HeunC(-4*b, 4*b, 2, 8*b^2, -8*b^2 - _c[2] + 1, (-r + 2*b)/(2*b));

HeunC(-4*b, 4*b, 2, 8*b^2, -8*b^2-_c[2]+1, (1/2)*(-r+2*b)/b)

> 

alpha, beta, gamma, delta, eta := op(1..5, Hc);

-4*b, 4*b, 2, 8*b^2, -8*b^2-_c[2]+1

A necessary condition for a polynomial solution of degree n is (see the HeunC help page)

> 

poly_cond := delta = -(n + (gamma + beta + 2)/2)*alpha;

8*b^2 = 4*(n+2+2*b)*b

So b has to be (for any n)

> 

solve(poly_cond, b);

0

But the last argument has b in the denominator, so b=0 does not work, and there can be no polynomial solutions

> 

eval(Hc,b=0);

Error, numeric exception: division by zero

> 

NULL

Download polysols.mw

There doesn't seem to be a solution, assuming the data you give is exactly correct. The presence of 1+A*3.7653333333333333333*10^(-18)*nu needs very large A for the types of nu values you have, otherwise the A*c*nu term is too small compared to 1 to make a difference.

MLsolve.mw

[Edit: I give a more general answer below.]

> 

restart;

> 

local gamma;

gamma

If we know that a is in the field generated by b and we have each of them written as RootOf  a single polynomial

> 

a := RootOf(4*_Z^3 - 12*_Z^2 - 24*_Z - 25, index = 1); # which actually equals “(b + 1)^2”, but assuming we do not know this beforehand
b := RootOf(2*_Z^3 - 3, index = 1);

RootOf(4*_Z^3-12*_Z^2-24*_Z-25, index = 1)

RootOf(2*_Z^3-3, index = 1)

> 

alias(alpha=a,beta=b);

alpha, beta

Then we "just" evaluate a specifying b as a field extension. Unfortunately, many of the evala subcommands such as Normal deduce the field extension from the coefficients, and don't allow us to specify it. Factor allows us to specify it, and even though we don't need to actually factor here it will serve the purpose:

> 

evala(Factor(a,b));

beta^2+2*beta+1

We don't need Primfield, which just confirms that a can be written in terms of b, so the field extension can just be b

> 

evala(Primfield({a,b}));

[[beta = beta], [beta = beta]]

Algfield works similarly, but gives the relationship between a and b

> 

evala(Algfield({a,b}));

[[alpha = beta^2+2*beta+1], [], {beta}, true, {}]

In the second case, Maple already simplifies each of these in terms of RootOf(_Z^2-2, index = 1) and RootOf(_Z^2+1, index = 1), which doesn't help us here.

> 

c:= convert((-1)^(1/4), RootOf);
evalf(%);
d:= convert((sqrt(2) - 1)*I, RootOf); # which actually equals c*(c^2 - c + 1), but assuming we do not know this beforehand
evalf(%);

(1/2)*RootOf(_Z^2-2, index = 1)+(1/2)*RootOf(_Z^2-2, index = 1)*RootOf(_Z^2+1, index = 1)

.7071067810+.7071067810*I

RootOf(_Z^2+1, index = 1)*(RootOf(_Z^2-2, index = 1)-1)

.4142135620*I

Convert these to forms that are RootOf a single polynomial

> 

mp:=evala(Minpoly(c,x)):
c:=select(q->evala(c-q)=0,[seq(RootOf(mp,index=i),i=1..degree(mp,x))])[];
evalf(%);
mp:=evala(Minpoly(d,x)):
d:=select(q->evala(d-q)=0,[seq(RootOf(mp,index=i),i=1..degree(mp,x))])[];
evalf(%);

RootOf(_Z^4+1, index = 1)

.7071067812+.7071067812*I

RootOf(_Z^4+6*_Z^2+1, index = 1)

.4142135624*I

> 

alias(gamma=c, delta=d);

alpha, beta, gamma, delta

> 

evala(Factor(d,c));

gamma^3-gamma^2+gamma

In this example, we can write RootOf(_Z^4+6*_Z^2+1) in terms of  RootOf(_Z^4+1, index = 1) The method doesn't work the other way around - it just returns c:

> 

evala(Factor(c,d));

gamma

Same problem for the first one if we want to do it the other way around.

So this method also does not work on the example in the original question

> 

evala(Factor(b,a));

beta

> 

NULL

Download evala.mw

Perhaps I missed something since age 21 in 1968 seems beyond a child. Not a very elegant solution.

"During a birthday party, the birthday child realizes: In 1968, I was the same age as the sum of the digits of my birth year. How old will I be now at the end of 2025?"

> 

restart

Birth year and its sum of digits

> 

b := 1900+10*d1+d2; s := 1+9+d1+d2

1900+10*d1+d2

10+d1+d2

Condition

> 

eq := 1968-b = s

68-10*d1-d2 = 10+d1+d2

Just search

> 

for d1 from 0 to 9 do for d2 from 0 to 9 do if eq then print(d1, d2, eval(2025-b)) end if end do end do

4, 7, 78

So born in 1947, age 21 in 1968 = 1+9+4+7; age in 2025:

> 

2025-1947

78

> 

NULL

Download Birthday.mw

It works using the simplify command rather than the context menu:

restartNULL

 

  k := 10.*Unit('N'/'m')

m := 2.*Unit('kg')NULL

omega := sqrt(k/m)"(->)"2.236067977*Units:-Unit(1/s)

 

"x(t):=A*cos(omega*t)+B*sin(omega*t):"

"v(t):=-A*omega*sin(omega*t)+B*omega*cos(omega*t):"

NULL

L1 := x(0*Unit('s')) = .15*Unit('m')

L2 := v(0*Unit('s')) = 0*Unit('m'/'s')

NULL

s := fsolve({L1, L2}, {A, B}) = {A = .15*Units:-Unit(m), B = -0.}NULL

assign(s)

 

x(t)"(=)".15*Units:-Unit(m)*cos(2.236067977*t*Units:-Unit(1/s))

v(t)

-.3354101966*Units:-Unit(m)*Units:-Unit(1/s^2)^(1/2)*sin(2.236067977*Units:-Unit(1/s^2)^(1/2)*t)

simplify(v(t))

-.3354101966*Units:-Unit(m/s)*sin(2.236067977*t*Units:-Unit(1/s))

v(t)"(=)"-.3354101966*Units:-Unit(m)*sin(2.236067977*t*Units:-Unit(1/s))

``

Download Error_with_units.mw

The change to upright font for D and O appears to be hardcoded into the last step of preparing the output, so diificult to fix (at least for users). 

> 

restart

> 

_local(D, O)

O

> 

a*D*O/sin(D)

a*D*O/sin(D)

The typesetting code for this.
Notice that both a and D are inside Typesetting:-mi(), which is for italics

> 

Typesetting:-Typeset(a*D*O/sin(D)); lprint(%)

a*D*O/sin(D)

Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mi("a"),Typesetting:-mo(
"⁢"),Typesetting:-mi("D"),Typesetting:-mo("⁢"),
Typesetting:-mi("O")),Typesetting:-mrow(Typesetting:-mi("sin",fontstyle =
"normal"),Typesetting:-mo("⁡"),Typesetting:-mfenced(Typesetting:-
mi("D"))))

Note different output here, suggesting that only in the final stage is there code to force D and O upright.

> 

Typesetting:-mi("D"); Typesetting:-mi("O"); Typesetting:-mi(a)

D

O

a

> 

NULL

Download Typesetting.mw

adding "assuming real" is a workaround for this case.

> 

interface(version);

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

> 

Physics:-Version();

`The "Physics Updates" package is not available for the version of Maple under development`

> 

SupportTools:-Version();

`The Customer Support Updates package is not installed. To install it, please run the command SupportTools:-Update().`

> 

 

Example fail when using _EnvAllSolutions := true:

 
> 

restart;

> 

kernelopts('assertlevel'=2):

> 

ode:=diff(y(x),x)=ln(1+y(x)^2);
IC:=y(0)=0;
x0:=0;
sol:=y(x) = -sqrt(-1 + exp(RootOf(-Intat(-1/(2*tau*sqrt(-1 + exp(tau))*exp(-tau)), tau = _Z) + x + _C2)));

eq:=0=eval(rhs(sol),x=x0);
_EnvAllSolutions := true:
_EnvExplicit := true:
sol_C:=_C2=solve(eq,_C2);
sol:=eval(sol,sol_C);
odetest(sol,[ode,IC]) assuming real;
simplify(%, symbolic);

diff(y(x), x) = ln(1+y(x)^2)

y(0) = 0

0

y(x) = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+x+_C2)))^(1/2)

0 = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+_C2)))^(1/2)

_C2 = Intat(-(1/2)*exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = (2*I)*Pi*_Z1)

y(x) = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+x+Intat(-(1/2)*exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = (2*I)*Pi*_Z1))))^(1/2)

[RootOf(Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = _Z)+2*x-Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = (2*I)*Pi*_Z1))-ln(exp(RootOf(Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = _Z)+2*x-Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = (2*I)*Pi*_Z1)))), 0]

[0, 0]

> 

 

Example Works  when using _EnvAllSolutions := false:

 
> 

restart;

> 

kernelopts('assertlevel'=2):

> 

ode:=diff(y(x),x) = ln(1+y(x)^2);
IC:=y(0)=0;
x0:=0;
sol:=y(x) = -sqrt(-1 + exp(RootOf(-Intat(-1/(2*tau*sqrt(-1 + exp(tau))*exp(-tau)), tau = _Z) + x + _C2)));

eq:=0=eval(rhs(sol),x=x0);
_EnvAllSolutions := false:
_EnvExplicit := true:
sol_C:=_C2=solve(eq,_C2);
sol:=eval(sol,sol_C);
odetest(%,[ode,IC])

diff(y(x), x) = ln(1+y(x)^2)

y(0) = 0

0

y(x) = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+x+_C2)))^(1/2)

0 = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+_C2)))^(1/2)

_C2 = Intat(-(1/2)*exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = 0)

y(x) = -(-1+exp(RootOf(-Intat(-(1/2)/(tau*(-1+exp(tau))^(1/2)*exp(-tau)), tau = _Z)+x+Intat(-(1/2)*exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = 0))))^(1/2)

[RootOf(Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = _Z)+2*x-Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = 0))-ln(exp(RootOf(Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = _Z)+2*x-Intat(exp(tau)/(tau*(-1+exp(tau))^(1/2)), tau = 0)))), 0]

> 

 

> 

 

Download odetest_fail_when_using_envAllsol_maple_2025_1_oct_21_2025.mw

Is this what you want? Note that when using collect, collecting wrt G{xi) handles negative powers so you do not need to collect wrt 1/G(xi).

> 

restart

> 

with(PDEtools)

> 

undeclare(prime, quiet)

> 

declare(U(xi), quiet)

> 

NULL

> 

ode := -2*(diff(diff(diff(diff(U(xi), xi), xi), xi), xi))*c-(108*c*(1/49)+8*U(xi)^2+1)*(diff(diff(U(xi), xi), xi))+18*U(xi)*(diff(U(xi), xi))^2+2*U(xi)*(-17*U(xi)^2*(1/49)+243*c*(1/2401)+9/98+w)

-2*(diff(diff(diff(diff(U(xi), xi), xi), xi), xi))*c-((108/49)*c+8*U(xi)^2+1)*(diff(diff(U(xi), xi), xi))+18*U(xi)*(diff(U(xi), xi))^2+2*U(xi)*(-(17/49)*U(xi)^2+(243/2401)*c+9/98+w)

> 

S := alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi)

alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi)

> 

WW := diff(G(xi), `$`(xi, 2)) = -lambda*G(xi)+mu

diff(diff(G(xi), xi), xi) = -lambda*G(xi)+mu

> 

NULL

> 

D1 := diff(S, xi)

> 

D11 := subs(WW, D1)

> 

D2 := diff(D11, xi)

> 

D22 := subs(WW, D2)

> 

D3 := diff(D22, xi)

> 

D33 := subs(WW, D3)

> 

D4 := diff(D33, xi)

> 

D44 := subs(WW, D4)

> 

K := U(xi) = S

> 

K1 := diff(U(xi), xi) = D11

> 

K2 := diff(diff(U(xi), xi), xi) = D22

> 

K3 := diff(diff(U(xi), xi), xi, xi) = D33

> 

K4 := diff(diff(U(xi), xi), xi, xi, xi) = D44

> 

P := eval(ode, {K, K1, K2, K3, K4})

-2*(alpha[1]*lambda^2*(diff(G(xi), xi))/G(xi)+15*alpha[1]*lambda*(-lambda*G(xi)+mu)*(diff(G(xi), xi))/G(xi)^2-20*alpha[1]*lambda*(diff(G(xi), xi))^3/G(xi)^3-60*alpha[1]*(-lambda*G(xi)+mu)*(diff(G(xi), xi))^3/G(xi)^4+30*alpha[1]*(-lambda*G(xi)+mu)^2*(diff(G(xi), xi))/G(xi)^3+24*alpha[1]*(diff(G(xi), xi))^5/G(xi)^5+24*beta[0]*(diff(G(xi), xi))^4/G(xi)^5-36*beta[0]*(-lambda*G(xi)+mu)*(diff(G(xi), xi))^2/G(xi)^4-8*beta[0]*lambda*(diff(G(xi), xi))^2/G(xi)^3+6*beta[0]*(-lambda*G(xi)+mu)^2/G(xi)^3+beta[0]*lambda*(-lambda*G(xi)+mu)/G(xi)^2)*c-((108/49)*c+8*(alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi))^2+1)*(-alpha[1]*lambda*(diff(G(xi), xi))/G(xi)-3*alpha[1]*(-lambda*G(xi)+mu)*(diff(G(xi), xi))/G(xi)^2+2*alpha[1]*(diff(G(xi), xi))^3/G(xi)^3+2*beta[0]*(diff(G(xi), xi))^2/G(xi)^3-beta[0]*(-lambda*G(xi)+mu)/G(xi)^2)+18*(alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi))*(alpha[1]*(-lambda*G(xi)+mu)/G(xi)-alpha[1]*(diff(G(xi), xi))^2/G(xi)^2-beta[0]*(diff(G(xi), xi))/G(xi)^2)^2+2*(alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi))*(-(17/49)*(alpha[0]+alpha[1]*(diff(G(xi), xi))/G(xi)+beta[0]/G(xi))^2+(243/2401)*c+9/98+w)

We want to find coefficients G^(n)/G^m. collect(P, {1/G(xi), diff(G(xi), xi)}, distributed)collects the way we want, but coeffs does not pick them off correctly (probably because it wasn't intended to work with diff.)

Write Gp for the derivative, and then collect (distributed) on the result.

> 

P2 := subs(diff(G(xi), xi) = Gp, P)

> 

P3 := collect(P2, {Gp, G(xi)}, distributed)

> 

eqs := {coeffs(P3, {Gp, G(xi)}, 'monomials')}

{2*beta[0]^3, -14*alpha[0]*alpha[1]^2, 6*beta[0]^2*alpha[1], 8*beta[0]^3*mu, 4*beta[0]^2*alpha[1]*mu, 4*alpha[0]*alpha[1]^2*lambda-(102/49)*alpha[0]*alpha[1]^2, 18*alpha[0]*alpha[1]^2*lambda^2+2*alpha[0]*(-(17/49)*alpha[0]^2+(243/2401)*c+9/98+w), -80*alpha[1]*lambda*c-2*((108/49)*c+8*alpha[0]^2+1)*alpha[1]+20*alpha[1]^3*lambda-(34/49)*alpha[1]^3, -32*alpha[1]*lambda^2*c-2*((108/49)*c+8*alpha[0]^2+1)*alpha[1]*lambda+18*alpha[1]^3*lambda^2-(68/49)*alpha[0]^2*alpha[1]+2*alpha[1]*(-(17/49)*alpha[0]^2+(243/2401)*c+9/98+w), -60*alpha[1]*mu^2*c+4*beta[0]^2*alpha[1]*lambda+28*beta[0]*alpha[0]*alpha[1]*mu+18*alpha[1]^3*mu^2-(102/49)*beta[0]^2*alpha[1], -56*beta[0]*lambda*c-2*((108/49)*c+8*alpha[0]^2+1)*beta[0]+32*beta[0]*alpha[1]^2*lambda+12*alpha[0]*alpha[1]^2*mu-(102/49)*beta[0]*alpha[1]^2, -12*beta[0]*mu^2*c-8*beta[0]^3*lambda+16*beta[0]^2*alpha[0]*mu+18*beta[0]*alpha[1]^2*mu^2-(34/49)*beta[0]^3, 90*alpha[1]*lambda*mu*c-12*beta[0]*alpha[0]*alpha[1]*lambda+3*((108/49)*c+8*alpha[0]^2+1)*alpha[1]*mu-36*alpha[1]^3*mu*lambda-(204/49)*beta[0]*alpha[0]*alpha[1], -10*beta[0]*lambda^2*c-((108/49)*c+8*alpha[0]^2+1)*beta[0]*lambda+18*beta[0]*alpha[1]^2*lambda^2-36*alpha[0]*alpha[1]^2*mu*lambda+2*beta[0]*(-(17/49)*alpha[0]^2+(243/2401)*c+9/98+w)-(68/49)*alpha[0]^2*beta[0], 22*beta[0]*mu*lambda*c-16*beta[0]^2*alpha[0]*lambda+((108/49)*c+8*alpha[0]^2+1)*beta[0]*mu-36*beta[0]*alpha[1]^2*mu*lambda+18*alpha[0]*alpha[1]^2*mu^2-(102/49)*beta[0]^2*alpha[0], 2*alpha[1]^3-48*c*alpha[1], 6*alpha[1]^2*beta[0]-48*c*beta[0], -16*mu*alpha[1]^2*beta[0]+72*c*mu*beta[0]-14*alpha[0]*beta[0]^2, -12*mu*alpha[1]^3+120*c*mu*alpha[1]-28*alpha[0]*alpha[1]*beta[0]}

> 

monomials

1, Gp^4/G(xi)^4, Gp^4/G(xi)^5, Gp^3/G(xi)^3, Gp^3/G(xi)^4, Gp^3/G(xi)^5, Gp^2/G(xi)^2, Gp^2/G(xi)^3, Gp^2/G(xi)^4, Gp^2/G(xi)^5, Gp/G(xi), Gp/G(xi)^2, Gp/G(xi)^3, Gp/G(xi)^4, 1/G(xi), 1/G(xi)^2, 1/G(xi)^3, 1/G(xi)^4, Gp^5/G(xi)^5

> 

``

NULL

Download system.mw

You had several omissions/errors in the PDE and H,E,X (in red). As usual there doesn't seem a universal method; here one needs to think about the sech/tanh relationship.

Download pp.mw

try params = {k}. If it's not that simple, please upload your worksheet using the Mapleprimes green up-arrow,

4 5 6 7 8 9 10 Last Page 6 of 89