This is a short tutorial on how to use ode_solver in SageMath. For a more detailed reference, please read Solving ODE numerically by GSL In order to use SageMath, you can either install it on your computer, or you can use it on CoCalc.

ode_solver is a Python class, so in order to use it, we need to make an instance of it:

In [1]:
T = ode_solver()

The basic routine is then:

  1. set T.function to a function that describes the differential equation system.
  2. call T.ode_solve(y_0, t_span, num_points) to solve the system numerically, where y_0 is a list of initial values at the point t = t_span[0], t_span is the interval on which we would like to solve the system, and num_points is the number of points we want to compute in the interval t_span. In general, more points take longer to compute but give more precise results.
  3. call T.plot_solution() to plot the solution. If needed, we can access the computational results in T.solution, which is a list of tuples of the form [(t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])].

First Order Equation

Let's say we would like to numerically solve

$$\left\{\begin{align*} & t^2y' = y + 3t \\ & y(1) = -2 \end{align*}\right.$$

over $1 \le t \le 4$.

We set $y_0 = y$, then $$y_0^\prime = \frac{y_0}{t^2} + \frac{3}{t}.$$

We try to use $y_0$ instead of $y$ in order to make T.function clearer.

In [2]:
def f1(t, y):
    return [y[0]/(t^2) + 3/t]

f1 is a function to describe the system. It takes 2 arguments t and y, where t is the independent variable and y is a list of variables. In general, for an $n$-th order differential equation, y should be a list of $n$ elements $[y_0, y_1, \dots, y_{n-1}]$, where $y_i$ is representing the $i$-th derivative $y^{(i)}$ of $y$. The return of f1 is a list of derivatives $[y_0^\prime, y_1^\prime, \dots, y_{n-1}^\prime]$, computed using $t$ and $[y_0, y_1, \dots, y_{n-1}]$. This explains why we use $y_0$ in place of $y$ above.

In [3]:
T.function = f1
T.ode_solve(y_0=[-2], t_span=[1, 4], num_points=100)
T.plot_solution()
Out[3]:

If we want to know the computed valued at $t = 4$, we can access T.solution to look at its last element.

In [4]:
T.solution[-1]
Out[4]:
(3.999999999999988, [1.4613570428866713])

Sometime, we don't get exact $t$ value because of the way how computers represent float numbers. We can see that $y(4) \approx 1.46$.

Second Order Equation

For example, we try to solve numerically $$\left\{\begin{align*} & y^{\prime\prime} + e^t y^\prime + ty = 3 \sin(2t) \\ & y(0) = 2, y^\prime(0) = 8 \end{align*}\right.$$ over $0 \le t \le 4$.

We set $y_0 = y$ and $y_1 = y^\prime$, then the above differential equation is equivalent to the system $$\left\{\begin{align*} & y_0^\prime = y_1 \\ & y_1^\prime = -ty_0 - e^t y_1 + 3\sin(2t) \\ & y_0(0) = 2, y_1(0) = 8 \end{align*}\right.$$

Rewriting the differential equation in the above equivalent form is necessary for us to define T.function correctly.

In [5]:
def f2(t, y):
    return [y[1], -t*y[0]-exp(t)*y[1] + 3*sin(2*t)]

T.function = f2
T.ode_solve(y_0=[2, 8], t_span=[0, 4], num_points=100)
T.plot_solution()
Out[5]:
In [6]:
T.solution[-1]
Out[6]:
(4.000000000000003, [3.641260138622823, -0.21523301500916253])

At $t = 4$, we have $y(4) = y_0(4) \approx 3.64$ and $y^\prime(4) = y_1(4) \approx -0.22$.

Differential Equation with Parameters

Sometime when we would like to study a differential equation with parameters and try to understand the behavior of the equation when the parameters vary. There is a way in SageMath to define T.function with parameters. Let us look at the following example: $$\left\{\begin{align*} & y^\prime + \epsilon y^2 = t \\ & y(0) = 0 \end{align*}\right.$$ over $0 \le t \le 10$. Analyze $y(10)$ of the equation as $\epsilon \to \infty$.

In [7]:
def f3(t, y, params):
    return [t - params[0] * (y[0]^2)]

We notice an extra argument params in the definition of f above. It is a list of parameters in the defferential equation. If we want to set $\epsilon = 0.2$ in the equation, we just need to pass a list params to T.ode_solve as follows.

In [8]:
T.function = f3
T.ode_solve(y_0=[0], t_span=[0, 10], num_points=100, params=[0.2])
T.plot_solution()
Out[8]:

Since we want to analyze $y(10)$ as $\epsilon$ goes to $\infty$, we can try $\epsilon = 1, 10, 100, 1000, 10000, \dots$ and list the corresponding $y(10)$.

In [9]:
T.function = f3
print 'epsilon \t y(10)'
for i in range(5):
    eps = 10^i
    T.ode_solve(y_0=[0], t_span=[0, 10], num_points=100, params=[eps])
    print eps, '\t\t', T.solution[-1][1][0]

epsilon 	 y(10)
1 		3.13675822167
10 		0.997484135111
100 		0.315977269577
1000 		0.099974984368
10000 		0.0316202761155

We might conclude from this experiment that $y(10) \to 0$ as $\epsilon \to \infty$.