Matthew L. Wright
Visiting Assistant Professor, St. Olaf College

Mathematica Examples

Math 230: Differential Equations

The following examples show how to use Mathematica to perform calculations that are common in Math 230. For more information about any of the Mathematica commands used below, see the Wolfram language documentation.

1. Slope Fields

The VectorPlot command can draw a slope field for the differential equation \( \frac{dy}{dt} = f(t,y) \). Since VectorPlot draws a vector field (not the same as a slope field) by default, it is necessary to specify some optional parameters to obtain a slope field. The general syntax is as follows:

VectorPlot[ {1, \(f(t,y) \)}, {t,\(t_{min}\),\(t_{max}\)}, {y,\(y_{min}\),\(y_{max}\)}, VectorStyle->"Segment", VectorScale->{Tiny, Tiny, None} ]

For example, the following command draws a slope field for \( \frac{dy}{dt} = ty \) in the window \(-3 \le t \le 3\), \(-3 \le y \le 3\):

VectorPlot[ {1, t*y}, {t,-3,3}, {y,-3,3}, VectorStyle->"Segment", VectorScale->{Tiny, Tiny, None} ]

The output of the previous command should look like this:

slope field

See also: Mathematica documentation for VectorPlot

2. Euler's Method

Using a loop

For example, suppose we want to approximate \(y(1)\), where \(y(t)\) is the solution to \(\frac{dy}{dt} = 2y - \sin(t)\), with \(y(0)=3\), using Euler's method with step size \(\Delta t = 0.01\). One way to do this is to use a loop to compute the successive approximations, as shown here:

f[ t_, y_ ] = 2*y - Sin[t];
y[0] = 3;
Δt = 0.01;
Do[ y[ n+1 ] = f[ Δt * n, y[n] ] *Δt + y[n], {n,0,99} ]

The semicolons in the code above surpress the output of the first three lines. However, after the above code is run, y[100] constains the value 20.7345, which is our approximation of \(y(1)\) using Euler's method with 100 steps.

To see the approximations of \(y\) at the intermediate steps of Euler's method, use the following command:

Table[ {n, y[n]}, {n, 0, 100} ]

The above command produces a list of pairs \( (n, y_n) \) such that \( y_n \approx y(n \cdot \Delta t) \) is the approximation found by step \(n\) Euler's method.

Using NDSovle

The NDSolve command finds numerical solutions to differential equations. It can do so via Euler's method, among other methods.

Suppose you have the initial-value problem \( y'(t) = f(t,y) \) with \(y(0) = y_0\), and you want to approximate \(y(t_1)\) using Euler's method with step size \( \Delta t\). The general syntax to do this with NDSolve is as follows:

NDSolve[ {y'[t] = \(f(t,y) \), y[0] == \(y_0\)}, y, {t, 0, \(t_1\)}, StartingStepSize -> \( \Delta t \), Method -> "ExplicitEuler" ]

Note that NDSolve returns a nested list containing an InterpolatingFunction object.

For example, if \(\frac{dy}{dt} = 2y - \sin(t)\), with \(y(0)=3\). To approximate \(y(1)\) using Euler's method with step size \(\Delta t = 0.01\), use the following command:

NDSolve[ {y'[t] == f[t, y[t]], y[0] == 3}, y, {t, 0, 1}, StartingStepSize -> 0.01, Method -> "ExplicitEuler" ]

The output of the above line looks like this:

NDSolve output

Notice that the NDSolve command has returned an InterpolatingFunction inside of a nested list. To extract the function from the nested list, use the /. (the short form of ReplaceAll) as follows:

y = y /. First[ NDSolve[ {y '[t] == f[t, y[t]], y[0] == 3}, y, {t, 0, 1}, StartingStepSize -> 0.01, Method -> "ExplicitEuler" ] ]

The solution (in the form of an InterpolatingFunction) is now stored in the variable y. For example, y[1] gives you the value of the approximation of \( y(1) \). You can also list the approximations of \( y \) at the intermediate steps of Euler's method with this command:

Table[ {t, y[t]}, {t, 0, 1, 0.01} ]

...or plot the graph of the approximate solution like this:

Plot[ y[t], {t, 0, 1} ]

See also: Mathematica documentation for NDSolve

3. Finding Analytic Solutions

The DSolve command can find analytic solutions to some differential equations and initial-value problems.

For example, to find the general solution to \( \frac{dy}{dt} = ty \):

DSolve[ y'[t] == t*y[t], y[t], t ]

To find the particular solution to \( \frac{dy}{dt} = ty \) with initial condition \( y(0)=5 \):

DSolve[ {y'[t] == t*y[t], y[0] == 5}, y[t], t ]

The previous command produces the following output:


Notice that solution is returned as a mathematical rule inside of a nested list. We often want to extract the solution from the nested list—for example, to plot the solution. Here is one way to do this:

sol = DSolve[ {y'[t] == t*y[t], y[0] == 5}, y[t], t ];
Plot[ Evaluate[ y[t] /. sol ], {t, -3,3} ]

The first line above solves the initial-value problem and stores the solution (which is a nested list containng a rule) in the variable sol. We then use /. (the short form of ReplaceAll) to apply the rule contained in the nested list. Since we want to evaluate the expression y[t] /. sol inside of the Plot function, we must use Evaluate[ y[t] /. sol ]. The output is the following plot:

solution plot

See also: Mathematica documentation for DSolve