aroche

Dr. Austin Roche

395 Reputation

7 Badges

12 years, 31 days
Maplesoft
Waterloo, Ontario, Canada
I am a Software Architect in the Math Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, in various areas including differential equations, integration, mathematical functions, simplification, root finding, and logic. I completed a Master's degree from McGill University with a thesis in Differential Geometry, and a PhD from Simon Fraser University with a thesis on Differential Equations.

MaplePrimes Activity


These are answers submitted by aroche

Hi @nm,

Thanks for this important report. Your analysis was spot on: The Cauchy-Euler solver was attempted despite the ODE not being in that class. The conditions were not being checked properly and that has now been fixed. The new result is an error telling that the ODE is not supported by the Student:-ODEs package:

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1840 and is the same as the version installed in this computer, created 2024, December 26, 12:58 hours Pacific Time.`

ode:=diff(diff(y(x),x),x)*sin(x)^2 = 2*y(x);
Student:-ODEs:-ODESteps(ode);

(diff(diff(y(x), x), x))*sin(x)^2 = 2*y(x)

Error, (in Student:-ODEs:-OdeSolveOrder2) ODE is not supported

Download NotEuler.mw


By the way, regarding the remark by @Alfred_F, dsolve gets an equivalent solution that needs a bit of massaging to get the desired form:

dsolve(ode);

y(x) = c__1*cot(x)+c__2*(I*ln(cos(2*x)+sin(2*x)*I)*sin(2*x)-2*cos(2*x)+2)/(-1+cos(2*x))

simplify(convert(%, exp), symbolic);

y(x) = (2*x*c__2+c__1)*cot(x)-2*c__2

Download dsolveNotEuler.mw


The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1840 or newer. To install the Updates, open Maple and input Physics:-Version(latest).

Austin Roche
Maplesoft

Hi @nm,

Thank you for this helpful report! The bug has been fixed and we now get the expected answer of [0,0] for the result of odetest comparing the solution with the ode and the IC:

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1839 and is the same as the version installed in this computer, created 2024, December 2, 10:11 hours Pacific Time.`

ode:=1+x*y(x)*(1+y(x)^2*x)*diff(y(x),x) = 0:
IC:=y(1) = 0:
sol:=x = 1/(3*exp(y(x)^2/2) - y(x)^2 - 2):

odetest(sol,ode);

0

odetest(sol,[ode,IC]);

[0, 0]

Download odetest_with_IC.mw

A note about what was going wrong: In the presence of an IC, implicit solutions in certain forms were wrongly being classified as explicit, and then the explicit tester was being used which was not effective. So in fact the order of the results was not switched! The 0 result was indeed a verification that the IC was satisfied, and the non-0 result was actually a failure to verify that the ODE itself was satisfied. This can be checked by providing a different IC as part of the test, and noticing that in that case both results are now non-0 (note that I am removing the fix to make this test):

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1839. The version installed in this computer is 1838 created 2024, December 2, 10:11 hours Pacific Time, found in the directory C:\Users\Austin\maple\toolbox\2024\Physics Updates\lib\`

ode:=1+x*y(x)*(1+y(x)^2*x)*diff(y(x),x) = 0:
IC:=y(1) = 0:
sol:=x = 1/(3*exp(y(x)^2/2) - y(x)^2 - 2):

odetest(sol,ode);

0

wrong_IC := y(1) = 1;

y(1) = 1

odetest(sol,[ode,wrong_IC]);

[(y(x)^4-y(x)^2*y(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))*(D(y))(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))+y(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))^3*(D(y))(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))-6*y(x)^2*exp((1/2)*y(x)^2)+3*exp((1/2)*y(x)^2)*y(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))*(D(y))(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))+4*y(x)^2-2*y(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))*(D(y))(-1/(y(x)^2-3*exp((1/2)*y(x)^2)+2))+9*exp(y(x)^2)-12*exp((1/2)*y(x)^2)+4)/(y(x)^2-3*exp((1/2)*y(x)^2)+2)^2, 1-(-2*LambertW(_Z2, -(3/2)*exp(-3/2))-3)^(1/2)]

Download previous_odetest_with_bad_IC.mw

The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1839 or newer. To install the Updates, open Maple and input Physics:-Version(latest).

Austin Roche
Maplesoft

Hi @nm,

Thanks again for finding this bug. The bug has been fixed and now after finding a polynomial solution via the series method, this particular solution is used to reduce the order and the reduced order ODE can also be solved, leading to the general solution of the original problem.

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1836 and is the same as the version installed in this computer, created 2024, December 2, 10:11 hours Pacific Time.`

