Desolve System Rk4

Description

We can also use the RK4 method to numerically solve systems of differential equations. Consider the predator-prey system

(1)
\begin{align} x' &= 2x - xy\\ y' &= -5y + xy \end{align}

As a nonlinear system, it is difficult, if not impossible, to find an analytical solution. However, we can use desolve_system_rk4 to approximate the solution curves. Since the command only outputs a list of points, we'll have to use a bit of Python to construct the plot ourselves. The arguments of the command are pretty much the same as those for desolve_system, except the command only takes the right-hand side of each equation and there are additional arguments end_points, which specifies the range over which we are approximating, and step, which specifies the step size for the RK4 method. A lower step size increases precision, but also computation time:good options are 0.1 and 0.01, depending on your approximation range.

In the first cell, we'll plot both our solution curves with respect to $t$:

Sage Cell

Code

var('x y t')
dx = 2*x - x*y
dy = -5*y + x*y
pts = desolve_system_rk4(des=[dx, dy], vars=[x, y], ivar=t, step = 0.01, end_points=[0, 10], ics=[0, 1, 1])
ptsx = [[i, j] for i, j, k in pts]
ptsy = [[i, k] for i, j, k in pts]
p = line(ptsx, color='blue') + line(ptsy, color='red')
p.show(axes_labels=['$t$', '$x, y$'])

In this cell, we plot the solution as a parametric curve in the $xy$ plane:

Code

var('x y t')
dx = 2*x - x*y
dy = -5*y + x*y
pts = desolve_system_rk4(des=[dx, dy], vars=[x, y], ivar=t, step = 0.01, end_points=[0, 10], ics=[0, 1, 1])
ptspar = [[j, k] for i, j, k in pts]
p = line(ptspar, color='blue')
p.show(axes_labels=['$x$', '$y$'], xmin=0, ymin=0)

In this cell, we do the same as in the previous cell, but superimpose the curve over the vector field of the system:

Code

var('x y t')
dx = 2*x - x*y
dy = -5*y + x*y
mag = sqrt(dx^2 + dy^2)
vf = vector([dx, dy])
pts = desolve_system_rk4(des=[dx, dy], vars=[x, y], ivar=t, step = 0.01, end_points=[0, 10], ics=[0, 1, 1])
ptspar = [[j, k] for i, j, k in pts]
p = line(ptspar, color='blue') + plot_vector_field(vf/mag, (x, 0, 15), (y, 0, 10))
p.show(axes_labels=['$x$', '$y$'], xmin=0, ymin=0)

Options

none

Tags

Primary Tags:

Secondary Tags:

A list of possible tags can be found at The WeBWorK Open Problem Library. For linear algebra tags see the Curated Courses Project.

Related Cells

Any related cells go here. Provide a link to the page containing the information about the cell.

Attribute

Permalink:

Author: Thomas Judson

Date: 14 May 2020 05:06

Submitted by: Zane Corbiere

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License