Solve an ordinary differential equation - Maple Help
For the best experience, we recommend viewing online help using Google Chrome or Microsoft Edge.

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Getting Started : How Do I... : Solve an ordinary differential equation

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:

1. 

Defining your system of ODEs.

2. 

Finding solutions.

3. 

Assigning the solutions to names (thus making them easier to plot and investigate).

4. 

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

Solving an ODE Symbolically

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=ay.

Define a variable called basicODE and set it equal to the ODE.

restart;

basicODEdiffyx,x=ayx;

basicODEⅆⅆxyx=ayx

(1.1.1)

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.

basicSolutiondsolvebasicODE;

basicSolutionyx=_C1ⅇax

(1.1.2)

 

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.

basicODEdiffyx,x=ayx, y0=5;

basicODEⅆⅆxyx=ayx,y0=5

(1.1.3)

basicSolutiondsolvebasicODE;

basicSolutionyx=5ⅇax

(1.1.4)

 

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

(1.1.5)

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;Yevalyx,basicSolution;

Y5ⅇax

Y5ⅇax

(1.1.6)

After assigning a value to the parameter a, use eval on Y to get an exact value for the solution at a point.

a13;

a13

(1.1.7)

evalY,  x=5;

5ⅇ53

(1.1.8)

Use the evalf function to get a floating-point value.

evalfevalY, x=5;

0.9443780140

(1.1.9)

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.

bfsolveYx=2,x;

b2.748872196

(1.1.10)

evalY,x=b;

2.000000000

(1.1.11)

Taking Derivatives and Integrals of Symbolic Solutions

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).

restart;

equationsdiffvt,t=1+vt2, v0=0

equationsⅆⅆtvt=1+vt2,v0=0

(1.2.1)

Run dsolve(equations) to get a solution for v(t).

solutiondsolveequations;

solutionvt=tanht

(1.2.2)

odetestsolution,equations;

0

(1.2.3)

Assign the solution for v(t) to a variable, V:

Vevalvt,solution;

Vtanht

(1.2.4)

You can use the diff function to take the derivative of V and get the acceleration of the object, A.

AdiffV,t;

A1+tanht2

(1.2.5)

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.)

SintV,t=0..T;

SlncoshT

(1.2.6)

If you are only interested in solutions for T > 0, you can specify this by using assuming on the int expression.

SintV,t=0..T assuming T>0;

Sln2lnⅇ2T+1+T

(1.2.7)

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:

SevalS,T=t;

Sln2lnⅇ2t+1+t

(1.2.8)

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;

Solving a System of ODEs

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.

restart;

odeSystem diffxt,t=xt yt,  diffyt,t=  xtyt

odeSystemⅆⅆtxt=xtyt,ⅆⅆtyt=xtyt

(1.3.1)

Run dsolve(odeSystem) to get a solution. The result is a set of equations that satisfies the system.

systemSoldsolveodeSystem;

systemSolxt=ⅇt_C2cost+_C1sint,yt=ⅇtcost_C1sint_C2

(1.3.2)

odetestsystemSol,odeSystem;

0

(1.3.3)

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

(1.3.4)

Yevalyt,systemSol;

Yⅇtcost_C1sint_C2

(1.3.5)

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=  xtyt,  x0=1,  y0=2

odeSystemⅆⅆtxt=xtyt,ⅆⅆtyt=xtyt,x0=1,y0=2

(1.3.6)

systemSoldsolveodeSystem;

systemSolxt=ⅇtcost2sint,yt=ⅇt2costsint

(1.3.7)

Xevalxt,systemSol;

Xⅇtcost2sint

(1.3.8)

Yevalyt,systemSol;

Yⅇt2costsint

(1.3.9)

Use eval and evalf to analyze your solution.

evalX, t=Pi; evalfevalX, t=Pi; 

ⅇπ

−0.04321391825

(1.3.10)

evalY, t=Pi; evalfevalY, t=Pi;

2ⅇπ

−0.08642783650

(1.3.11)

The following command finds the first time x(t) = y(t) with t ≥ 1.5.

fsolveXt=Yt,t=1.5;

2.819842099

(1.3.12)

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.

restart;

odediffyx,x,x+yx=0;

odeⅆ2ⅆx2yx+yx=0

(1.4.1)

soldsolveode;

solyx=_C1sinx+_C2cosx

(1.4.2)

The general solution has two constants, _C1 and _C2. To get a specific solution, include two initial value equations along with ode in dsolve.

