The material below is a small excerpt/appetizer from the online course Introduction to Differential Equations - which is just one of many online courses available at Wolfram U.

To carry on with the full course or any of the other course, simply copy the Wolfram Language/Mathematica expressions into Jupyter notebooks with a Wolfram Language kernel of your own.

Introduction to Differential Equations

What Is a Differential Equation?

See https://www.wolframcloud.com/obj/online-courses/introduction-to-differential-equations/what-is-a-differential-equation.html

Direction Fields

Overview

When studying a differential equation describing a mathematical model, it is often desirable to make conclusions about the behavior of the solution. This can be done using a valuable tool called a direction field. One major advantage is that it is not necessary to solve the differential equation in order to construct a direction field. Direction fields are a useful first step in the investigation of a differential equation.

image.png

Direction Fields

Direction fields are a useful tool for differential equations of the form:

y'(x)=f(x,y(x))

To construct a direction field, plot a line segment at each (x,y) coordinate where the slope of the segment is the value of the derivative evaluated at that point. As an example, consider the differential equation:

y'(x)=sin(x)+y(x)

The Wolfram Language function VectorPlot can be used to visualize the direction field:

In [13]:
VectorPlot[{1,Sin[x]+y},{x,-3,3},{y,-3,3},FrameLabel->{Style["x",18],Style["y",18]}]
Out[13]:
Output

At different (x,y) coordinates, the derivative

y'(x)=sin(x)+y(x)

is evaluated and represented as a line segment. For example, when x=0,y=0 , the derivative evaluates to:

y'=sin(0)+0=0

Notice in the VectorPlot that solutions tend to diverge to positive and negative infinity. This is a direct result of the fact that the derivative is largely determined by the y value.

Population Model

Suppose a government agency needs to model a population in a closed environment. It is fair to assume that the change in the population is directly proportional to the number present. Using t as a variable for time, p(t) as the population and r as the proportionality constant, the relationship can be written as:

p'(t)=rp(t)

Use VectorPlot with r=0.5 to visualize the direction field:

image.png

Each segment on the direction field corresponds to the slope of the solution at that particular point. When there is a large population, there are more people giving birth, which increases the population at a faster rate. When there are no people, there is an equilibrium solution where the population does not change with time.

In [5]:
VectorPlot[{1,0.5 p},{t,0,5},{p,0,5},InterpretationBox[DynamicModuleBox[{Set[Typeset`open, False]}, TemplateBox[{"Expression", StyleBox["\"Options\"", "IconizedCustomName", StripOnInput -> False], GridBox[{{ItemBox[RowBox[{TagBox["\"Head: \"", "IconizedLabel"], "\:f360", TagBox["Rule", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Byte count: \"", "IconizedLabel"], "\:f360", TagBox["320", "IconizedItem"]}]]}}, GridBoxAlignment -> {"Columns" -> {{Left}}}, DefaultBaseStyle -> "Column", GridBoxItemSize -> {"Columns" -> {{Automatic}}, "Rows" -> {{Automatic}}}], Dynamic[Typeset`open]}, "IconizedObject"]], FrameLabel -> {Style["t", 18], Style["p(t)", 18]}, SelectWithContents -> True, Selectable -> False]]
Out[5]:
Output

Population Model with Predation

Consider now a population of rabbits. Assume that in addition to the population of rabbits, there is a population of wolves that are preying on the rabbits. The previous model can be improved by noting that the population is decreasing at a constant predation rate d:

p'(t)=rp(t)-d

Visualize the direction field with VectorPlot using

r=0.5,d=25:

In [6]:
VectorPlot[{1,0.5 p-25},{t,0,5},{p,45,55},InterpretationBox[DynamicModuleBox[{Set[Typeset`open, False]}, TemplateBox[{"List", StyleBox["\"Options\"", "IconizedCustomName", StripOnInput -> False], GridBox[{{ItemBox[RowBox[{TagBox["\"Head: \"", "IconizedLabel"], "\:f360", TagBox["List", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Length: \"", "IconizedLabel"], "\:f360", TagBox["1", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Byte count: \"", "IconizedLabel"], "\:f360", TagBox["384", "IconizedItem"]}]]}}, GridBoxAlignment -> {"Columns" -> {{Left}}}, DefaultBaseStyle -> "Column", GridBoxItemSize -> {"Columns" -> {{Automatic}}, "Rows" -> {{Automatic}}}], Dynamic[Typeset`open]}, "IconizedObject"]], {FrameLabel -> {Style["Time", 18], Style["Rabbit Population ", 18]}}, SelectWithContents -> True, Selectable -> False]]
Out[6]:
Output

When the initial population at t=0 is below 50, the rate of predation is higher than the birth rate, and the rabbits become extinct. There is a new equilibrium solution when the initial population starts at 50. Solutions tend to diverge from the equilibrium solution above that population. This is an unstable equilibrium point.

Falling Object Model

Consider an object dropped from a large height. There are two major forces acting on this object: • The force of gravity pulling it downward. • The force of air resistance resisting the movement in the downward direction. From physics, the sum of all forces acting on an object is equal to the object’s mass multiplied by the object’s acceleration:

mv'(t)=(forceofgravity)-(airresistance)

Assume that the air resistance is proportional to the velocity and that the force of gravity is constant:

mv'(t)=mg-αv(t)

With parameters m=50 and α=10 , the falling object model becomes:

50v'(t)=50*9.81-10v(t)

Observe the direction field of the falling object using VectorPlot:

