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);
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;
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
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);
Ask Maple for a structural analysis of the non-linear ODE.
>
odeadvisor( myDE2 );
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);
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 );
Now tell Maple to resolve the ODE using that symmetry information as a "hint".
>
dsolve( myDE2, HINT = [ mySymmetryInfo ] );
Second, we might like to know what integrating factor was found. Let's find out, using the
intfactor
function.
>
myIntegratingFactor := DEtools[ intfactor ]( myDE2 );
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;
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 );
Finally, let's solve this 1st order DE and verify that its solution is identical to that of the original.
>
dsolve( reducedDE );
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);
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;
Compute the solution.
>
sol := pdsolve(myPDEsystem);
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);
>
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);
>
bvpSolution := dsolve( {de1=0, de2=0,
x(0)=.418, y(0)=.573, x(1)=0, y(1)=0},
numeric):
>
bvpSolution( .3 );
>
plots[odeplot](bvpSolution, [[t,x(t)],[t,y(t)]],0..1,numpoints=25, scaling=constrained, thickness=2);
Symbolic Solution
>
bvpSolution := simplify( dsolve( {de1=0, de2=0, x(0)=.418, y(0)=.573, x(1)=0, y(1)=0}) );
>
assign(bvpSolution);
>
plot([x(t), y(t)], t=0..1, color=[red, green], scaling=constrained, thickness=2);
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 };
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 );
>
plots[odeplot]( populations, 0..100, thickness=2, numpoints=50, color=blue );
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);