iv1y0=y1; iv2y'0=y2;

iv1y0=y1

iv2Dy0=y2

(1.4.3)

soldsolveode,iv1,iv2;odetestsol,ode,iv1,iv2;

solyx=y2sinx+y1cosx

0

(1.4.4)

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.

bc1ay0+by'0=f;

bc1ay0+bDy0=f

(1.4.5)

bc2cyPi+dy'Pi=g;

bc2cyπ+dDyπ=g

(1.4.6)

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 ad - bc ≠ 0.)

solutiondsolveode,bc1,bc2;

solutionyx=ag+cfsinxadbc+bg+dfcosxadbc

(1.4.7)

odetestsolution,ode,bc1,bc2;

0

(1.4.8)

Yevalyx,solution;

Yag+cfsinxadbc+bg+dfcosxadbc

(1.4.9)

To analyze the solution, assign values to the parameters a, b, c, d, f, and g (ensuring, of course, that ad - bc ≠ 0).

a1; b2; c3; d4; f34; g2; bc1; bc2;

a1

b2

c3

d4

f34

g−2

y0+2Dy0=34

3yπ+4Dyπ=−2

(1.4.10)

Y;

sinx8+cosx2

(1.4.11)

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 ad - bc = 0. We do this by setting the boundary values to y(0) = 1 and y(π) = 0.

a1; b0; c1; d0; f1; g0; bc1; bc2;

a1

b0

c1

d0

f1

g0

y0=1

yπ=0

(1.4.12)

soldsolveode,bc1,bc2;

sol

(1.4.13)

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.

g1; bc1; bc2;

g−1

y0=1

yπ=−1

(1.4.14)

soldsolveode, bc1, bc2;

solyx=_C1sinx+cosx

(1.4.15)

odetestsol,ode,bc1,bc2;

0

(1.4.16)

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.

Implicit Solutions with RootOf

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.

restart;

odesysdiffyx,x=1cosyx+1.01, y0=1;

 

odesysⅆⅆxyx=1cosyx+1.01,y0=1

(1.5.1)

soldsolveodesys;

solyx=RootOf101_Z100x+100sin_Z101100sin1

(1.5.2)

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;

0

(1.5.3)

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).

sol1dsolveodesys,implicit;

sol1x101yx100sinyx+101100+sin1=0

(1.5.4)

The second is to use the DEtools:-remove_RootOf command.

withDEtools:

sol2remove_RootOfsol;

sol2101yx100x+100sinyx101100sin1=0

(1.5.5)

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).

eqevalsol2, yx=8;

eq707100x+100sin8100sin1=0

(1.5.6)

fsolveeq,x;

7.217887262

(1.5.7)

Solving an ODE Numerically

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

Basic Example: y' = -y

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.

restart;

basicOdediffyx,x = yx, y0=5;

basicOdeⅆⅆxyx=yx,y0=5

(2.1.1)

 

To get a numerical solution, use dsolve with the numeric parameter.

soldsolvebasicOde,numeric;

solprocx_rkf45...end proc

(2.1.2)

 

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

(2.1.3)

sol0.916;

x=0.916,yx=2.00058112757547

(2.1.4)

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.

 

withplots:

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

(2.1.5)

To solve this, use the output=listprocedure option with dsolve, and then assign the solution to a variable with the eval(expression,equations) command.

soldsolvebasicOde,numeric,output=listprocedure;

solx=procx...end proc,yx=procx...end proc

(2.1.6)

Yevalyx,sol;

Yprocx...end proc

(2.1.7)

afsolveYx=2,x;

a0.9162905219

(2.1.8)

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

(2.1.9)

Ya;

1.99999999996477

(2.1.10)

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;

Basic Example with Parameters: y' = a y

This example shows how to use dsolve get a numerical solution to ⅆⅆ x y=ay 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.

restart;

basicOde2diffyx,x = ayx, yx0=y0;

basicOde2ⅆⅆxyx=ayx,yx0=y0

(2.2.1)

 

Use dsolve to get a numerical solution. Also, because basicOde2 contains the parameters a, x0, and y0,use the parameters option.

sol2dsolvebasicOde2,numeric, parameters=a, x0, y0;

sol2procx_rkf45...end proc

(2.2.2)

 

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.

(2.2.3)

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;

x=5.0,yx=0.0336896836290933

(2.2.4)

The following graphs sol2 over x=0..5.

withplots:

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.

(2.2.5)