ode:=3*diff(y(x),x$2)+x*diff(y(x),x)-4*y(x)=0;

3*(diff(diff(y(x), x), x))+x*(diff(y(x), x))-4*y(x) = 0

Student:-ODEs:-ODESteps(ode,y(x));

"[[,,"Let's solve"],[,,3 ((ⅆ)^2)/(ⅆx^2) y(x)+x ((ⅆ)/(ⅆx) y(x))-4 y(x)=0],["•",,"Highest derivative means the order of the ODE is" 2],[,,((ⅆ)^2)/(ⅆx^2) y(x)],["•",,"Isolate 2nd derivative"],[,,((ⅆ)^2)/(ⅆx^2) y(x)=-(x ((ⅆ)/(ⅆx) y(x)))/3+(4 y(x))/3],["•",,"Group terms with" y(x) "on the lhs of the ODE and the rest on the rhs of the ODE; ODE is linear"],[,,((ⅆ)^2)/(ⅆx^2) y(x)+(x ((ⅆ)/(ⅆx) y(x)))/3-(4 y(x))/3=0],["•",,"Multiply by denominators"],[,,3 ((ⅆ)^2)/(ⅆx^2) y(x)+x ((ⅆ)/(ⅆx) y(x))-4 y(x)=0],["•",,"Assume series solution for" y(x)],[,,y(x)=(∑)a[k] x^k],["▫",,"Rewrite DE with series expansions"],[,"?","Convert" x*((ⅆ)/(ⅆx) y(x)) "to series expansion"],[,,[]=(∑)a[k] k x^k],[,"?","Convert" ((ⅆ)^2)/(ⅆx^2) y(x) "to series expansion"],[,,((ⅆ)^2)/(ⅆx^2) y(x)=(∑)a[k] k (k-1) x^(k-2)],[,"?","Shift index using" k "->" k+2],[,,((ⅆ)^2)/(ⅆx^2) y(x)=(∑)a[k+2] (k+2) (k+1) x^k],[,,"Rewrite DE with series expansions"],[,,(∑)(3 a[k+2] (k+2) (k+1)+a[k] (k-4)) x^k=0],["•",,"Each term in the series must be 0, giving the recursion relation"],[,,(3 k^2+9 k+6) a[k+2]+a[k] (k-4)=0],["•",,"Recursion relation; series terminates at" k=4],[,,a[k+2]=-(a[k] (k-4))/(3 (k^2+3 k+2))],["•",,"Apply recursion relation for" k=0],[,,a[2]=(2 a[0])/3],["•",,"Apply recursion relation for" k=2],[,,a[4]=(a[2])/18],["•",,"Express in terms of" a[0]],[,,a[4]=(a[0])/27],["•",,"Terminating series solution of the ODE. Use reduction of order to find the second linearly independent solution"],[,,y(x)=a[0] (1+2/3 x^2+1/27 x^4)],["•",,"Use this particular solution to reduce the order of the ODE and find a complete solution"],[,,y(x)=(1+2/3 x^2+1/27 x^4) (∫((e)^(-(x^2)/6+c__1))/((x^4+18 x^2+27)^2) ⅆx+c__2)]]6""

Download steps.mw
The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1836 or newer. To install the Updates, open Maple and input Physics:-Version(latest).

Austin Roche
Maplesoft

 

 

Hi @nm,

This appears to be a problem in the RegularChains library and is currently under investigation.

Thank you for the report!

Austin Roche
Maplesoft

Hi Dmitry,

Thank you very much for this helpful report!

As mentioned by @ecterrab, there was an issue in `simplify/RootOf` that was encountered during the solving process for this problem. This issue has now been fixed, and the ode is now solvable via dsolve, returning a large answer. The answer is verifiable with odetest, and can be simplified to something of a more reasonable size by calling simplify:

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 1817 and is the same as the version installed in this computer, created 2024, October 3, 12:7 hours Pacific Time.`

sys := diff(S(t), t) = -b*(-alpha*u__1+c)*S(t)*K(t)/N-u__2*a*M, diff(K(t), t) = b*(-alpha*u__1+c)*S(t)*K(t)/N

inc := S(0) = S0, K(0) = K0

ans := dsolve({inc, sys})

`[Length of output exceeds limit of 1000000]`

odetest(ans, {inc, sys})

{0}

simplified_ans := simplify(ans)

length(simplified_ans)

77415

odetest(simplified_ans, {inc, sys})

{0}

Download dsolve_solve_RootOf_issue2.mw

Moreover, it should be noted that the alternative approach mentioned by @dharr actually leads to an even more convenient answer.

The fix is distributed for everybody using Maple 2024.1 within the Maplesoft Physics Updates v.1816 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

Hi,

Thank you for finding this bug! It's been fixed, but as mentioned by @Thomas Richard, the ODE is beyond the scope of what the Student package can handle, so the new result is just an error message saying that this ODE is not supported:

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1811 and is the same as the version installed in this computer, created 2024, September 21, 11:10 hours Pacific Time.`

