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≔ⅆⅆxy⁡x=a⁢y⁡x
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≔y⁡x=_C1⁢ⅇa⁢x
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≔ⅆⅆxy⁡x=a⁢y⁡x,y⁡0=5
basicSolution≔y⁡x=5⁢ⅇa⁢x
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⁢ⅇa⁢x
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≔ⅆⅆtv⁡t=−1+v⁡t2,v⁡0=0
Run dsolve(equations) to get a solution for v(t).
solution≔dsolveequations;
solution≔v⁡t=−tanh⁡t
odetestsolution,equations;
Assign the solution for v(t) to a variable, V:
V≔evalvt,solution;
V≔−tanh⁡t
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+tanh⁡t2
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≔−ln⁡cosh⁡T
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≔ln⁡2−ln⁡ⅇ2⁢T+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≔ln⁡2−ln⁡ⅇ2⁢t+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≔ⅆⅆtx⁡t=−x⁡t−y⁡t,ⅆⅆty⁡t=x⁡t−y⁡t
Run dsolve(odeSystem) to get a solution. The result is a set of equations that satisfies the system.
systemSol≔dsolveodeSystem;
systemSol≔x⁡t=ⅇ−t⁢_C2⁢cos⁡t+_C1⁢sin⁡t,y⁡t=−ⅇ−t⁢cos⁡t⁢_C1−sin⁡t⁢_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⁢_C2⁢cos⁡t+_C1⁢sin⁡t
Y≔evalyt,systemSol;
Y≔−ⅇ−t⁢cos⁡t⁢_C1−sin⁡t⁢_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≔ⅆⅆtx⁡t=−x⁡t−y⁡t,ⅆⅆty⁡t=x⁡t−y⁡t,x⁡0=1,y⁡0=2
systemSol≔x⁡t=ⅇ−t⁢cos⁡t−2⁢sin⁡t,y⁡t=−ⅇ−t⁢−2⁢cos⁡t−sin⁡t
X≔evalxt,systemSol;
X≔ⅇ−t⁢cos⁡t−2⁢sin⁡t
Y≔−ⅇ−t⁢−2⁢cos⁡t−sin⁡t
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ⅆx2y⁡x+y⁡x=0
sol≔dsolveode;
sol≔y⁡x=_C1⁢sin⁡x+_C2⁢cos⁡x
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≔y⁡0=y1
iv2≔D⁡y⁡0=y2
sol≔dsolveode,iv1,iv2;odetestsol,ode,iv1,iv2;
sol≔y⁡x=y2⁢sin⁡x+y1⁢cos⁡x
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≔a⁢y⁡0+b⁢D⁡y⁡0=f
bc2≔c⋅yPi+d⋅y'Pi=g;
bc2≔c⁢y⁡π+d⁢D⁡y⁡π=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≔y⁡x=−a⁢g+c⁢f⁢sin⁡xa⁢d−b⁢c+b⁢g+d⁢f⁢cos⁡xa⁢d−b⁢c
odetestsolution,ode,bc1,bc2;
Y≔evalyx,solution;
Y≔−a⁢g+c⁢f⁢sin⁡xa⁢d−b⁢c+b⁢g+d⁢f⁢cos⁡xa⁢d−b⁢c
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
y⁡0+2⁢D⁡y⁡0=34
3⁢y⁡π+4⁢D⁡y⁡π=−2
Y;
sin⁡x8+cos⁡x2
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
y⁡0=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≔y⁡x=_C1⁢sin⁡x+cos⁡x
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≔ⅆⅆxy⁡x=1cos⁡y⁡x+1.01,y⁡0=1
sol≔dsolveodesys;
sol≔y⁡x=RootOf⁡101⁢_Z−100⁢x+100⁢sin⁡_Z−101−100⁢sin⁡1
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−101⁢y⁡x100−sin⁡y⁡x+101100+sin⁡1=0
The second is to use the DEtools:-remove_RootOf command.
withDEtools:
sol2≔remove_RootOfsol;
sol2≔101⁢y⁡x−100⁢x+100⁢sin⁡y⁡x−101−100⁢sin⁡1=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−100⁢x+100⁢sin⁡8−100⁢sin⁡1=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≔ⅆⅆxy⁡x=−y⁡x,y⁡0=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,y⁡x=0.0336896836290933
sol0.916;
x=0.916,y⁡x=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;
fsolve⁡sol⁡x=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,y⁡x=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;
x⁡0.9162905219=0.9162905219,y⁡x⁡0.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≔ⅆⅆxy⁡x=a⁢y⁡x,y⁡x0=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,y⁡x=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≔ⅆⅆtx⁡t=a⁢x⁡t−b⁢x⁡t⁢y⁡t,ⅆⅆty⁡t=−c⁢y⁡t+d⁢x⁡t⁢y⁡t,x⁡t0=x0,y⁡t0=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.,x⁡t=20.0599816840952,y⁡t=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,x⁡t=proct...end proc,y⁡t=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≔ⅆⅆtv⁡t=−1+v⁡t2,v⁡0=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≔a⁡t=ⅆⅆtv⁡t,ⅆⅆts⁡t=v⁡t,ⅆⅆtv⁡t=−1+v⁡t2,s⁡0=s0,v⁡0=v0
Generate the solution with parameters for v0 and s0.
solution≔dsolveequations,numeric,parameters=v0,s0,output=listprocedure;
solution≔t=proct...end proc,a⁡t=proct...end proc,s⁡t=proct...end proc,v⁡t=proct...end proc
Set the initial velocity and initial position both to zero.
solutionparameters=0,0;
t⁡parameters=0,0=v0=0.,s0=0.,a⁡t⁡parameters=0,0=v0=0.,s0=0.,s⁡t⁡parameters=0,0=v0=0.,s0=0.,v⁡t⁡parameters=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≔ⅆⅆtv⁡t=−1+v⁡t2,v⁡0=0
sol2≔dsolveequations2,numeric,output=listprocedure;
sol2≔t=proct...end proc,v⁡t=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⁡∫0TV2⁡tⅆ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↦ⅆⅆtV2⁡tt=T|ⅆⅆtV2⁡tt=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ⅆx2y⁡x+y⁡x2⁢ⅆⅆxy⁡x+cos⁡y⁡x=x
bc1≔1+y0⋅y'0−3⋅y0=0;
bc1≔1+y⁡0⁢D⁡y⁡0−3⁢y⁡0=0
bc2≔ 1+y1⋅y'1−3⋅y1−y12+0.5=0;
bc2≔1+y⁡1⁢D⁡y⁡1−3⁢y⁡1−y⁡12+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+h⁡a,b[1]*h⁡a,b[2] − 3*h⁡a,b[1] − h⁡a,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,3⁢x−2,−10⁢x,x⁢1−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ⅆt2h⁡t=−9.81+6⁢ⅆⅆth⁡t2ⅆⅆth⁡t≤0−9.81−6⁢ⅆⅆth⁡t2otherwise
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≔h⁡0=0
initialVelocity≔D⁡h⁡0=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ⅆt2h⁡t=−9.81−6⁢ⅆⅆth⁡t2,h⁡0=h0,D⁡h⁡0=v0
sol≔dsolveodes,numeric,parameters=h0,v0,output=listprocedure;
sol≔t=proct...end proc,h⁡t=proct...end proc,ⅆⅆth⁡t=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;
t⁡parameters=−0.001,30=h0=−0.001,v0=30.,h⁡t⁡parameters=−0.001,30=h0=−0.001,v0=30.,ⅆⅆth⁡t⁡parameters=−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;
fsolve⁡H⁡t=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;
fsolve⁡H⁡t=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,h⁡t=7.11507675693612⁢10−20,ⅆⅆth⁡t=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.