odeplotsol2,x=0..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.

sol2dsolvebasicOde2,numeric, parameters=a, x0, y0,output=listprocedure;

sol2x=procx...end proc,yx=procx...end proc

(2.2.6)

Yevalyx,sol2;

Yprocx...end proc

(2.2.7)

Before you can use the function, you must set values for the parameters a, x0, and y0.

Yparameters=1,0,5;

a=−1.,x0=0.,y0=5.

(2.2.8)

You can now use fsolve on the solution. The following finds the value of x for which the solution equals 2.

bfsolveYx=2,x;

b0.9162905219

(2.2.9)

Yb;

1.99999999996477

(2.2.10)

Solving a System of ODEs Numerically

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.

restart;

desdiffxt,t=a xtb xtyt,diffyt,t=c yt+d xtyt, xt0=x0, yt0=y0;

desⅆⅆtxt=axtbxtyt,ⅆⅆtyt=cyt+dxtyt,xt0=x0,yt0=y0

(2.3.1)

 

Generate the solution with a list of parameters.

desSoldsolvedes,numeric,parameters=a,b,c,d,t0,x0,y0;

desSolprocx_rkf45...end proc

(2.3.2)

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.

(2.3.3)

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

(2.3.4)

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.

withplots:

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.

desSoldsolvedes,numeric,parameters=a,b,c,d,t0,x0,y0,output=listprocedure;

desSolt=proct...end proc,xt=proct...end proc,yt=proct...end proc

(2.3.5)

Xevalxt,desSol;

Xproct...end proc

(2.3.6)

Yevalyt,desSol;

Yproct...end proc

(2.3.7)

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;

a=0.1,b=0.01,c=0.025,d=0.001,t0=0.,x0=50.,y0=15.

(2.3.8)

Yparameters;

a=0.1,b=0.01,c=0.025,d=0.001,t0=0.,x0=50.,y0=15.

(2.3.9)

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.

afsolveXt=Yt,t;

a16.82819661

(2.3.10)

Xa;Ya;

16.5912730797044

16.5912730836271

(2.3.11)

Taking Derivatives and Integrals of Numeric Solutions

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.

restart;

velocityEqsdiffvt,t=1+vt2,v0=v0;

velocityEqsⅆⅆtvt=1+vt2,v0=v0

(2.4.1)

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.

equationsvelocityEqs, at=diffvt,t, diffst,t=vt, s0=s0

equationsat=ⅆⅆtvt,ⅆⅆtst=vt,ⅆⅆtvt=1+vt2,s0=s0,v0=v0

(2.4.2)

Generate the solution with parameters for v0 and s0.

solutiondsolveequations,numeric,parameters=v0,s0,output=listprocedure;

solutiont=proct...end proc,at=proct...end proc,st=proct...end proc,vt=proct...end proc

(2.4.3)

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.

(2.4.4)

Assign the velocity, acceleration, and position to the variables V, A, and S, respectively.

Vevalvt,solution;

Vproct...end proc

(2.4.5)

Aevalat,solution;

Aproct...end proc

(2.4.6)

Sevalst,solution;

Sproct...end proc

(2.4.7)

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.

equations2diffvt,t=1+vt2, v0=0;

equations2ⅆⅆtvt=1+vt2,v0=0

(2.4.1.1)

sol2dsolveequations2,numeric,output=listprocedure;

sol2t=proct...end proc,vt=proct...end proc

(2.4.1.2)

V2evalvt,sol2;

V2proct...end proc

(2.4.1.3)

Next, we'll define SInt by integrating our numerical solution, V2.

SIntTevalfIntV2t,t=0..T;

SIntTevalf0TV2tⅆt

(2.4.1.4)

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;evalflncosh6;

 

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

(2.4.1.5)

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.

ADiffTevalDiffV2t,t,t=T;

ADiffTⅆⅆtV2tt=T|ⅆⅆtV2tt=T

(2.4.1.6)

The following commands compare ADiff with A and the exact solution, tanh(t)2−1.

CodeTools:-UsageevalfADiff6; CodeTools:-UsageA6; evalftanh621;

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

(2.4.1.7)

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.

Techniques for Solving a Nonlinear Two-Point BVP

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).

restart;

ode  diffyx,x,x+yx2diffyx,x+cosyx=x;

odeⅆ2ⅆx2yx+yx2ⅆⅆxyx+cosyx=x

(2.5.1)

bc11+y0y'03y0=0;

bc11+y0Dy03y0=0

