Instruction of ode.py

---

Name: Feng Chen
Student number: 2013301020145
March 24, 2016

Abstract

This is the instruction file for the first order ordinary differential equation algorithms package ode.py. There are totally four kinds of numerical algorithms in this package temporarily, including simple Euler method, modified Euler method, Second-order Runge-Kutta method, Third-order Runge-Kutta method.
On the other hand, this package can be also used to solve differential equation systems, which I will talk about in the latter part of insturction.

Numerical Algorithms Introduction

For initial value problem of first-order ordinary differential equation:

$$\begin{equation} \begin{cases} \frac{dy}{dx}=f(x,y) & a\leq x \leq b \\ y(a) = y_0 \end{cases} \end{equation}$$ where a and b are real constant.
Therefore, when we make a Taylor expression of y at position x, we will get:
$$\begin{equation} y(x+\Delta x)\approx y(x)+\frac{dy}{dx}\Delta x=y(x)+f(x,y)\Delta x \end{equation}$$ Changing the form of this formula, we can get:
$$\begin{equation} y_{n+1} = y_n +f(x_n,y_n)\Delta x \end{equation}$$ This is the recursion equation of simple Euler method. If we try to use different recursion equation, which includes the higher order terms' contribution, we will get a more accurate result than simple Euler method.
Here are some useful recursion equations I used in this programme.

1. Simple Euler method

$$\begin{equation} y_{n+1}=y_n+\Delta xf(x_n,y_n) \end{equation}$$

2. Modified Euler method

$$\begin{equation} \begin{cases} t=y_n+\Delta xf(x_n,y_n) \\ y_{n+1}=y_n+\frac{\Delta x}{2}[f(x_n,y_n)+f(x_{n+1},t)] \end{cases} \end{equation}$$

3. Second-order Runge-Kutta method

$$\begin{equation} y_{n+1}=y_n+\Delta x[f(x_n+\frac{\Delta x}{2},y_n+\frac{\Delta x}{2}f(x_n,y_n))] \end{equation}$$

4. Third-order Runge-Kutta method

$$\begin{equation} \begin{cases} K_1=f(x_n,y_n) \\ K_2=f(x_n+\frac{\Delta x}{2},y_n+\frac{\Delta x}{2}K_1)\\ K_3=f(x_n+\Delta x,y_n-\Delta xK_1+2\Delta xK_2)\\ y_{n+1}=y_n+\frac{\Delta x}{6}(K_1+4K_2+K_3) \end{cases} \end{equation}$$

Usage and Examples

Firstly, we should set the function on the right side and the name of independent variable and dependent variable.

   ode(step,left limit,right limit,initial condition)
   set_fx(f_expression,independet_variable,dependent_variable) 

For example:

import ode
A = ode.ode(0.05,0,1,10)
ode.set_fx('y','x','y')

Then we can use different algorithm to solve the problem.

ode.euler()      #### simple euler method
ode.eulerplus()  #### modified euler method
ode.rgkt_2()     #### 2nd Runge-Kutta method
ode.rgkt_3()     #### 3rd Runge-Kutta method

Example:

import ode as solving_ode
import matplotlib.pyplot as plt
A = solving_ode.ode(0.1,0,2,5)  ##### setting the initial population N(0) = 10
A.set_fx('N',['t','N'])
eulerplus_record = A.euler()
rgkt_3_record = A.rgkt_3()
plt.figure(figsize = (10,6))
plt.plot(eulerplus_record[0],eulerplus_record[1][0],'*',label = 'simple euler method')
plt.plot(rgkt_3_record[0],rgkt_3_record[1][0],label = '3rd Runge-Kutta method')
plt.legend()
plt.savefig('chapter1_1.6.png',dpi = 144)
plt.show()

Algorithms Testing

Solving equation: $y'=y$, setting the initial value: $y_0=5$, step:$0.1$

The bigger, the faster
As the figures showed above, though Simple Euler method is the fastest algorithm, it makes more mistakes than other three kinds of algorithms.
It shows that Modified Euler method is a very good approach to solve ordinary differential equation. You can choose an appropriate method to solve your own problem.

