1532-Numerically Solving Systems of Differential Equations.sagews

Numerically Solving Systems of Differential Equations

Later in the semester, we will take a closer look at numerical solvers for systems of differential equations.  As it will turn out, Euler's method does work for systems, but not very well for many systems with a periodic solution (like the model of an ideal pendulum).  So, for now, we will treat numerical solvers as a "black box": inputs go in, completely undetectable stuff goes on, and then outputs come out.  The following notebook pulls together some of the SAGE package's tools for solving and plotting numerical solutions for 2-dimensional systems.

So, let's suppose that you want to solve$\frac{dY}{dt}=\left[\begin{array}{c}f(t,x,y)\\ g(t,x,y)\end{array}\right],\ \ \ \ Y(t_0)=\left[\begin{array}{c}x_0\\ y_0\end{array}\right]$You can type:

solution=SolveSystem([f(t,x,y),g(t,x,y)],[t,t0,t1],[x0,y0])

and SAGE will produce a numerical solution, stored in the variable solution.  For example, to solve the system$\frac{dY}{dt}=\left[\begin{array}{c}-x+2y\\ -2y\end{array}\right]\ \ \ \ Y(0)=\left[\begin{array}{c}4\\ 5\end{array}\right]$enter

solution=SolveSystem([-x+2*y,-2*y],[t,0,10],[4,5])

At this point, SAGE has stored a numerical solution as a list of triples $(t,x,y)$.  We will probably never want to look at the coordinates.  Instead, we would like to see some plots of the solution.  You can produce two plots easily: the phase plane portrait (the graph of the solution $(x(t),y(t))$ as a parametric curve in the $(x,y)$ plane), and the coordinate plots.

solution.PhasePlane()

solution.Coordinates()

The x coordinate is shown in blue, the y coordinate in red.

We would frequently like to see the "flow" of the differential equation.  One way to do this is to randomly choose a bunch of initial points in a specified rectangle, solve the system of differential equations starting at each of the initial conditions, and then plot all of the solutions you just computed.  The syntax for this is similar to the SolveSystem command:

solution=SolveFlow([f(t,x,y),g(t,x,y)],[t,t0,t1],[[a,b],[c,d]])

where $[[a,b],[c,d]]$ describes the rectangle from which we will draw the initial conditions.  As an example, let's plot the flow for the system we just solved.

solution=SolveFlow([-x+2*y,-2*y],[t,0,10],[[3.5,4.5],[4.5,5.5]])

This will very likely take a little longer than the first example did.  By default, we are choosing 20 initial conditions and solving the differential equation for each.  So, we should expect the approximation process to take about twenty times as long.  (The number of randomly chosen initial conditions can be changed; see the section on options, below.  For more initial conditions, more time is required.)

Now we plot the solutions in the phase plane and in coordinate space as we did before.

solution.PhasePlane()

solution.Coordinates()

The x coordinate is shown in blue, the y coordinate in red.

Setting Options for the Solver

At this time, you can tweak only a few options in the solver.  Typing

solution=SolveSystem([f(t,x,y),g(t,x,y)],[t,t0,t1],[x0,y0],stepsize=decimal number)