(2.5.2)

bc2 1+y1y'13y1y12+0.5=0;

bc21+y1Dy13y1y12+0.5=0

(2.5.3)

We want to find all of the solutions to this system in the interval 0x1.

A first attempt at using dsolve/numeric fails.

solutiondsolveode,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.

solutiondsolveode,bc1,bc2,numeric, approxsoln=yx=x;

solutionprocx_bvp...end proc

(2.5.4)

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+y0y'03y0=0

1+y1y'13y1y12+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+ab3a=0

1+y1y'13y1y12+0.5=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.

hproca,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.          Fdsolveode,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+ab3a=0

1+ha,b1ha,b23ha,b1ha,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.

BC1proca,b      1+ab3a end proc;

BC1proca,b1+a*b3*aend proc

(2.5.5)

BC2proca,b      1+ha,b1ha,b23ha,b1ha,b12+0.5 end proc;

BC2proca,b1+ha,b[1]*ha,b[2]3*ha,b[1]ha,b[1]^2+0.5end proc

(2.5.6)

 

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).

withplots:

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.

S1fsolveBC1,BC2,4..3,3..5; S2fsolveBC1,BC2,2.5..1.5,5..7;S3fsolveBC1,BC2,1..0,11..9;S4fsolveBC1,BC2,0..1,0..1;

S1−3.5359390185871604,4.1829937463052169

S2−2.0022359204985014,5.9933072030663520

S3−0.77342410180917261,−10.240596303289645

S40.11555841921815774,0.31076387545660007

(2.5.7)

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.

sol1dsolveode,y0=S11,y'0=S12,numeric:

odeplotsol1,x=0..1;

sol2dsolveode,y0=S21,y'0=S22,numeric:

odeplotsol2,x=0..1;

sol3dsolveode,y0=S31,y'0=S32,numeric:

odeplotsol3,x=0..1;

sol4dsolveode,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.

Lx4,3x2, 10x, x1x;

Lx4,3x2,10x,x1x

(2.5.8)

sol1dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L1:

odeplotsol1,x=0..1;

sol2dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L2:

odeplotsol2,x=0..1;

sol3dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L3:

odeplotsol3,x=0..1;

sol4dsolveode,bc1,bc2,yx,numeric,approxsoln=yx=L4:

odeplotsol4,x=0..1;

Solving Piecewise ODEs and an Introduction to Event Handling

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.

restart;

accelerationdiffht,t,t=piecewisediffht,t0,9.81+6diffht,t2, 9.816diffht,t2

accelerationⅆ2ⅆt2ht=9.81+6ⅆⅆtht2ⅆⅆtht09.816ⅆⅆtht2otherwise

(2.6.1)

 

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.

initialPositionh0=0;initialVelocityDh0=20;

initialPositionh0=0

initialVelocityDh0=20

(2.6.2)

Now we'll find the solution and plot it.

solutiondsolveacceleration,initialPosition,initialVelocity,numeric;

solutionprocx_rkf45...end proc

(2.6.3)

withplots:

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.)

solutiondsolveacceleration,initialPosition,initialVelocity,numeric,events=ht=0, diffht,t=diffht,t;

solutionprocx_rkf45...end proc

(2.6.4)

odeplotsolution,0..1;

 

Finding Roots Using Events and the Halt Action

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.

restart;

odesdiffht,t,t=9.816diffht,t2,h0=h0,Dh0=v0;

odesⅆ2ⅆt2ht=9.816ⅆⅆtht2,h0=h0,Dh0=v0

(2.7.1)

soldsolveodes,numeric,parameters=h0,v0,output=listprocedure;

solt=proct...end proc,ht=proct...end proc,ⅆⅆtht=proct...end proc

(2.7.2)

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.

(2.7.3)

Plotting the solution shows that we can expect to find two roots between t=0 and t=0.4.

Hevalht,sol;

Hproct...end proc

(2.7.4)

plotH,0..0.4;

Now try using the fsolve function to find the roots.

fsolveHt=0,t;

0.3983488244

(2.7.5)

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

(2.7.6)

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

(2.7.7)

fsolveHt=0,t=0..0.01;

0.00003343371658

(2.7.8)

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.

sol2dsolveodes,numeric,parameters=h0,v0,events=ht=0,halt;

sol2procx_rkf45...end proc

(2.7.9)

sol2parameters=0.001,30;

h0=−0.001,v0=30.

(2.7.10)

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

(2.7.11)

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.