ode:=sin(x)*diff(y(x),x$2)+cos(x)*diff(y(x),x)+(sin(x)-cos(x))*y(x)=0:

Student:-ODEs:-ODESteps(ode);

Error, (in Student:-ODEs:-OdeSolveOrder2) ODE is not supported

Download Physicsv1811.mw

The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1811 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

Hi,

Thanks for this report. A fix has been implemented internally and is expected to be avaiable in the next Physics update (edit: it's there already in v.1793, see comment below), as well as in the next full release of Maple.

The new answer is 'true'.

Austin Roche
Software Architect
Mathematical Software
Maplesoft

Hi @nm,

Just to let you know that I've made some changes and now the derivative is being solved for first, so for this example y(x)*y'(x) = 0 the solution is now the simpler y(x) = c as expected:

Student:-ODEs:-ODESteps(y(x)*diff(y(x),x)=0);

"[[,,"Let's solve"],[,,y(x) ((ⅆ)/(ⅆx) y(x))=0],["•",,"Highest derivative means the order of the ODE is" 1],[,,(ⅆ)/(ⅆx) y(x)],["•",,"Separate variables"],[,,(ⅆ)/(ⅆx) y(x)=0],["•",,"Integrate both sides with respect to" x],[,,∫((ⅆ)/(ⅆx) y(x)) ⅆx=∫0 ⅆx+`c__1`],["•",,"Evaluate integral"],[,,y(x)=`c__1`],["•",,"Solve for" y(x)],[,,y(x)=`c__1`]]"

Download 238382_Ans1.mw

As usual, the fix is included within the Maplesoft Physics Updates v.1758 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

Hi,

Thanks again for finding this bug. The problem here was again the fact that Student:-ODEs was relying on odeadvisor to check if the ODE was exact. In this case odeadvisor was saying that it was, but that was only true if it was first solved for the derivative:

ode:=diff(y(x), x)^3 = (x - 2)^2:

DEtools:-odeadvisor(ode, [exact]);

[_exact]

solved_ode := op(solve(ode, {diff(y(x),x)})[1]);

diff(y(x), x) = (x-2)^(2/3)

DEtools:-odeadvisor(solved_ode, [exact]);

[_exact]

Download StudentODEs_odeadvisor_beforefix.mw


I fixed this problem previously in the case that the input ODE was linear but not exact. I didn't realize it was also a problem in the non-linear case, but that should be solved now as well using the latest update:

ode:=diff(y(x), x)^3 = (x - 2)^2:

DEtools:-odeadvisor(ode, [exact]);

[NONE]

solved_ode := op(solve(ode, {diff(y(x),x)})[1]);

diff(y(x), x) = (x-2)^(2/3)

DEtools:-odeadvisor(%, [exact]);

[_exact]

ic:=y(2)=1:
Student:-ODEs:-ODESteps([ode,ic]);

"[[,,"Let's solve"],[,,[((ⅆ)/(ⅆx) y(x))^3=(x-2)^2,y(2)=1]],["•",,"Highest derivative means the order of the ODE is" 1],[,,(ⅆ)/(ⅆx) y(x)],["•",,"Separate variables"],[,,(ⅆ)/(ⅆx) y(x)=(x-2)^(2/3)],["•",,"Integrate both sides with respect to" x],[,,∫((ⅆ)/(ⅆx) y(x)) ⅆx=∫(x-2)^(2/3) ⅆx+`c__1`],["•",,"Evaluate integral"],[,,y(x)=(3 (x-2)^(5/3))/5+`c__1`],["•",,"Solve for" y(x)],[,,y(x)=(3 (x-2)^(5/3))/5+`c__1`],["•",,"Use initial condition" y(2)=1],[,,1=`c__1`],["•",,"Solve for" `c__1`],[,,`c__1`=1],["•",,"Substitute" `c__1`=1 "into general solution and simplify"],[,,y(x)=1+((3 x-6) (x-2)^(2/3))/5],["•",,"Solution to the IVP"],[,,y(x)=1+((3 x-6) (x-2)^(2/3))/5]]"

   

Download StudentODEs_odeadvisor.mw

The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1755 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

Hi,

Thanks again for finding this bug! The code was just attempting to evaluate the given solution at the given initial condition, which leads directly to an error in this case.  I've updated it to catch this type of error and instead declare that the given solution does not satisfy the given initial condition. Furthermore, in this particular case I've also added an extra step,  converting the arctanh to ln, that makes solving the general solution for y(x) easier. Once that is done, finding the value of the constant of integration becomes feasible. The new result is essentially the one expected:

ode:=diff(y(x), x) - 2*y(x) = 2*sqrt(y(x)):
ic:=y(0) = 1:
Student:-ODEs:-ODESteps([ode,ic]);

"[[,,"Let's solve"],[,,[(ⅆ)/(ⅆx) y(x)-2 y(x)=2 sqrt(y(x)),y(0)=1]],["•",,"Highest derivative means the order of the ODE is" 1],[,,(ⅆ)/(ⅆx) y(x)],["•",,"Separate variables"],[,,((ⅆ)/(ⅆx) y(x))/(2 y(x)+2 sqrt(y(x)))=1],["•",,"Integrate both sides with respect to" x],[,,∫((ⅆ)/(ⅆx) y(x))/(2 y(x)+2 sqrt(y(x))) ⅆx=∫1 ⅆx+`c__1`],["•",,"Evaluate integral"],[,,(ln(y(x)-1))/2+arctanh(sqrt(y(x)))=x+`c__1`],["•",,"Convert" arctanh "to 'ln'"],[,,(ln(y(x)-1))/2+(ln(sqrt(y(x))+1))/2-(ln(1-sqrt(y(x))))/2=x+`c__1`],["•",,"Solve for" y(x)],[,,{y(x)=-2 sqrt(-(e)^(2 x+2 `c__1`))-(e)^(2 x+2 `c__1`)+1,y(x)=2 sqrt(-(e)^(2 x+2 `c__1`))-(e)^(2 x+2 `c__1`)+1}],["•",,"Use initial condition" y(0)=1],[,,1=-2 sqrt(-(e)^(2 `c__1`))-(e)^(2 `c__1`)+1],["•",,"Solve for" `c__1`],[,,`c__1`=ln(2)+(ⅈ Pi)/2],["•",,"Substitute" `c__1`=ln(2)+(ⅈ Pi)/2 "into general solution and simplify"],[,,y(x)=-4 sqrt((e)^(2 x))+4 (e)^(2 x)+1],["•",,"Use initial condition" y(0)=1],[,,1=2 sqrt(-(e)^(2 `c__1`))-(e)^(2 `c__1`)+1],["•",,"Solve for" `c__1`],[,,`c__1`=()],["•",,"Solution does not satisfy initial condition"],["•",,"Solution to the IVP"],[,,y(x)=-4 sqrt((e)^(2 x))+4 (e)^(2 x)+1]]"

   

Download StudentODEsAnswered.mw
The fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1753 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

 

Hi,

Thanks for finding this bug!

As mentioned by other commenters, simplify always just maps into the operands of non-algebraic objects like equations. So there is no way to use simplify to manipulate an equation as you are trying to do here, but again as mentioned, convert could be a better option.

Nevertheless, the simplify extension mechanism was broken in the case where the name of the extension doesn't occur as a function in the input expression. So for example, this worked:

`simplify/F` := u -> G(u):
simplify(F(x), F);

G(F(x))

but this didn't:

simplify(x, F);

x

This issue has now been fixed. After installing the latest Physics Updates (v.1749), we now have:

`simplify/F` := u -> G(u):
eq := (a*x + b)/(c*x + d) = 11/3:
simplify(eq, F);

G((a*x+b)/(c*x+d)) = 11/3

eq := (a*x + b)/(c*x + d) = sqrt(2):
simplify(eq, F);

G((a*x+b)/(c*x+d)) = G(2^(1/2))

eq := (a*x + b)/(c*x + d) = y:
simplify(eq, F);

G((a*x+b)/(c*x+d)) = G(y)

Note that simplify also avoids applying any simplification to expressions which it deems are already simple enough (in particular that means anything of type complex(extended_numeric)), in this case the rhs of the first input 11/3 because it's numeric.

As an aside, I'll note that this problem with the extension mechanism first happened due to an improvement in simplify that was put in for Maple2023. Before that, if you called simplify with multiple extra arguments specifying sub-algorithms to be applied, it would just apply them once each, and this sometimes meant that by applying the same simplify command twice in succession, the second result would be an improvement:

 

e := RootOf(_Z^2 - exp(2*RootOf(_Z^2+1,index=1)*x) - 1):

simplify( diff(e-subs(RootOf(_Z^2+1,index=1) = I, e), x), RootOf, exp);

-I*(RootOf(_Z^2-exp((2*I)*x)-1)^2*exp((2*I)*x)-exp((2*I)*x)-exp((4*I)*x))/(RootOf(_Z^2-exp((2*I)*x)-1)*(exp((2*I)*x)+1))

simplify(%, RootOf, exp);

0

After that improvement was put in from Maple 2023, simplify will continue to alternate through the requested simplifications until the result stabilizes, which in this example means that simplify only needs to be called once to achieve the desired answer:

 

When multiple particular simplification procedures are requested via extra arguments to simplify, they will each be retried as many times as needed to obtain the full simplification. For example, the extra arguments exp and RootOf can now be given in either order to achieve the full simplification:

e := RootOf(_Z^2 - exp(2*RootOf(_Z^2+1,index=1)*x) - 1):

simplify( diff(e-subs(RootOf(_Z^2+1,index=1) = I, e), x), RootOf, exp);

0

Download simplify_ext.mw
Download simplify_ext2.mw
Download simplify_ext3.mw
Download simplify_ext4.mw

As usual, the fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1749 or newer. To install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

 

Hi,

Thank you for finding this bug!

The general form of the ODE where this problem occurred is this one:

ode := F(x,y(x))*diff(y(x),x) = 0;

F(x, y(x))*(diff(y(x), x)) = 0

It so happens that Student:-ODEs was relying on odeadvisor to classify the input ODE. Now, while a call to odeadvisor with no extra arguments produces the expected answer,

DEtools:-odeadvisor(ode);

[_quadrature]

Student:-ODEs was actually trying the 'exact' method first, and odeadvisor - for this very specific, trivially solvable class and when provided the extra argument [exact] - was asserting that indeed the ODE is exact.

DEtools:-odeadvisor(ode, [exact]);

[_exact]

This is of course not an exact ODE unless F has no dependence on x.

In any case these bugs are now fixed. The new result from odeadvisor will be [NONE] signifying that the ODE is not exact:

DEtools:-odeadvisor(ode, [exact]);

[NONE]

 and the previously wrong results from Student:-ODEs are now replaced by the expected answer

Student:-ODEs:-Solve(ode);

y(x) = c__1

With regard to the result for y(x)*diff(y(x),x)=0, there was in my opinion no bug - the result (which has not changed after this fix) was correct, although not fully simplified.

Student:-ODEs:-ODESteps(y(x)*diff(y(x),x)=0);

"[[,,"Let's solve"],[,,y(x) ((ⅆ)/(ⅆx) y(x))=0],["•",,"Highest derivative means the order of the ODE is" 1],[,,(ⅆ)/(ⅆx) y(x)],["•",,"Integrate both sides with respect to" x],[,,∫y(x) ((ⅆ)/(ⅆx) y(x)) ⅆx=∫0 ⅆx+`c__1`],["•",,"Evaluate integral"],[,,((y(x))^2)/2=`c__1`],["•",,"Solve for" y(x)],[,,{y(x)=sqrt(2) sqrt(`c__1`),y(x)=-sqrt(2) sqrt(`c__1`)}]]"

The y(x) = 0 is fully contained in the given answer when c__1 is set to 0, so the only issue really is that the two answers could be combined into one simpler answer y(x) = c__1, but this kind of polishing is not really the main goal of the Student package.

Download Answer6.mw

The issue is fixed, and the fix is distributed for everybody using Maple 2024 within the Maplesoft Physics Updates v.1747 or newer. As usual, to install the Updates, open Maple and input Physics:-Version(latest)

Austin Roche
Maplesoft

Apart from being available in the Physics updates as mentioned by @ecterrab, this fix should also appear in the next patch release (2024.1).

This bug was fixed for Maple 2024.0:

expr := -1/7 - (-1/7*7^(5/7)*exp(2/7*Pi*I)*sin(1/7*Pi)*I - 1/7*cos(1/7*Pi)*7^(5/7)*exp(2/7*Pi*I))^(7/2);
simplify(expr);

-1/7-(-((1/7)*I)*7^(5/7)*exp(((2/7)*I)*Pi)*sin((1/7)*Pi)-(1/7)*cos((1/7)*Pi)*7^(5/7)*exp(((2/7)*I)*Pi))^(7/2)

-1/7-(1/7)*(-1)^(2/7)*(-(-1)^(3/7))^(1/2)

Download mp_tmp.mw

This weakness was fixed for Maple2024.0:

expr:=c[2]*sin(sqrt(3)*x/2) + c[3]*cos(sqrt(3)*x/2);
expr:=convert(expr,tan);

simplify(expr);

1 2 Page 1 of 2