How Do I
Solve an Ordinary Differential Equation?
This topic introduces you to the commands and techniques used to solve ordinary differential equations (ODEs) in Maple. Unless noted otherwise, each example goes through the following steps:
Defining your system of ODEs.
Finding solutions.
Assigning the solutions to names (thus making them easier to plot and investigate).
Plotting your solution.
Before you start looking at the examples on this page, you should consider what type of solution you need: symbolic or numeric.
Symbolic solutions return equations that contains your independent variables. These solutions are exact and can be easily manipulated to find, for example, a value at a point, a derivative of the solution, or an integral of the solution. Note that not every ODE or system of ODEs has a symbolic solution. However, if you are taking an introductory course in ODEs, it is most likely that you'll be interested in finding symbolic solutions.
Numeric solutions return a procedure that integrates your ODEs. Since they are integrated numerically, numeric solutions are approximations. Most ODEs can be solved numerically even if they don't have a symbolic solution. If you need to find derivatives or integrals of your solution, you should include equations for the derivatives and integrals in your system of ODEs.
This topic contains separate sections on how to find symbolic and numeric solutions since the techniques for solving those problems differ slightly. The examples in each section are arranged from simple to more complex.
Finally, if you do not find the information you are looking for in this topic, the Additional Resources section contains links to other pages that go into more detail about solving ODEs.
Solving an ODE Symbolically
Solving an ODE Numerically
Additional Resources
Related Topics
The default behavior of dsolve is to find a symbolic solution for your ODE. After you have a symbolic solution, you can assign it to a variable so that you can evaluate and graph it.
Basic Example: y' =a y
Taking Derivatives and Integrals of Symbolic Solutions
Solving a System of ODEs
Solving a Linear Two-Point BVP for a Second-order ODE
Implicit Solutions with RootOf
Basic Example: y' = a y
The following example shows how to find a symbolic solution to the ODE ⅆⅆ x y=a⋅y.
Define a variable called basicODE and set it equal to the ODE.
restart;
basicODE≔diffyx,x=a⋅yx;
basicODE≔ⅆⅆxyx=ayx
Use dsolve(basicODE) to find the solution. Here, we assign the result to a variable called basicSolution. The result is an equation that satisfies the ODE. Terms of the form _Cn (where n is an integer) are constants of integration.
basicSolution≔dsolvebasicODE;
basicSolution≔yx=_C1ⅇax
To find the constant for a particular solution, include an initial value equation with the ODE in a set or list and then pass the set/list to dsolve. The following expression finds a solution that satisfies the condition y=5 when x=0.
basicODE≔diffyx,x=a⋅yx, y0=5;
basicODE≔ⅆⅆxyx=ayx,y0=5
basicSolution≔yx=5ⅇax
Test your solution with the odetest function. The odetest function takes the symbolic solution from dsolve and your ODEs as parameters. It returns 0 if the solution satisfies the ODEs.
odetestbasicSolution, basicODE;
0
To find the value of y at a particular value of x or to graph the solution, you can use either the rhs function or the eval(expression, equations) command to assign the right-hand side of the solution to a variable. Of the two methods, eval is more flexible (especially when dealing with a set of solutions), and will be used in the rest of the examples in this section.
Y≔ rhsbasicSolution;Y≔evalyx,basicSolution;
Y≔5ⅇax
After assigning a value to the parameter a, use eval on Y to get an exact value for the solution at a point.
a≔−13;
a≔−13
evalY, x=5;
5ⅇ−53
Use the evalf function to get a floating-point value.
evalfevalY, x=5;
0.9443780140
To visualize the solution, use the plot function. The following graphs Y over x=0..5.
plotY,x=0..5, labels=x, y(x);
Finally, the function Y can be used in fsolve to find, for example, the value of x for which y(x)=2.
b≔fsolveYx=2,x;
b≔2.748872196
evalY,x=b;
2.000000000
Symbolic solutions can be differentiated and integrated using the commands diff and int.
In this example, we have a first-order differential equation in v(t), the velocity of an object. Once we find v(t), we want to take its derivative to find the object's acceleration, and then take the integral of v(t) to find the distance traveled by the object. Finally, we want to plot the velocity, acceleration, and distance on a single plot.
We start by defining equations, a set that contains the ODE and the initial condition (we assume the initial velocity to be 0).
equations≔diffvt,t=−1+vt2, v0=0
equations≔ⅆⅆtvt=−1+vt2,v0=0
Run dsolve(equations) to get a solution for v(t).
solution≔dsolveequations;
solution≔vt=−tanht
odetestsolution,equations;
Assign the solution for v(t) to a variable, V:
V≔evalvt,solution;
V≔−tanht
You can use the diff function to take the derivative of V and get the acceleration of the object, A.
A≔diffV,t;
A≔−1+tanht2
Similarly, you can use the int command to integrate V and get S, the distance traveled by the object. Since V is an expression in t, we need to integrate V over t from our initial time, 0, to some final time, T, giving us an expression in T. (Note that we assume the initial position of the object to be 0.)
S≔intV,t=0..T;
S≔−lncoshT
If you are only interested in solutions for T > 0, you can specify this by using assuming on the int expression.
S≔intV,t=0..T assuming T>0;
S≔ln2−lnⅇ2T+1+T
Note that this expression is equal to our initial expression for S, -ln(cosh(T)).
To get an expression in t, evaluate S at t:
S≔evalS,T=t;
S≔ln2−lnⅇ2t+1+t
The following command plots our solution, V, its derivative, A, and its integral, S.
plotV,A,S,t=0..3,legend=V(t) - velocity,A(t) - acceleration, S(t) - distance;
When solving a system of ODEs, enclose the system in a set or list. The following is an example of a system with two ODEs enclosed in a set.
odeSystem≔ diffxt,t=−xt− yt, diffyt,t= xt−yt
odeSystem≔ⅆⅆtxt=−xt−yt,ⅆⅆtyt=xt−yt
Run dsolve(odeSystem) to get a solution. The result is a set of equations that satisfies the system.
systemSol≔dsolveodeSystem;
systemSol≔xt=ⅇ−t_C2cost+_C1sint,yt=−ⅇ−tcost_C1−sint_C2
odetestsystemSol,odeSystem;
To access the equations in the solution, use the eval(expression, equations) command. Enter either x(t) or y(t) for expression; enter systemSol for equations.
X≔ evalxt,systemSol;
X≔ⅇ−t_C2cost+_C1sint
Y≔evalyt,systemSol;
Y≔−ⅇ−tcost_C1−sint_C2
Include the initial conditions in the set of ODEs to find a particular solution for the system. The following specifies that the solution should go through the point (1,2) at t=0.
odeSystem≔ diffxt,t=−xt− yt, diffyt,t= xt−yt, x0=1, y0=2
odeSystem≔ⅆⅆtxt=−xt−yt,ⅆⅆtyt=xt−yt,x0=1,y0=2
systemSol≔xt=ⅇ−tcost−2sint,yt=−ⅇ−t−2cost−sint
X≔evalxt,systemSol;
X≔ⅇ−tcost−2sint
Y≔−ⅇ−t−2cost−sint
Use eval and evalf to analyze your solution.
evalX, t=Pi; evalfevalX, t=Pi;
−ⅇ−π
−0.04321391825
evalY, t=Pi; evalfevalY, t=Pi;
−2ⅇ−π
−0.08642783650
The following command finds the first time x(t) = y(t) with t ≥ 1.5.
fsolveXt=Yt,t=1.5;
2.819842099
The X and Y variables can be used to graph the solutions using the plot function. The following is a graph of x(t) and y(t) versus t.
plotX, Y, t=0..5, legend=x(t), y(t);
The following is a parametric plot of the solution over the same time interval as the preceding graph.
plotX,Y,t=0..5,labels=x,y;
Solving a Linear Two-Point BVP for a Second-Order ODE
Finding solutions to a two-point boundary value problem (BVP) is more involved than solving an initial value problem. For example, let's say you have an ODE that can always be solved as an initial value problem. When that same ODE is cast as a two-point BVP, there is no guarantee that a solution can be found. In fact, we will see that to a two-point BVP can result in any of the following:
A unique solution
No solution
An infinite number of solutions
Note: The theory that explains the behavior of the two-point BVP for linear second-order ODEs is attributed to Sturm and Liouville. Such problems are called Sturm-Liouville problems.
Let's start off by finding the solution for an ODE given some initial conditions. After that, we'll look at solving the same ODE but as a two-point BVP.
ode≔diffyx,x,x+yx=0;
ode≔ⅆ2ⅆx2yx+yx=0
sol≔dsolveode;
sol≔yx=_C1sinx+_C2cosx
The general solution has two constants, _C1 and _C2. To get a specific solution, include two initial value equations along with ode in dsolve.
iv1≔y0=y1; iv2≔y'0=y2;
iv1≔y0=y1
iv2≔Dy0=y2
sol≔dsolveode,iv1,iv2;odetestsol,ode,iv1,iv2;
sol≔yx=y2sinx+y1cosx
This verifies that we can always find a solution for this ODE given the initial values we defined.
Let's now look at the solutions we get from the same ODE as part of a two-point BVP. First, define two generic boundary value conditions at x = 0 and x = π. Our boundary values are linear combinations of y(x) and y'(x) at the boundaries.
bc1≔a⋅y0+b⋅y'0=f;
bc1≔ay0+bDy0=f
bc2≔c⋅yPi+d⋅y'Pi=g;
bc2≔cyπ+dDyπ=g
When we solve the system consisting of ode, bc1, and bc2, we get the following. (Note that the solution exists and is unique as long as a⋅d - b⋅c ≠ 0.)
solution≔dsolveode,bc1,bc2;
solution≔yx=−ag+cfsinxad−bc+bg+dfcosxad−bc
odetestsolution,ode,bc1,bc2;
Y≔evalyx,solution;
Y≔−ag+cfsinxad−bc+bg+dfcosxad−bc
To analyze the solution, assign values to the parameters a, b, c, d, f, and g (ensuring, of course, that a⋅d - b⋅c ≠ 0).
a≔1; b≔2; c≔3; d≔4; f≔34; g≔−2; bc1; bc2;
a≔1
b≔2
c≔3
d≔4
f≔34
g≔−2
y0+2Dy0=34
3yπ+4Dyπ=−2
Y;
sinx8+cosx2
plotY,x=0..Pi
We have found the case where we get a unique solution to the system. To see the other two possibilities (that is, no solution and an infinite number of solutions), let's examine what happens when we set a, b, c, and d such that a⋅d - b⋅c = 0. We do this by setting the boundary values to y(0) = 1 and y(π) = 0.
a≔1; b≔0; c≔1; d≔0; f≔1; g≔0; bc1; bc2;
b≔0
c≔1
d≔0
f≔1
g≔0
y0=1
yπ=0
sol≔dsolveode,bc1,bc2;
sol≔
We get no solution from dsolve. To understand why, recall that the general solution to our ODE is:
y(x) = _C1⋅sin(x) + _C2⋅cos(x)
When we apply our boundary values to the general solution, we get the following equations for the constants _C1 and _C2.
_C1⋅sin(0) + _C2⋅cos(0) = 1 ⇒ _C2 = 1 _C1⋅sin(π) + _C2⋅cos(π) = 1 ⇒ _C2 = 0
Our boundary values, therefore, do not allow for a solution to our system.
Finally, let's see what happens when we set the boundary values to y(0) = 1 and y(π) = -1.
g≔−1; bc1; bc2;
g≔−1
yπ=−1
sol≔dsolveode, bc1, bc2;
sol≔yx=_C1sinx+cosx
odetestsol,ode,bc1,bc2;
We get a solution from dsolve, but it does not look complete. One of the constants in the general solution was found, but the other, _C1, remains in the solution. We therefore have infinitely many solutions to this BVP since any multiple of sin(x) can be added to cos(x). To understand why this happens, apply the boundary values to the general solution to get the following equations.
_C1⋅sin(0) + _C2⋅cos(0) = 1 ⇒ _C2 = 1 _C1⋅sin(π) + _C2⋅cos(π) = -1 ⇒ _C2 = 1
We can see that our boundary values enable us to solve for _C2, but give no criteria at all for _C1. This means that _C1 can take on any value, which is consistent with the solution found by dsolve.
The dsolve function always tries to return an explicit solution. If the solution cannot be expressed explicitly, dsolve returns a solution that contains a RootOf expression.
Consider the following ODE and initial condition.
odesys≔diffyx,x=1cosyx+1.01, y0=1;
odesys≔ⅆⅆxyx=1cosyx+1.01,y0=1
sol≔dsolveodesys;
sol≔yx=RootOf101_Z−100x+100sin_Z−101−100sin1
The RootOf expression in the solution for this ODE tells us that an explicit solution for the ODE could not be found. (Note that odetest verifies that the solution does satisfy the ODE.)
odetestsol, odesys;
There are two ways to solve this. The first is to call dsolve with the implicit option to force dsolve to return an implicit equation for y(x).
sol1≔dsolveodesys,implicit;
sol1≔x−101yx100−sinyx+101100+sin1=0
The second is to use the DEtools:-remove_RootOf command.
withDEtools:
sol2≔remove_RootOfsol;
sol2≔101yx−100x+100sinyx−101−100sin1=0
To graph the solution, use the plots:-implicitplot command.
withplots:
implicitplotsol2,x=−10..10, y=−10..15;
For this example we can use eval and fsolve to find the value for x given a value of y(x).
eq≔evalsol2, yx=8;
eq≔707−100x+100sin8−100sin1=0
fsolveeq,x;
7.217887262
Use the numeric parameter with dsolve to get a numerical solution for your ODE or ODE system. The numerical solution is a procedure that you execute to analyze.
Basic Example: y' = -y
Basic Example with Parameters: y' = a y
Solving a System of ODEs Numerically
Taking Derivatives and Integrals of Numeric Solutions
Techniques for Solving a Nonlinear Two-Point BVP
Solving Piecewise ODEs and an Introduction to Event Handling
Finding Roots Using Events and the Halt Action
An Advanced Example with Piecewise ODEs and Event Detection
This example shows how to use dsolve get a numerical solution to ⅆⅆ x y= −y.
Define a variable called basicOde, a set that contains the ODE and the initial condition, y(0)=5.
basicOde≔diffyx,x = −yx, y0=5;
basicOde≔ⅆⅆxyx=−yx,y0=5
To get a numerical solution, use dsolve with the numeric parameter.
sol≔dsolvebasicOde,numeric;
sol≔procx_rkf45...end proc
The output is a procedure that you run to analyze the solution to your ODE. In the output expression, the x_rkf45 term tells you the numerical method used in the solution (in this case rkf45, a Fehlberg fourth-fifth order Runge-Kutta).
To get a value for your solution at a point, run the procedure with the point as the parameter. The following command finds the value for y when x=5.0. The output is a list of equations. The first element in the list echoes your input. The second element gives you the value of y(x).
sol5.0;
x=5.0,yx=0.0336896836290933
sol0.916;
x=0.916,yx=2.00058112757547
Use the plots:-odeplot command to graph your solution. The odeplot command accepts the solution from dsolve/numeric as its first parameter. The following commands plot sol over x=0..5.
odeplotsol,x=0..5;
Note that this solution is very specific in terms of the initial conditions and the terms in the ODE. The next section shows how you can generalize the solution by using the parameters option with dsolve/numeric.
Assigning the Solution to a Variable
The procedure that dsolve returns is useful for finding the value of the solution at a point and plotting the solution with odeplot, but it cannot be used in fsolve. For example, let's say you want to find the value of x for which y(x)=2. Using sol directly in fsolve returns the same expression unevaluated.
fsolvesolx=2,x;
fsolvesolx=2,x
To solve this, use the output=listprocedure option with dsolve, and then assign the solution to a variable with the eval(expression,equations) command.
sol≔dsolvebasicOde,numeric,output=listprocedure;
sol≔x=procx...end proc,yx=procx...end proc
Y≔evalyx,sol;
Y≔procx...end proc
a≔fsolveYx=2,x;
a≔0.9162905219
The following commands verify that the output from sol and Y are identical (and that the result from fsolve is correct).
sola;
x0.9162905219=0.9162905219,yx0.9162905219=1.99999999996477
Ya;
1.99999999996477
With the solution assigned to a variable, you can use the standard plot command to graph your solution instead of plots:-odeplot.
plotY,0..5;
This example shows how to use dsolve get a numerical solution to ⅆⅆ x y=a⋅y with the initial condition y(x0) = y0. The terms a, x0, and y0 are parameters that you will be able to set after the solution has been generated.
Define a variable called basicOde2, a set that contains the ODE and the initial condition, y(x0) = y0.
basicOde2≔diffyx,x = a⋅yx, yx0=y0;
basicOde2≔ⅆⅆxyx=ayx,yx0=y0
Use dsolve to get a numerical solution. Also, because basicOde2 contains the parameters a, x0, and y0,use the parameters option.
sol2≔dsolvebasicOde2,numeric, parameters=a, x0, y0;
sol2≔procx_rkf45...end proc
Before you can work with the solution, you must initialize the parameters a, x0, and y0. To initialize the parameters, run the output procedure with values inserted in the parameters list. The following sets a to -1, x0 to 0, and y0 to 5.
sol2parameters=−1,0,5;
a=−1.,x0=0.,y0=5.
To get a value for your solution at a point, run the procedure with the point given as the argument. The following command finds the value for y when x=5.0. The output is a list of equations. The first element in the list echoes your input. The second element gives you the value of y(x).
sol25.0;
The following graphs sol2 over x=0..5.
odeplotsol2,x=0..5;
At any time you can change the values for the parameters to get a different solution for your system.
sol2parameters=−13,0,5;
a=−0.333333333333333,x0=0.,y0=5.
Finally, if you want to analyze the solution using fsolve, you must generate a solution using the output=listprocedure option, and then assign the solution to a variable.
sol2≔dsolvebasicOde2,numeric, parameters=a, x0, y0,output=listprocedure;
sol2≔x=procx...end proc,yx=procx...end proc
Y≔evalyx,sol2;
Before you can use the function, you must set values for the parameters a, x0, and y0.
Yparameters=−1,0,5;
You can now use fsolve on the solution. The following finds the value of x for which the solution equals 2.
b≔0.9162905219
Yb;
Define the system of ODEs along with their initial conditions. The system of ODEs in this example is the Lotka-Volterra predator-prey system, where x(t) represents the population of the prey and y(t) the population of the predator.
des≔diffxt,t=a ⋅xt−b ⋅xt⋅yt,diffyt,t=−c ⋅yt+d⋅ xt⋅yt, xt0=x0, yt0=y0;
des≔ⅆⅆtxt=axt−bxtyt,ⅆⅆtyt=−cyt+dxtyt,xt0=x0,yt0=y0
Generate the solution with a list of parameters.
desSol≔dsolvedes,numeric,parameters=a,b,c,d,t0,x0,y0;
desSol≔procx_rkf45...end proc
Before working with the solution, set values for the parameters.
desSolparameters=0.1,0.01,0.025,0.001,0,50,15;
a=0.1,b=0.01,c=0.025,d=0.001,t0=0.,x0=50.,y0=15.
To find the populations at a particular time, run the solution with the time given as the argument. The result is a list of equations for the independent variable (t), the prey population (x(t)), and the predator population (y(t)).
desSol90;
t=90.,xt=20.0599816840952,yt=5.29175138992018
To graph multiple curves with odeplot, use a nested list that contains a list for each pair of coordinates you would like to graph. For example, for the following odeplot command graph, both the prey and predator populations are displayed on the same graph.
odeplotdesSol,t,xt, t, yt, 0..300, labels=t,Population,legend=x(t) - Prey, y(t) - Predator;
If three terms are given for the variables list, a 3-D graph is produced.
odeplotdesSol,t,xt,yt,0..300;
Finally, a parametric graph of the solution, with the prey population on the horizontal axis and the predator population on the vertical axis.
odeplotdesSol, xt,yt, 0..150,labels=x - Prey,y - Predator;
Assigning Solutions to Variables
If you want to analyze your solution using fsolve or other commands, use the output=listprocedure option with dsolve and then assign the solutions to variables.
desSol≔dsolvedes,numeric,parameters=a,b,c,d,t0,x0,y0,output=listprocedure;
desSol≔t=proct...end proc,xt=proct...end proc,yt=proct...end proc
X≔evalxt,desSol;
X≔proct...end proc
Y≔evalyt,desSol;
Y≔proct...end proc
Assigning parameter values for either X or Y automatically assigns the parameter values for the other variable.
Xparameters=0.1,0.01,0.025,0.001,0,50,15;
Yparameters;
plotXt,Yt,t=0..150,labels=t,Population,legend=x(t) - Prey,y(t) - Predator;
The following command uses fsolve to find the first time the populations x(t) and y(t) are equal.
a≔fsolveXt=Yt,t;
a≔16.82819661
Xa;Ya;
16.5912730797044
16.5912730836271
To find the derivative or integral of a numeric solution, you should include ODEs for the derivative or integral in your system of equations.
For this example, we have a first-order differential equation in v(t), the velocity of an object, along with an initial condition.
velocityEqs≔diffvt,t=−1+vt2,v0=v0;
velocityEqs≔ⅆⅆtvt=−1+vt2,v0=v0
Let's suppose that in addition to the velocity we want to find the object's acceleration, a(t), and position, s(t). We could use the preceding equations to find v(t), and then take the derivative and integral of v(t) to find the acceleration and position. However, taking the derivative or integral of a numeric solution is not recommended (see Why You Should Not Use Int or Diff on a Numeric Solution).
Instead, you should define an augmented system of ODEs that includes the ODEs and initial conditions for the object's velocity, acceleration, and position. In the following command, we define equations, an augmented system that contains the equations for the acceleration and position.
equations≔velocityEqs, at=diffvt,t, diffst,t=vt, s0=s0
equations≔at=ⅆⅆtvt,ⅆⅆtst=vt,ⅆⅆtvt=−1+vt2,s0=s0,v0=v0
Generate the solution with parameters for v0 and s0.
solution≔dsolveequations,numeric,parameters=v0,s0,output=listprocedure;
solution≔t=proct...end proc,at=proct...end proc,st=proct...end proc,vt=proct...end proc
Set the initial velocity and initial position both to zero.
solutionparameters=0,0;
tparameters=0,0=v0=0.,s0=0.,atparameters=0,0=v0=0.,s0=0.,stparameters=0,0=v0=0.,s0=0.,vtparameters=0,0=v0=0.,s0=0.
Assign the velocity, acceleration, and position to the variables V, A, and S, respectively.
V≔proct...end proc
A≔evalat,solution;
A≔proct...end proc
S≔evalst,solution;
S≔proct...end proc
Finally, use plot to plot the velocity, acceleration, and position.
plotVt,At,St,t=0..3,legend=V(t) - velocity,A(t) - acceleration,S(t) - position;
Why You Should Not Use Int or Diff on a Numeric Solution
Although it is possible to use the Int or Diff commands to numerically integrate or differentiate the procedure from a numeric solution, doing so will lead to results that will either take much longer to calculate, be inaccurate, or both. The following example illustrates this.
First, we'll define a second set of equations with only the velocity, and then we'll assign the solution for the velocity to V2.
equations2≔diffvt,t=−1+vt2, v0=0;
equations2≔ⅆⅆtvt=−1+vt2,v0=0
sol2≔dsolveequations2,numeric,output=listprocedure;
sol2≔t=proct...end proc,vt=proct...end proc
V2≔evalvt,sol2;
V2≔proct...end proc
Next, we'll define SInt by integrating our numerical solution, V2.
SInt≔T→evalfIntV2t,t=0..T;
SInt≔T↦evalf∫0TV2tⅆt
The following commands compare the results from SInt with S and the exact solution, -ln(cosh(T)) (which was found in the example on taking integrals and derivatives of symbolic solutions).
CodeTools:-UsageSInt6;CodeTools:-UsageS6;evalf−lncosh6;
memory used=18.88MiB, alloc change=36.00MiB, cpu time=296.00ms, real time=728.00ms, gc time=31.20ms
−5.30685989107477
memory used=8.88KiB, alloc change=0 bytes, cpu time=0ns, real time=2.00ms, gc time=0ns
−5.30685896501487
−5.306858964
Notice that SInt is takes longer to calculate the integral and is also less accurate when compared to the exact solution.
We will now use the Diff function to find ADiff, the derivative of V2, and compare that with A and the exact solution, tanh(t)2−1.
ADiff≔T→evalDiffV2t,t,t=T;
ADiff≔T↦ⅆⅆtV2tt=T|ⅆⅆtV2tt=T
The following commands compare ADiff with A and the exact solution, tanh(t)2−1.
CodeTools:-UsageevalfADiff6; CodeTools:-UsageA6; evalftanh62−1;
memory used=110.88KiB, alloc change=0 bytes, cpu time=0ns, real time=27.00ms, gc time=0ns
−0.00002442061000
memory used=7.10KiB, alloc change=0 bytes, cpu time=0ns, real time=0ns, gc time=0ns
−0.0000245764945147024
−0.0000245764
The result from ADiff is less accurate than the result from A.
This should convince you that it is always better to augment your system of ODEs with equations for the integral or derivative of your solution.
This example looks at a nonlinear two-point boundary value problem (BVP). The system consists of the following nonlinear ODE and two nonlinear boundary conditions (bc1 and bc2).
ode ≔ diffyx,x,x+yx2⋅diffyx,x+cosyx=x;
ode≔ⅆ2ⅆx2yx+yx2ⅆⅆxyx+cosyx=x
bc1≔1+y0⋅y'0−3⋅y0=0;
bc1≔1+y0Dy0−3y0=0
bc2≔ 1+y1⋅y'1−3⋅y1−y12+0.5=0;
bc2≔1+y1Dy1−3y1−y12+0.5=0
We want to find all of the solutions to this system in the interval 0≤x≤1.
A first attempt at using dsolve/numeric fails.
solution≔dsolveode,bc1,bc2,numeric;
Error, (in dsolve/numeric/bvp) unable to store '0.562014431217500e-1-0.175146559896416e-14*I' when datatype=float[8]
Despite this, we will show that this system can be solved and we will give two techniques for finding the solutions.
The first technique finds a solution by giving dsolve an approximate solution (that is, by using the approxsoln option).
The second technique finds all of the solutions for the nonlinear BVP by turning the BVP into an initial value problem.
Getting One Solution Using approxsoln
For some nonlinear boundary conditions it is necessary to give dsolve an initial approximation to the solution. To provide the initial approximation, use the approxsoln option. This forces dsolve to try finding a solution that is close to the initial approximation.
For our system, let's try using yx=x as our initial approximation.
solution≔dsolveode,bc1,bc2,numeric, approxsoln=yx=x;
solution≔procx_bvp...end proc
plots:-odeplotsolution;
We succeeded in finding a solution, but we do not know if this is the only solution or if there are more solutions to our system. The next section demonstrates how to find all of the solutions that exist for the system.
Finding all possible solutions
One option to find other solutions for our nonlinear two-point BVP is to try different initial approximations. This, however, would be tedious and would not guarantee that every solution is found.
A better approach is to look at the boundary conditions we have and use them to turn our nonlinear two-point BVP into an initial value problem and from there try finding all of the solutions to the ODE. That is, we will implement the shooting method.
Our boundary conditions consist of two equations with four unknowns (y(0), y'(0), y(1), and y'(1)):
1+y0⋅y'0−3⋅y0=0
1+y1⋅y'1−3⋅y1−y12+0.5=0
We want to rewrite these equations in terms of two unknowns, y(0) and y'(0), so that we can solve for these unknowns and then use them as the initial values for our ODE.
Start by making the substitutions y(0) = a and y'(0) = b in the first equation:
1+a⋅b−3⋅a=0
Next, define a procedure, h( a, b), that, given the initial values a and b, numerically integrates our ODE and returns the endpoint values, y(1) and y'(1), in a list.
h≔proca,b local F; Digits ≔10; if typea,numeric and typeb,numeric then # Find a solution to the ODE, but now as an initial value # problem with y(0)=a and y'(0)=b. F≔dsolveode,y0=a,Dy0=b,yx,numeric; # Return a list with two values, [y(1), y'(1)]. rhsF12,rhsF13; else 'h'a,b; end if; end proc:
We can now substitute h(a, b)[1] for y(1) and h(a, b)[2] for y'(1) in our second equation. This gives us two equations with two unknowns.
1+ha,b1⋅ha,b2−3⋅ha,b1−ha,b12+0.5=0
Before we can solve for the (a, b) pairs, we need to redefine our boundary conditions as procedures. This makes it easier to analyze our boundary conditions with the implicitplot and fsolve commands.
BC1≔proca,b 1+a⋅b−3⋅a end proc;
BC1≔proca,b1+a*b − 3*aend proc
BC2≔proca,b 1+ha,b1⋅ha,b2−3⋅ha,b1−ha,b12+0.5 end proc;
BC2≔proca,b1+ha,b[1]*ha,b[2] − 3*ha,b[1] − ha,b[1]^2+0.5end proc
In the following implicit graph of BC1(a, b) = 0 and BC2(a, b) = 0, the points where the curves intersect are the (a, b) pairs that solve our BVP. There are four solutions to our BVP. Note that the solution we found using approxsoln should correspond to the point of intersection in the first quadrant (since that is the only solution with a > 0).
implicitplotBC1a,b=0,BC2a,b=0, a=−5..1,b=−11..7,color=black,red,legend=BC1(a,b)=0,BC2(a,b)=0;
To find the values for (a, b), use the fsolve command with an interval around each of the points of intersection.
S1≔fsolveBC1,BC2,−4..−3,3..5; S2≔fsolveBC1,BC2,−2.5..−1.5,5..7;S3≔fsolveBC1,BC2,−1..0,−11..9;S4≔fsolveBC1,BC2,0..1,0..1;
S1≔−3.5359390185871604,4.1829937463052169
S2≔−2.0022359204985014,5.9933072030663520
S3≔−0.77342410180917261,−10.240596303289645
S4≔0.11555841921815774,0.31076387545660007
We can now find all of the solutions for the BVP by using the (a, b) pairs as initial values. As expected, S4 (the (a, b) pair in the first quadrant) corresponds to the solution we found using the approxsoln option.
sol1≔dsolveode,y0=S11,y'0=S12,numeric:
odeplotsol1,x=0..1;
sol2≔dsolveode,y0=S21,y'0=S22,numeric:
odeplotsol2,x=0..1;
sol3≔dsolveode,y0=S31,y'0=S32,numeric:
odeplotsol3,x=0..1;
sol4≔dsolveode,y0=S41,y'0=S42,numeric:
odeplotsol4,x=0..1;
We also could have arrived at the solutions if we had used appropriate initial solutions with approxsoln. The following is a list of initial approximations, along with the solutions derived from them.
L≔x−4,3⋅x−2, −10⋅x, x⋅1−x;
L≔x−4,3x−2,−10x,x1−x
sol1≔dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L1:
sol2≔dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L2:
sol3≔dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L3:
sol4≔dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L4:
The dsolve/numeric command handles ODEs that are defined using the piecewise command.
The following example describes the motion of an object that is undergoing a constant downward acceleration due to gravity and a drag force that is proportional to the square of its velocity. The direction of the drag force is opposite to the direction of the object's velocity.
acceleration≔diffht,t,t=piecewisediffht,t≤0,−9.81+6⋅diffht,t2, −9.81−6⋅diffht,t2
acceleration≔ⅆ2ⅆt2ht=−9.81+6ⅆⅆtht2ⅆⅆtht≤0−9.81−6ⅆⅆtht2otherwise
The following assignments set the initial conditions. At t = 0 the object is at the origin (h = 0) and has an initial velocity of 20.
Note: For the initial velocity, we use the differential operator (D) to represent the derivative of h evaluated at a point.
initialPosition≔h0=0;initialVelocity≔Dh0=20;
initialPosition≔h0=0
initialVelocity≔Dh0=20
Now we'll find the solution and plot it.
solution≔dsolveacceleration,initialPosition,initialVelocity,numeric;
solution≔procx_rkf45...end proc
plots:-odeplotsolution,0..1;
The first thing we notice is that the solution handled the piecewise ODE properly. However, it would be nice if we could have the object's velocity reverse when it hits the ground to simulate a bouncing object.
To do this, we'll use the events=[[trigger, action]] option. The events option makes the numerical solver look for the trigger condition between integration steps. When the trigger is found, the solver performs the associated action.
In the following command, the trigger is set to h(t)=0: this fires an event whenever the solution detects the object hitting the ground. The action term is set to diff(h(t),t)=-diff(h(t),t): this reverses the direction of the velocity so that the object bounces up. (For more information on how to set triggers and actions, see dsolve,numeric,Events.)
solution≔dsolveacceleration,initialPosition,initialVelocity,numeric,events=ht=0, diffht,t=−diffht,t;
odeplotsolution,0..1;
In most cases fsolve is sufficient to find all of the roots of your solution. However, there are cases where fsolve will either fail to find a root or you will need to specify a range in order for fsolve to find a root. In such cases, you can use the halt action to stop the solution every time a root is found.
We start by defining an ODE and then plotting the solution to see where the roots for our solution are.
odes≔diffht,t,t=−9.81−6⋅diffht,t2,h0=h0,Dh0=v0;
odes≔ⅆ2ⅆt2ht=−9.81−6ⅆⅆtht2,h0=h0,Dh0=v0
sol≔dsolveodes,numeric,parameters=h0,v0,output=listprocedure;
sol≔t=proct...end proc,ht=proct...end proc,ⅆⅆtht=proct...end proc
We will set the initial position to a very small negative value and the initial velocity to a positive value so that the solution has to cross zero very close to t=0.
solparameters=−0.001,30;
tparameters=−0.001,30=h0=−0.001,v0=30.,htparameters=−0.001,30=h0=−0.001,v0=30.,ⅆⅆthtparameters=−0.001,30=h0=−0.001,v0=30.
Plotting the solution shows that we can expect to find two roots between t=0 and t=0.4.
H≔evalht,sol;
H≔proct...end proc
plotH,0..0.4;
Now try using the fsolve function to find the roots.
fsolveHt=0,t;
0.3983488244
We only found the second of the two roots. To find the first root, we can give fsolve an interval to search in.
fsolveHt=0,t=0..0.1;
fsolveHt=0,t,0..0.1
Again, fsolve failed to find the first root. We will try again, using smaller and smaller intervals.
fsolveHt=0,t=0..0.05;
fsolveHt=0,t,0..0.05
fsolveHt=0,t=0..0.01;
0.00003343371658
It took several tries (and some knowledge of what we were looking for) to find both roots. Obviously, this isn't ideal.
Let's see what happens when we use the events option to find the roots. For this, we'll set the trigger to h(t)=0 and the action to halt. The halt action stops the solution every time the event triggers.
sol2≔dsolveodes,numeric,parameters=h0,v0,events=ht=0,halt;
sol2parameters=−0.001,30;
h0=−0.001,v0=30.
To find the roots, we'll execute the solution at a time beyond the zero-crossings.
sol20.4;
Warning, cannot evaluate the solution further right of .33433716e-4, event #1 triggered a halt
t=0.0000334337167628194,ht=7.1150767569361210−20,ⅆⅆtht=29.8202119178629
The first root was found using the event. Notice that the value of h is not exactly zero. This is because the event is triggered when the solver detects that the event occurred during an integration step, and not necessarily at the exact moment h=0.
Before we can continue to the next root, we have to clear the event using eventclear.
sol2eventclear;
Executing the solution at a time beyond the second zero-crossing should give us the second root for our solution.
Warning, cannot evaluate the solution further right of .39834882, event #1 triggered a halt
t=0.398348824350896,ht=−4.8707104972553510−16,ⅆⅆtht=−29.8202310162904
We have seen from a previous example that dsolve/numeric can solve a piecewise ODE. However, if your piecewise ODE has a condition that depends on your solution being equal to a certain value, the solution could fail. This is because it is extremely unlikely that the solution will be equal to the value at the end of an integration step, and therefore the condition in the piecewise ODE will never be checked.
In this example, dsolve/numeric fails to solve a piecewise ODE because of a condition that needs to be checked when the derivative equals zero. We will see how using an event to force the derivative to be zero fixes this.
This model consists of a spring attached to a block that slides over a table. The table exerts static and kinetic friction on the block. As the block oscillates about equilibrium, it loses energy because of friction and its oscillations decay over time. The block eventually stops at a point where the restoring force from the spring cannot overcome the static friction from the table.
The following piecewise function defines the forces on the block when it comes to rest (see Model description for more information).
forceAtRest ≔ piecewiseabsxt < μs ⋅m⋅gk,0, −k⋅xt+signumxt⋅μk⋅m⋅g;
forceAtRest≔0xt<μsmgk−kxt+signumxtμkmgotherwise
The following piecewise function combines the forces on the block when it is in motion with the forceAtRest function. This gives us the ODE for all the forces acting on the block.
Force≔m⋅diffxt,t,t=piecewisediffxt,t>0, −k⋅xt−μk⋅m⋅g, diffxt,t<0,−k⋅xt+μk⋅m⋅g,diffxt,t=0, forceAtRest;
Force≔mⅆ2ⅆt2xt=−kxt−μkmg0<ⅆⅆtxt−kxt+μkmgⅆⅆtxt<00xt<μsmgk−kxt+signumxtμkmgotherwiseⅆⅆtxt=0
Define equations for the initial position and initial velocity of the block.
initialPosition≔x0=x0;
initialPosition≔x0=x0
initialVelocity≔Dx0=v0;
initialVelocity≔Dx0=v0
Initialize values for the mass (m), the acceleration due to gravity (g), the spring constant (k), and the coefficients of kinetic and static friction (mk and ms, respectfully).
m≔2.0; g≔9.81; k≔50.0; μk≔0.35; μs≔0.40;
m≔2.0
g≔9.81
k≔50.0
μk≔0.35
μs≔0.40
Generate the numerical solution for the ODE with parameters for x0 and v0.
sol1≔dsolveForce, initialPosition, initialVelocity,numeric,parameters= x0,v0;
sol1≔procx_rkf45...end proc
Set the initial position to 2.0 and the initial velocity to zero.
sol1parameters=2.0, 0.0;
x0=2.0,v0=0.
The following command generates a graph of the solution for the range t = 0..5 using plots:-odeplot. The graph includes two dashed lines that show the region where the restoring force of the spring is not strong enough to overcome the static friction from the table. We expect the graph to show the following behavior:
The block will oscillate about the spring's equilibrium point.
The amplitude of the oscillations will decay linearly in time.
The block will remain at rest the first time it comes to a stop between the dashed lines.
odeplotsol1, t,xt, t, μs⋅m⋅gk,color=black,linestyle=dash,t, −μs⋅m⋅gk,color=black,linestyle=dash,0..5;
Warning, cannot evaluate the solution further right of 4.3983772, maxfun limit exceeded (see ?dsolve,maxfun for details)
The solution behaves as expected except for when the block comes to rest between the dashed lines. At this point, the net force on the block is zero and the block should remain at rest.
The reason for this failure is that the numerical solver evaluates the solution (and its derivatives) at discrete integration steps. It is very unlikely that the velocity will be exactly zero at an integration step. Due to this, the zero velocity condition is never checked and the force is never set to zero.
To fix this, use the events=[[trigger, action]] option. Set the trigger to diff(x(t),t)=0.0: this fires an event whenever the solution detects that diff(x(t),t)=0 between integration steps. Set the action term to diff(x(t),t)=0.0: this sets the velocity to zero when an event fires and guarantees that the condition for zero velocity is checked.
sol2≔dsolveForce, initialPosition, initialVelocity, numeric, parameters = x0, v0, events = diffxt, t = 0.0, diffxt, t = 0.0;
sol2parameters=2.0, 0.0;
A graph of the solution now shows the expected behavior when the block comes to rest between the dashed lines.
odeplotsol2, t,xt,t, μs⋅m⋅gk,color=black,linestyle=dash,t, −μs⋅m⋅gk,color=black,linestyle=dash,0..5;
Model description
The block has a mass, m, and is attached to a spring with spring constant, k. The block sits on a level table that exerts friction on it. The static and kinetic friction forces are proportional to the magnitude of the normal force on the block, m⋅g (where g is the acceleration due to gravity). The coefficients of static and kinetic friction are, respectively, ms and mk, with ms ≥ mk .
The block moves in one dimension, x. Take the origin (x=0) to be the position of the block when the spring is at equilibrium. When the spring is stretched, it exerts a restoring force on the block that is given by -k⋅x(t).
The dynamics are similar to that of a block moving on a frictionless table with the following exceptions.
Kinetic friction opposes the motion of the block. The magnitude of the kinetic friction force is:
μk⋅m⋅g
The amplitude of the oscillations decrease linearly over time (see Marchewka, Abbott, and Beichner or Bottega). This is due to the non-conservative kinetic friction force acting on the block.
The block stops momentarily at the end of each oscillation. At this point, the block stays at rest if the spring cannot overcome the force of static friction between the block and the table. Specifically, when the block stops, we must check the following condition:
x(t)<μs⋅m⋅gk
If the preceding condition is true, the block stays at rest. If the preceding condition is not true, then the block starts to move under the effects of the spring's restoring force and the kinetic friction force from the table.
The motion of the block stops completely after a finite amount of time and within the region where the spring cannot overcome the static friction force (see Marchewka, Abbott, and Beichner or Bottega).
The following equation gives the forces acting on the block when it is moving. Note that the sign of the constant-magnitude kinetic friction force is opposite to the sign of the block's velocity.
F__moving={−k⋅xt−μk⋅m⋅g0< ⅆⅆ t xt−k⋅xt+μk⋅m⋅gⅆⅆ t xt< 0
The forces acting on the block when it comes to rest are given by the following equation:
F__rest={0xt<μs⋅m⋅gk−k⋅xt+signumxt⋅μk⋅m⋅gotherwise
The first condition handles the case where the spring cannot overcome the maximum static friction force from the table. Here, the net force on the block is zero and the block stays at rest.
The second condition handles the case where the restoring force from the spring overcomes the static friction from the table and the block starts to move. The signum function returns the sign of the block's position (x(t)) and is opposite in direction to the spring's restoring force.
Combining the equations for Fmoving and Frest gives the complete ODE for the force of the block.
Force={−k⋅xt−μk⋅m⋅g0< ⅆⅆ t xt−k⋅xt+μk⋅m⋅gⅆⅆ t xt< 0{0xt<μs⋅m⋅gk−k⋅xt+signumxt⋅μk⋅m⋅gotherwiseⅆⅆ t xt=0
Return to the example.
References
Marchewka, A., Abbott, D. S., & Beichner, R. J. (2004). Oscillator damped by a constant-magnitude friction force, Am. J. Phys., 72, 477-483. doi: 10.1119/1.1624113
Bottega, W. J. (2006). Engineering vibrations. Boca Raton, FL: CRC Press.
The ODE Analyzer Assistant is a graphical tool for solving ODEs either symbolically or numerically. With it, you can try different options for your ODE, view graphs of the solution, and see the commands used to generate the solution. For more information, see Using the Interactive ODE Analyzer Assistant and the Ordinary Differential Equations (ODEs) section in Chapter 4 of the User Manual.
There are a number of example worksheets in the Maple Help System related to ODEs. To access them, browse to Mathematics > Differential Equations > Example Worksheets in the Maple Help System.
This topic covered some of the basics of using dsolve, but did not cover all the available options. For more information on these options, see the dsolve, dsolve/details, and dsolve/numeric topics.
The numerical examples in this topic used the default numerical solver, rkf45. For your particular ODE, you may need to use a different numerical solver. For more information on the available numerical solvers, see dsolve/numeric, rkf45, ck45, rosenbrock, bvp, rkf45_dae, ck45_dae, rosenbrock_dae, dverk78, lsode, gear, taylorseries, mebdfi, or classical.
A class of IVPs not covered in this topic are "stiff" IVPs. These IVPs usually fail with the default numerical solver, rkf45. For information how to solve a stiff IVP, see Numerical Solution of Stiff IVPs in Maple.
The How Do I... topics cover the essentials for doing mathematics in Maple. Learn more about available tools and features, such as palettes and the context panel.
How Do I...
...Enter a Complex Number?
...Enter a Function?
...Enter a Matrix?
...Enter a Simple Expression?
...Evaluate an Expression?
...Import Tabular Data?
...Plot a Function?
...Plot a Straight Line?
...Plot Multiple Functions?
...Work with Random Generators?
Tools and Features
Palettes
Context Panel
Command Completion
Equation Labels
Assistants
Maple Help
Plotting Guide
Applications
Example Worksheets
Manuals
Refer to Help>Quick Reference for basic Getting Started tips.
See Also
MaplePortal
dsolve
dsolve/numeric
dsolve/numeric/BVP
dsolve/Events
plots/odeplot
Using the Interactive ODE Analyzer Assistant
Download Help Document