In [8]:
VectorPlot[{1,9.81-v/5},{t,0,40},{v,45,55},InterpretationBox[DynamicModuleBox[{Set[Typeset`open, False]}, TemplateBox[{"List", StyleBox["\"Options\"", "IconizedCustomName", StripOnInput -> False], GridBox[{{ItemBox[RowBox[{TagBox["\"Head: \"", "IconizedLabel"], "\:f360", TagBox["List", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Length: \"", "IconizedLabel"], "\:f360", TagBox["1", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Byte count: \"", "IconizedLabel"], "\:f360", TagBox["376", "IconizedItem"]}]]}}, GridBoxAlignment -> {"Columns" -> {{Left}}}, DefaultBaseStyle -> "Column", GridBoxItemSize -> {"Columns" -> {{Automatic}}, "Rows" -> {{Automatic}}}], Dynamic[Typeset`open]}, "IconizedObject"]], {FrameLabel -> {Style["Time", 18], Style["Velocity", 18]}}, SelectWithContents -> True, Selectable -> False]]
Out[8]:
Output

There is an equilibrium solution near v(t)=50 . Solutions in the falling object model tend to converge to the equilibrium solution. This describes a stable equilibrium point.

The equilibrium solution can be exactly solved for by setting v'(t)=0 and solving for v(t) :

In [9]:
Solve[0\[Equal]50*9.81-10*v[t],v[t]]
Out[9]:
Output

From this result, conclude that there exists an equilibrium solution when v(t)=49.05 :

In [10]:
Show[Plot[49.05,{t,0,40},InterpretationBox[DynamicModuleBox[{Set[Typeset`open, False]}, TemplateBox[{"List", StyleBox["\"Options\"", "IconizedCustomName", StripOnInput -> False], GridBox[{{ItemBox[RowBox[{TagBox["\"Head: \"", "IconizedLabel"], "\:f360", TagBox["List", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Length: \"", "IconizedLabel"], "\:f360", TagBox["3", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Byte count: \"", "IconizedLabel"], "\:f360", TagBox["776", "IconizedItem"]}]]}}, GridBoxAlignment -> {"Columns" -> {{Left}}}, DefaultBaseStyle -> "Column", GridBoxItemSize -> {"Columns" -> {{Automatic}}, "Rows" -> {{Automatic}}}], Dynamic[Typeset`open]}, "IconizedObject"]], {PlotRange -> {{0, 40}, {40, 60}}, AxesLabel -> {Style["t", 18], Style["v(t)", 18]}, PlotStyle -> Thickness[Large]}, SelectWithContents -> True, Selectable -> False]],VectorPlot[{1,9.81-v/5},{t,0,40},{v,40,60}]]
Out[10]:
Output

Leaking Tank Model

Consider a tank with 100 L of salt water at a concentration of 3 g/L. A mixture of salt water (5 g/L) is being poured into the tank at a rate of 2 L/min. The tank has a hole where it is leaking at a rate of 0.15 L/min. Assuming that the salt is always evenly distributed inside the tank, determine the amount of salt over time.

image.png

First, construct a model to describe the amount of salt s(t) in the tank at any time t :

s'(t) = (rate of salt entering) - (rate of salt exiting)

To determine the rate of salt entering, look at the concentration of salt in the mixture multiplied by the volume of mixture entering:

(rate of salt entering) = concentration volume = 5 g / 1 L 2 L / 1 min = 10 g / 1min

Thus there is 10 g of salt entering per minute. In order to calculate the rate of salt exiting, again look at the concentration of salt in the mixture multiplied by the volume of mixture exiting:

concentration = (amount of salt in tank) / (volume of mixture)

The volume of mixture starts at 100 L, with 2 L being poured in and 0.15 L leaking out:

Volume = 100 + (2-.15)t

The concentration of salt in the tank can be calculated as:

concentration = s(t) / (100+1.85t)

Multiply the concentration by the volume of mixture exiting to get an expression for the rate of salt leaving the tank:

rate of salt exiting = (concentration) (volume) = s(t) / (100+1.85t) 0.15

Combining the salt entering and salt exiting, the model can be written as:

s'(t) = 10 - 0.15 s(t) / (100+1.85t)

The model under consideration is:

s'(t) = 10 - 0.15 s(t) / (100+1.85t)

Construct a direction field to see if there can be any conclusions about the behavior of solutions to the model:

In [15]:
VectorPlot[{1,10-(0.15s)/(100+1.85t)},{t,0,2},{s,0,10},InterpretationBox[DynamicModuleBox[{Set[Typeset`open, False]}, TemplateBox[{"List", StyleBox["\"Options\"", "IconizedCustomName", StripOnInput -> False], GridBox[{{ItemBox[RowBox[{TagBox["\"Head: \"", "IconizedLabel"], "\:f360", TagBox["List", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Length: \"", "IconizedLabel"], "\:f360", TagBox["1", "IconizedItem"]}]]}, {ItemBox[RowBox[{TagBox["\"Byte count: \"", "IconizedLabel"], "\:f360", TagBox["368", "IconizedItem"]}]]}}, GridBoxAlignment -> {"Columns" -> {{Left}}}, DefaultBaseStyle -> "Column", GridBoxItemSize -> {"Columns" -> {{Automatic}}, "Rows" -> {{Automatic}}}], Dynamic[Typeset`open]}, "IconizedObject"]], {FrameLabel -> {Style["t", 18], Style["s(t)", 18]}}, SelectWithContents -> True, Selectable -> False]]
Out[15]:
Output

Solutions tend to diverge to positive infinity. This behavior occurs because the model ignores the fact that the volume of the tank is finite.

Summary

Conclusions about the behavior of solutions to differential equations can largely be made from the direction field. Solutions that do not change over time are called equilibrium solutions. When solutions in a direction field tend to converge to the equilibrium solution, this describes a stable equilibrium point. The next lesson will cover the solutions to differential equations.