This activity shows how to use Sage to solve differential equations.

Let’s first consider the problem of finding the general solution of an ODE symbolically. As an example, let’s consider the differential equation:

The Sage function is used to find symbolic solutions of ODEs. The following Sage cell shows how to solve the ODE above:

Clicking the Evaluate button will execute the cell and compute the general solution of the differential equation.

In the answer box below, enter the general solution of the differential equation

Use the variable for the constant of integration.

The general solution of the DE is: .

Hint: If you are getting an error, you are probably trying to enter as your constant of integration. Use the variable instead.

We can also use to find the solution of initial value problems. Let’s, for example, consider the problem:

The following Sage cell demonstrates how to solve this problem:

The solution of the initial value problem is: .

If we need a single numerical value of the solution, we can compute it as shown in the cell below:

Use the Sage cell above to compute the following values for the solution of the IVP:

More often than not, we need several values of the solution of the ODE. The following Sage cell demonstrates how to

The output of the previous cell is a list, with the values of the solution in a specified range. The range is defined by the Sage expression: This expression defines the following range of values in Sage: Notice that the final value of the range is , not as one might expect. This is a convention that Sage inherits from Python, where the rightmost element of a range expression is never actually included in the range. On way to think of it is that, in the expression , the set of values specified is a subset of the semi-open interval .

Let’s now plot the solution of the differential equation:

Answer the following question based on the plot generated by the Sage cell above:

The value of is in the interval: , , , , , , ,

In an applied problem, a symbolic solution will rarely be available, and one must resort to numerical solutions. As an example, let’s consider the following model for a population with harvesting: In this equation, the parameter represents a harvesting rate (in units of number of individuals per unit of time). Let’s first try to solve this equation symbolically in Sage, for the parameter values , , and the initial condition :

Evaluating the cell above, we can see that the symbolic solution looks quite complicated. Notice that the function that represents the solution, , has not been isolated in the left-hand side of the equation. In other words, this the solution is given in implicit form. In cases like this, it may be preferable to use a numerical method. This can be done in Sage as shown in the following cell, using the function . Notice that the interface for this function differs in a significant way from the function

From the plot generated by the Sage cell above, we can conclude that the population :
Grows without bound as . Stabilizes near a certain positive value as . Becomes extinct at a certain finite value of . Oscillates around some value as .
Change the value of in the cell above to 1.5 and run the cell again. Experiment with the maximum value of specified in the statement: Notice that, as you increase the maximum value of , you may also want to increase the value of the step-size. Answer the following question based on your experiments:

For the harvesting rate the population :

Grows without bound as . Stabilizes near a certain positive value as . Becomes extinct at a certain finite value of . Oscillates around some value as .
It is of interest to determine what is the largest value of the harvesting rate we can choose so that we have a stable population in the model. Experiment with different values in the Sage cell above to determine this value approximately to three decimal places.

Notice that you may also have to adjust the maximum value of the time variable , since the differential equation becomes unstable for certain values of , which leads to numerical errors in the computation.

The maximum value of the harvesting rate, to three decimal places is: