Application Center - Maplesoft

App Preview:

Enhancements in differential equations in Maple 7

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Differential Equations.mws

Maple 7 Enhancements in Differential Equations

by Waterloo Maple, 2001

Maple 7 delivers a sweeping range of new symbolic and numerical tools for solving and analyzing both ordinary and partial differential equations. It offers new algorithms for analytical solutions to PDE systems and many new classes of ODEs, especially non-linear ODEs.

Maple 7 also offers powerful new numerical routines for linear and non-linear BVPs and vastly improved speed and precision for numerical solutions of ODEs, including new implementations of the rkf45 and Rosenbrock algorithms for IVPs by Dr. Larry Shampine and Dr. Rob Corless.

These improvements are the result of years of fruitful collaboration between Waterloo Maple and universities and research institutions all over the world.

Symbolic ODEs and PDEs

Linear ODEs

The following is a an example of a 3rd order linear ODE of some complexity. In addition to demonstrating solution power, we demonstrate how to display the detailed steps Maple took in finding the solution. This feature is useful in education or research settings in which understanding the steps to solution is as important as the solution itself.

Access the suite of advanced tools required to solve and manage DE problems. Maple returns a list of available functions.

> with(DEtools);

[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...
[DEnormal, DEplot, DEplot3d, DEplot_polygon, DFacto...

Define the ODE

> myDE := x^2*diff(y(x),x,x,x)-2*(n+1)*x*diff(y(x),x,x)+(a*x^2+6*n)*diff(y(x),x)-2*a*x*y(x)=0;

myDE := x^2*diff(y(x),`$`(x,3))-2*(n+1)*x*diff(y(x)...

Increase the information level to display solution steps Maple takes in searching for a solution to the DE.

> infolevel[dsolve] := 3:

Solve the ODE. Solution steps will be diplayed. The solution will appear at the end.

> dsolve(myDE);

Methods for third order ODEs:

Trying to isolate the derivative d^3y/dx^3...

Successful isolation of d^3y/dx^3

--- Trying classification methods ---

trying a quadrature

checking if the LODE has constant coefficients

checking if the LODE is of Euler type

trying high order exact linear fully integrable

trying to convert to a linear ODE with constant coefficients

Equation is the LCLM of   -2*x/(2*(2*n-1)/a+x^2)*y(x)+diff(y(x),x), a*y(x)-2/x*n*diff(y(x),x)+diff(diff(y(x),x),x)

checking if the LODE is of Euler type

-> Attemtping a differential factorization

   trying exponential solutions

   checking if the LODE is of Euler type

   1   exponential solutions found

   exponential solutions successful

   "Partial success" case. Entering dsolve with a lower order ODE

   -> Tackling the linear ODE "as given":

      trying a quadrature

      checking if the LODE has constant coefficients

      checking if the LODE is of Euler type

      trying a symmetry of the form [xi=0, eta=F(x)]

      checking if the LODE is missing 'y' 

      trying a Liouvillian solution using Kovacic's algorithm

         No Liouvillian solutions exists

      trying a solution in terms of special functions:

         -> Bessel

         <- Bessel successful

      <- solution in terms of special functions successful

y(x) = _C1*x^(1/2+n)*BesselJ(-n-1/2,sqrt(a)*x)+_C2*...

Reset infolevel to default

> infolevel[dsolve] := 0:

Non-linear ODEs

This example solves a complicated non-linear ODE. It also illustrates the function odeadvisor() , which identifies structural elements of the equation that it will exploit to compute a solution. Finally, we show how to use the tools in the DEtools package to gain deeper insight into an ODE.

> myDE2:=diff(y(x),x,x) = -3*diff(y(x),x)^2/y(x)-2*diff(y(x),x)/x+1/(x^6*diff(y(x),x)*y(x)^6);

myDE2 := diff(y(x),`$`(x,2)) = -3*diff(y(x),x)^2/y(...

Ask Maple for a structural analysis of the non-linear ODE.

> odeadvisor( myDE2 );

[[_2nd_order, _with_linear_symmetries], [_2nd_order...

The advisor's analysis tells us that (1) the ODE is 2nd order (which is obvious), (2) the ODE has linear symmetries (not as obvious), (3) the ODE is reducible to first order, and (4) Maple is able to compute an integrating factor polynomial in y'. Maple can exploit all of this information in its solution process, a unique feature among DE solvers in computer algebra systems.

But let's first try solving the differential equation directly, pretending we don't know anything about the ODE's structure. The default method returns all 8 solutions explicitly as functions of x .

> dsolve(myDE2);

y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...

How did Maple exploit the structural information to obtain this answer? We can examine Maple's thinking more deeply using the DEtools package.

First let's use the symgen function to ask exactly what "linear symmetries" were found.

> mySymmetryInfo := DEtools[symgen]( myDE2 );

mySymmetryInfo := [_xi = x^2, _eta = 0], [_xi = -8/...

Now tell Maple to resolve the ODE using that symmetry information as a "hint".

> dsolve( myDE2, HINT = [ mySymmetryInfo ] );

-1/x-1/8*2^(2/3)*(-3*y(x)^4+_C1)^(2/3)-_C2 = 0, -1/...
-1/x-1/8*2^(2/3)*(-3*y(x)^4+_C1)^(2/3)-_C2 = 0, -1/...

Second, we might like to know what integrating factor was found. Let's find out, using the intfactor function.

> myIntegratingFactor := DEtools[ intfactor ]( myDE2 );

myIntegratingFactor := diff(y(x),x)*x^4*y(x)^6

If you were teaching a course on DEs, you could ask students to verify the integrating factor. That is, form a new DE by multiplying the original by the integrating factor and trying to reduce the order of the resulting DE.

> newDE := myIntegratingFactor * myDE2;

newDE := diff(y(x),x)*x^4*y(x)^6*diff(y(x),`$`(x,2)...

If the integrating factor is valid, then the new DE will be exact, and thus reducible from 2nd order to 1st order. (Recall the odeadvisor told us the original DE was reducible!) This reduction is accomplished using the firint function.

> reducedDE := DEtools[ firint ]( newDE );

reducedDE := 1/2*diff(y(x),x)^2*x^4*y(x)^6+1/x+_C1 ...

Finally, let's solve this 1st order DE and verify that its solution is identical to that of the original.

> dsolve( reducedDE );

y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...
y(x) = 1/3*3^(3/4)/x*((-8*sqrt(-2*x-2*_C1*x^2)-8*x*...

We can shorten the output by re-solving the DE with the implicit option. Notice the output below expresses y(x) as two implicit functions of x instead of eight explicit functions.

> dsolve( reducedDE, implicit);

y(x)^4-4/3*1/x^3*(-2*x-2*_C1*x^2)^(3/2)-_C2 = 0, y(...

Systems of PDEs

Maple 7 is the only commercial system that can solve systems of partial differential equations symbolically. The following example illustrates the solution of 3 simultaneous PDEs.

Define the PDE system

> myPDEsystem := [-y*diff(f(x,y,z,t),x) + z^2*diff(f(x,y,z,t),z) + 3*t*z*diff(f(x,y,z,t),t) - 3*t^2-4*f(x,y,z,t)*z = 0,
-y*diff(f(x,y,z,t),y) - z*diff(f(x,y,z,t),z) - t*diff(f(x,y,z,t),t) + f(x,y,z,t) = 0,
-x*diff(f(x,y,z,t),y) - diff(f(x,y,z,t),z) = 0]:

for _eq in myPDEsystem do
_eq;
od;

-y*diff(f(x,y,z,t),x)+z^2*diff(f(x,y,z,t),z)+3*t*z*...

-y*diff(f(x,y,z,t),y)-z*diff(f(x,y,z,t),z)-t*diff(f...

-x*diff(f(x,y,z,t),y)-diff(f(x,y,z,t),z) = 0

Compute the solution.

> sol := pdsolve(myPDEsystem);

sol := {f(x,y,z,t) = (_C1*(x*z-y)*t^(3/2)+3*t^2*x*s...

Numerical DEs

Boundary Value Problems (BVPs)

Boundary value problems (BVPs) for ODEs constrain the solution to specific values at the end points of the domain and in many ways are more difficult than ODE initial value problems (IVPs), in which a variable (or one or more derivatives) are specified when the independent variable (e.g. t ) is 0.

Maple 7 provides powerful new algorithms for solving both linear and non-linear BVPs. Maple 7 now solves more classes of numerical BVPs and with greater accuracy than any of our competitors.

This example illustrates Maple 7's ability to solve a non-linear 2x2 BVP system, first numerically, then symbolically. This BVP system comes from the Maxwell-Stefan equations, which are used to model multicomponent diffusional fluxes in chemical engineering.

Numerical Solution

> de1 := diff((1+x(t)/(1-x(t)-y(t)))*diff(x(t),t) +
x(t)/(1-x(t)-y(t))*diff(y(t),t), t);

de1 := (diff(x(t),t)/(1-x(t)-y(t))-x(t)/(1-x(t)-y(t...
de1 := (diff(x(t),t)/(1-x(t)-y(t))-x(t)/(1-x(t)-y(t...

> de2 := diff((y(t)/(1-x(t)-y(t)))*diff(x(t),t) +
(1+y(t)/(1-x(t)-y(t)))*diff(y(t),t), t);

de2 := diff(x(t),t)/(1-x(t)-y(t))*diff(y(t),t)-y(t)...
de2 := diff(x(t),t)/(1-x(t)-y(t))*diff(y(t),t)-y(t)...

> bvpSolution := dsolve( {de1=0, de2=0,
x(0)=.418, y(0)=.573, x(1)=0, y(1)=0},
numeric):

> bvpSolution( .3 );

[t = .3, x(t) = .406198039810033996, diff(x(t),t) =...

> plots[odeplot](bvpSolution, [[t,x(t)],[t,y(t)]],0..1,numpoints=25, scaling=constrained, thickness=2);

[Maple Plot]

Symbolic Solution

> bvpSolution := simplify( dsolve( {de1=0, de2=0, x(0)=.418, y(0)=.573, x(1)=0, y(1)=0}) );

bvpSolution := {y(t) = 573/991000*(1000*exp(-(3*ln(...
bvpSolution := {y(t) = 573/991000*(1000*exp(-(3*ln(...

> assign(bvpSolution);

> plot([x(t), y(t)], t=0..1, color=[red, green], scaling=constrained, thickness=2);

[Maple Plot]

Initial Value Problems (IVPs)

Maple 7 has vastly improved both the speed and accuracy of its IVP solvers, especially the rkf45 and dverk78 methods. Both offer 100% speed ups over their implementations in Maple 6, and both, for the first time, can compute solutions to arbitrary accuracy.

The improvement in the rkf45 method is due to the work of Larry Shampine and Rob Corless , renowned researchers in the field of DEs.

We use the Lotka-Volterra Predator-Prey model as our sample problem:

> restart;

> predatorPreyDE := { diff( x(t),t) = x(t)*(1-y(t)),
diff( y(t),t) = 0.3*y(t)*(x(t)-1),
x(0) = 0.5, y(0) = 1 };

predatorPreyDE := {diff(x(t),t) = x(t)*(1-y(t)), di...

100% speed up and arbitrary accuracy for the rkf45 method. Notice the option to specify the error tolerance, in this example set to 1e-8.

> populations := dsolve( predatorPreyDE, numeric,
method = rkf45,
abserr = 1e-8, relerr = 1e-8):

> populations( 1000 );

[t = 1000., x(t) = 1.29348611666789614, y(t) = .723...

> plots[odeplot]( populations, 0..100, thickness=2, numpoints=50, color=blue );

[Maple Plot]

200% speed up and arbitrary accuracy for the dverk78 method. Notice the option to specify the error tolerance, in this example set to 1e-8.

> populations := dsolve( predatorPreyDE, numeric,
method=dverk78,
abserr=1e-8, relerr=1e-8);

> populations(1000);

[t = 1000., x(t) = 1.29351555619824564, y(t) = .723...