Numerical Solutions and Direction Fields for Systems

Kurt Bryan and SIMIODE

A notebook to illustrate how to solve systems of ODEs numerically, and to draw a direction field for a pair of autonomous ODEs in Mathematica.

Example: Consider a pair of ODEs of the form x1'(t) = f(t,x1(t),x2(t)), x2'(t) = g(t,x1(t),x2(t)) where

f[x1_,x2_] = x1-x2^2
g[x1_,x2_] = x1*x2+x1

(in this example they are automous). And the ODEs themselves

de1 = x1'[t] == f[x1[t],x2[t]]
de2 = x2'[t] == g[x1[t],x2[t]]

Numerical Solution: We can solve the system numerically with data x1(0) = 4, x2(0) = 1 on the interval -2 < t < 2, with the command

sol = NDSolve[{de1,de2, x1[0]==4,x2[0]==1}, {x1,x2}, {t,-2,2}]

The solution at a given time t, say t = 1, can be obtained with

x1[1]/.sol
x2[1]/.sol

The resulting solution components can be plotted for -2 < t < 2 using the

plt1 = Plot[Evaluate[x1[t]/.sol],{t,-2,2},PlotRange->All,PlotStyle->{Red}];
plt2 = Plot[Evaluate[x2[t]/.sol],{t,-2,2},PlotRange->All,PlotStyle->{Blue}];
Show[plt1,plt2]

Or we can plot (x1(t), x2(t)) as a parametric curve:

ParametricPlot[Evaluate[{x1[t],x2[t]}]/.sol,{t,-2,2},PlotRange->All]

The numerical solution procedure of NDSolve works on systems of any dimension, as does the plotting of solution components xj(t) versus t. Parametric plotting of solution curves works in 2 and 3 dimensions.

Direction Fields: If the ODEs are autonomous, we can sketch a direction field on the range -5 < x1, x2 < 5 using the VectorPlot command, as

VectorPlot[{f[x1,x2],g[x1,x2]}, {x1,-5,5}, {x2,-5,5},VectorColorFunction->None,VectorStyle->Red]

Below is the same plot but with representative solution curves:

VectorPlot[{f[x1,x2],g[x1,x2]}, {x1,-5,5}, {x2,-5,5},VectorColorFunction->None,VectorStyle->Red,StreamPoints->Coarse,VectorColorFunction->None,StreamColorFunction->None,StreamStyle->Blue]