Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am trying to find a large system (>20) of coupled ordinary differential equations in order to approximate them numerically on the computer and check the efficiency and effectiveness of various steppers. Can someone point me to such a system? Thank you in advance

share|cite|improve this question
You can create arbitrarily large linear ODE systems of the form $y'=Ay+F(t)$ for an appropriate sized $A$. Alternatively, you can create non linear ODE systems by choosing the functions yourself. In this manner, you would have a benchmark to test the numerical error of your solver as you know the exact solution(s). – Daryl Feb 9 '13 at 23:23

If you don't care about the type of system (i.e. whether it's stiff or not) you could consider the [Heat equation][1]. $$\frac{\partial u }{\partial t} -\lambda \Delta^2 u = f$$ with $ \Delta u =\nabla^2 u = \sum_i \frac{\partial^2 u}{\partial x_i^2}$ and $\lambda > 0$.

Let us consider a 2D domain of length 1. We define $\Omega =[0,1]\times [0,1]$. If we want to approximate our domain with $n$ elements in each direction, we get $n^2$ elements in total. Therefore our equidistant stepsize $h=\frac{1}{n+1}$ if we exclude (count only inner points) the boundary. Let us assume $u=0$ on $\partial \Omega$, so we have homogenous boundary condition, and $u(t=0) = 0$ as our initial condition We approximate $\Delta u$ by $h^2\Delta u(i\cdot h, j\cdot h) \approx u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}$. With that we get: \begin{align} &u(t=0) = 0 \quad \text{as the initial condition}\\ &u(\partial \Omega) = 0 \quad \text{as our boundary condition}\\ &\frac{\partial}{\partial t} u_{i,j} - \frac{\lambda}{h^2} (u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}) = f_{ij}\\ \Leftrightarrow & \frac{\partial}{\partial t} u_{i,j} = \frac{\lambda}{h^2} (u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}) + f_{ij}\\ \end{align} Now we have to set up our grid. Let in the following $G$ denote the inner grid points, that means that $G$ does not include the boundary. We will enumerate our elements as follows: \begin{align} G = \begin{pmatrix} 1 & 2 & 3 & ... & n \\ n+1 & n+2 & n+3 &... & 2n\\ ... & ...& ... & ... &... &\\ ... & ...& ... & ... &n^2 \end{pmatrix} \end{align} This should present the numbering of $u$, so $u(2\cdot h, 3 \cdot h) = u_{n+3}$. Of course the same holds for $f$.

If we rewrite our ODE as $\tilde u_t =\frac{\lambda}{h^2} A \tilde u+F$ we get \begin{align} \tilde u &= (u_1,u_2,...,u_{n^2})^T\\ F &= (f_1,f_2,...,f_{n^2})^T \\ A &= \begin{pmatrix} \tilde{A} & I & 0 & 0 & ... & 0\\ I & \tilde{A} & I & 0 & ... & 0\\ 0 &I & \tilde{A} & I & ... & 0 \\ ... & ...& ...& ...& ...& ...\\ ... & ...& ...& ...& \tilde{A} & I \end{pmatrix} \in \mathbb{R}^{n^2 \times n^2} \end{align} with \begin{align} \tilde{A} &= \text{tridiag}(1,-4,1) \in \mathbb{R}^{n \times n}\\ I &= \text{diag(1)} \in \mathbb{R}^{n \times n}\\ 0 &\in \mathbb{R}^{n \times n} \end{align} This allows you to set up large systems easily, which has a physicall meaningfull background. For a simple example, you could choose $F$ to be zero everywhere except on one side of $\Omega$, heat should diffuse through the domain.

Note, that your system tends to get large very soon, as $A \in \mathbb{R}^{n^2\times n^2}$. Also note, that this is a sparse matrix, so there might be a method which has an advantage if it can handle sparse matrices.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.