Numerical Solution to Ordinary Differential Equations

Ordinary differential equations (ODEs) are commonly encountered in scientific and engineering analysis. The simplest form of ODEs is the 1st order ODEs:

y' = f(y, x)

Here, y(x) is the function dependent on x, and it is the function that we want to solve for. X is the independent variable. The left hand side of the equation is the first order derivative of the function y, while the right hand side of the equation is a function that links the derivative to some combination of y and x. Function y is also known in ODEs as the solution function, while function f(y, x) is known as the flow function.

In generally, ODEs can be more than 1st order. It can be readily shown that, any ODEs (of any order) can be rewritten as a system of 1st order ODEs. Take the motion equation for the harmonic oscillator as an example:

x(t)'' = -k/m*x(t)

This equation can be rewritten as two 1st order ODEs:

 x' = F1(x,v,t) =v

v' = F2(x,v,t)=-k/m*x

Here, the harmonic oscillator 2nd order ODE is rewritten as a system of two 1st order ODEs, where there are two solution variables, x and v, two flow functions, F1 and F2, and one independent variable, t.

Numerical solutions to ODEs are readily achievable through discreet integration. The idea is that given an initial condition, i.e. y(x=x0), one can approximately compute the value of y at the next step, i.e. y(x=x0+h), assuming h is sufficiently small. The most straightforward solution is the Euler method. Takes the simplest, 1st order ODE as an example, at any given step, one can approximate the next step as:

y(n+1) = y(n) +h*f(y(n),x(n))

The essence of the Euler equation shown above is that, the next y value can be approximated by the current y value y(n), plus the step size times the derivative (slope) of the current x location, x(n). This is quite a crude estimation, because the derivative of y could change as x changes from x(n) to x(n+1). However, as a first order approximation, sometimes the Euler method is good enough, and its advantage is that it's easy to understand.

The Euler method can be easily extended to solve any ODEs with more than 1st order. As shown earlier, any ODE can be rewritten as a system of 1st order ODEs:

Y' = F(Y, x)

Here, vector notion is used to simplified the presentation of an ODE system. Y is a vector of all the solution variable in the system, while F is the vector of all the flow functions in the system. Therefore, Y is also known as the solution vector, while F is known as the flow vector. For example, in the harmonic oscillator example, Y=(x,v), and F=(F1(x, v, t), F2(x, v, t)).      

Applying the Euler method to an 1st ODE system is very similar to a single 1st order ODE:

Y(n+1) = Y(n) + h*F(Y(n), x(n))

Here, Y and F are both vector, while, h and x are scalar. The equation above in fact represents a set of equations. The vector notation is just to simplify the presentation of them. 

The disadvantage of the Euler method, despite its simplicity and intuitiveness, is that it's quite unstable, and tends to produce large error. In the next blog entry, we will introduce another method, the Runge-Kutter method, which is the work house for most of numerical analysis of ODEs.

Block comment and uncomment codes in VIM

In VIM, it is quite easy to comment/uncomment a block of lines.

To comment a block of lines:
a. Ctrl-V and mark the block of lines to be commented out;
b. shift-I and write the text that you want to prepend to each line, e.g. #
c. press escape.

To uncomment a block of lines:
a. Ctrl-V and go down vertically until the last commented line;
b. press x.