will change the numerical step size the solver uses.  (This is the same as the step size in Euler's Method.)  By default, the stepsize is 0.05.  Decreasing it will make solutions more accurate, but will require more computing time to solve.

For computing flows, two options may be set:

solution=SolveFlow([f(t,x,y),g(t,x,y)],[t,t0,t1],[[a,b],[c,d]],stepsize=decimal number,flows=positive integer)

The stepsize option is as above.  The flows option determines the number of initial conditions to randomly choose in the specified rectangle $[[a,b],[c,d]]$.  Increasing the number of flows to draw can result in a more complete picture, but requires more computing time.  Setting this number above 100 is probably not helpful, since it will require a large amount of computing time and will clutter up the image with too much noise.  By default, flows is set to 20.

Warnings and Notes

First, always use $x$, $y$, and $t$ as your variables when defining your differential equation.  This restriction was made to simplify the command syntax.

Second, you can save the plots you make for printing or inclusion in homework writeups.  Just right-click on the plot and select "Save image as..." or the corresponding choice in your web browser.  (That is, treat the picture just like any picture in any web page.)

Third, the code at the bottom of this page is necessary for solving differential equations with this particular syntax.  Don't delete it!

Code:

#auto
from sage.calculus.desolvers import desolve_system_rk4
x,y,t=var('x y t')
class PlanarDE:
    def __init__(self,point_list):
        self.T=[i for (i,j,k) in point_list]
        self.X=[j for (i,j,k) in point_list]
        self.Y=[k for (i,j,k) in point_list]
    def PhasePlane(self,color='black'):
        return line(zip(self.X,self.Y),color=color)
    def Coordinates(self,info=True):
        if info: print "The x coordinate is shown in blue, the y coordinate in red."
        return line(zip(self.T,self.X),color='blue')+line(zip(self.T,self.Y),color='red')
class PlanarFlow():
    def __init__(self,de_list):
        self.des=de_list
    def PhasePlane(self):
        colorchoices=['black','red','blue','green','yellow','orange','brown']
        colors=[choice(colorchoices) for i in range(len(self.des))]
        return sum([z.PhasePlane(color=colors[i]) for i,z in enumerate(self.des)])
    def Coordinates(self,info=True):
        if info: print "The x coordinate is shown in blue, the y coordinate in red."
        return sum([z.Coordinates(info=False) for z in self.des])    
def SolveSystem(dynamics,TInterval,ICs,stepsize=0.05):
    try:
        soln=PlanarDE(desolve_system_rk4(dynamics,[x,y],ics=[TInterval[1],ICs[0],ICs[1]], \
        ivar=TInterval[0],end_points=TInterval[2],step=stepsize))        
    except:
        print "Oops!  The solver failed to numerically solve your differential equation."
        soln=[]
    return soln
def SolveFlow(dynamics,TInterval,rect,stepsize=0.05,flows=20):
    flowsfix=max(1,int(flows))
    a,b=rect[0];c,d=rect[1]
    def r(): return (a+(b-a)*random(),c+(d-c)*random())
    icdata=[r() for i in range(flowsfix)]
    return PlanarFlow([SolveSystem(dynamics,TInterval,z,stepsize=stepsize) for z in icdata])

Modeling with Systems (Section 2.1)

In this section we look at two models: a predator-prey model and an oscillating spring system.  The predator-prey model is a classic example of constructing more complex behaviors from simple models.

The Predator-Prey Example:

Begin with the two exponential growth models for the prey ($x$) and predator ($y$):$\frac{dx}{dt}=2x\ \ \ \frac{dy}{dt}=-y$We then add an interaction term to "couple" the populations together.$\frac{dY}{dt}=\left[\begin{array}{c}2x-1.2xy\\ -y+0.9xy\end{array}\right]$Let's take a look at the phase plane and coordinate plots.

F=[2*x-1.2*x*y,-y+0.9*x*y]
solution=SolveSystem(F,[t,0,20],[1,0.5])

show(solution.PhasePlane())

show(solution.Coordinates())

The x coordinate is shown in blue, the y coordinate in red.

We then modify the system to model logistic growth in the prey.$\frac{dY}{dt}=\left[\begin{array}{c}2x\left(1-\frac{x}{2}\right)-1.2xy\\ -y+0.9xy\end{array}\right]$

G=[2*x*(1-x/2)-1.2*x*y,-y+0.9*x*y]
solution2=SolveSystem(G,[t,0,20],[1,0.5])

show(solution2.PhasePlane())

show(solution2.Coordinates())

The x coordinate is shown in blue, the y coordinate in red.

The Mass-Spring Example:

Let $x$ denote position and $y$ denote velocity.  Then, we rewrite the second-order equation$\frac{d^2x}{dt^2}=-\frac{k}{m}x$as the first-order system$\frac{dY}{dt}=\left[\begin{array}{c}y\\ -\frac{k}{m}x\end{array}\right]$

H=[y,-x]
solution3=SolveSystem(H,[t,0,15.],[0.,1.])

plot1=solution3.PhasePlane()
show(plot1,aspect_ratio=1)

show(solution3.Coordinates())

The x coordinate is shown in blue, the y coordinate in red.

solution4=SolveSystem(H,[t,0,15.],[0.5,1.])

plot2=solution4.PhasePlane()
show(plot1+plot2,aspect_ratio=1)

solution5=SolveSystem([2*x-3*y,-8*x+12*y],[t,0,5],[2,5])

show(solution5.PhasePlane())

Vector fields and solutions

Plotting a vector field in SAGE is straightforward.  For example, we can plot a vector field for the system$\frac{dY}{dt}=\left[\begin{array}{c}2x-1.2xy\\ -y+0.9xy\end{array}\right]\ \ \ \ Y(0)=\left[\begin{array}{c}1\\ 0.5\end{array}\right]

p1=plot_vector_field([2*x-1.2*x*y,-y+0.9*x*y],(x,0,5),(y,0,5))
show(p1)

We can add a solution curve (or more!) on top of this.

soln=SolveSystem([2*x-1.2*x*y,-y+0.9*x*y],(t,0,10),[1,0.5])
soln2=SolveSystem([2*x-1.2*x*y,-y+0.9*x*y],(t,0,10),[4,3])

p2=soln.PhasePlane(color='red')
p3=soln2.PhasePlane(color='green')
show(p1+p2+p3)

Van der Pol Equation

F=[y,-x+(1-x^2)*y]

solution=SolveSystem(F,[t,0,100],[1,1])
solution2=SolveSystem(F,[t,0,100],[2.5,1.6])

vec_field=plot_vector_field(F,(x,-3,3),(y,-3,3))
show(solution.PhasePlane(color='blue')+solution2.PhasePlane(color='red')+vec_field,aspect_ratio=1)

p1=plot_vector_field([2*x+2*y,-4*x+6*y],(x,-3,3),(y,-3,3))
show(p1)

soln=SolveSystem([2*x+2*y,-4*x+6*y],(t,0,1.05),[0.05,0.1])
p2=soln.PhasePlane(color='red')
show(p1+p2)