Solving differential equation systems or Higher order differential equation

Equation systems:

Giving a set of equation systems：

$$\begin{equation} \begin{cases} \frac{dy_1}{dx}=f_1(x,y_1,y_2,...,y_n) \\ \frac{dy_2}{dx}=f_2(x,y_1,y_2,...,y_n) \\ .\\ .\\ .\\ \frac{dy_n}{dx}=f_n(x,y_1,y_2,...,y_n) \end{cases} \end{equation}$$
What we should do is to make Taylor expension on each independent variable $y_i$, then we will get:
$$\begin{equation} \begin{cases} y_1(x+\Delta x)\approx y_1(x)+f_1(x,y_1,y_2,...,y_n)\Delta x \\ y_2(x+\Delta x)\approx y_2(x)+f_2(x,y_1,y_2,...,y_n)\Delta x \\ .\\ .\\ .\\ y_n(x+\Delta x)\approx y_n(x)+f_n(x,y_1,y_2,...,y_n)\Delta x \end{cases} \end{equation}$$ which means:
$$\begin{equation} \begin{cases} y_{1(m+1)}=y_{1(m)}+f_1(x_m,y_{1(m)},y_{2(m)},...,y_{n(m)})\Delta x \\ y_{2(m+1)}=y_{2(m)}+f_2(x_m,y_{1(m)},y_{2(m)},...,y_{n(m)})\Delta x \\ .\\ .\\ .\\ y_{n(m+1)}=y_{n(m)}+f_n(x_m,y_{1(m)},y_{2(m)},...,y_{n(m)})\Delta x \end{cases} \end{equation}$$ These are recursion equations for Simple Euler Method.
Similarly, we can use other algorithms to solve this problem by improving these recursion equations just like what I showed in Algorithms Introduction

Higher-order equation:

Giving a n-order differential equation:

$$\begin{equation} y^{(n)}+p_1(x)y^{(n-1)}+...+p_{n-1}(x)y'+p_n(x)y = f(x) \end{equation}$$ Setting:
$$\begin{equation} y_1(x)=y', y_2(x)=y'',...,y_{n-1}(x)=y^{(n-1)} \end{equation}$$ Then the differential equation will be transferred into a set of equations:
$$\begin{equation} \begin{cases} \frac{dy}{dx}= y_1(x) \\ \frac{dy_1}{dx}=y_2(x) \\ .\\ .\\ .\\ \frac{dy_{n-1}}{dx} = -p_1(x)y_{n-1}-...-p_{n-1}(x)y_1-p_n(x)y + f(x) \end{cases} \end{equation}$$ Therefore, we can use the same methods as equation systems to solve this problem.

Example: solve $y''=-xy+e^xy'+3sin2x$

Firstly, simpilify this equation into:

$$\begin{equation} \begin{cases} \frac{dy}{dx}= y_1(x) \\ \frac{dy_1}{dx} =-xy+e^xy_1+3sin2x \end{cases} \end{equation}$$ where $3.8\leq x \leq 4.0, y(0)=2, y_1(0)=8$

import ode as solving_ode
import matplotlib.pyplot as plt
import math
A = solving_ode.ode(0.005,3.8,4,(2,8))  
A.set_fx(('y1','-x*y+y1*math.e**x+3*math.sin(2*x)'),['x','y','y1'])
rgkt_3_record = A.rgkt_3()
B = solving_ode.ode(0.005,3.8,4,(2,8))
B.set_fx(('y1','-x*y+y1*math.e**x+3*math.sin(2*x)'),['x','y','y1'])
euler_record = B.euler()
plt.figure(figsize = (10,6))
plt.plot(euler_record[0],euler_record[1][0],label = 'simple euler method')
plt.plot(rgkt_3_record[0],rgkt_3_record[1][0],label = '3rd Runge-Kutta method')
plt.legend(loc = 'left upper')
plt.savefig('example2.png',dpi = 144)
plt.show()

Package code

• ode.py
• example.py
• example2.py