I am trying to convert the diffusion equation to ODEs so that it can be programmed using Matlab's ODE solvers. The diffusion equation I'm using is:
$$
{\partial u \over \partial t} = D\,{\partial^{2}u \over \partial x^{2}}\quad
\mbox{where}\ D\ \mbox{is diffusivity}
$$
I am not a mathematician so have become confused by the mathematical notation used without the explanation to go along side it. I have been pointed in the direction of the method of lines to help solve it but haven't managed to understand it.
I have tried programming it using pdepe, but one of my initial boundary conditions needs to be different to the other and I can't get this to happen using it.
Any pointers or explanations on how to do this would be appreciated.
Best Answer
You may use the method of lines (MOL) in two different manners: by discretizing either the time variable, so you obtain a boundary value problem for each time step (for example applying a $\theta$-method discretization):
$$\frac{U^{n+1}(x) - U^{n}(x)}{\Delta t} = D \left( \theta \frac{d^2 U^{n+1}}{d x^2} + (1- \theta)\frac{d^2 U^n}{d x^2} \right), \quad 0 \leq \theta \leq 1, $$ which is to be solved in each time step with two boundary conditions for $U^n(x)$ (the solution in the step $t_n = t_0 + n \Delta t$) and specifying an initial value for $U^0(x)$; or by discretizing the spatial variable, fixing $t$:
$$ \frac{d U_i(t)}{d t} = \frac{D}{\Delta x^2} (U_{i+1} (t) - 2 U_i(t) + U_{i-1} (t)),$$ which is to be solved with a given initial condition $U_i(0)$ (the solution at the grid point $x_i = x_0 + i \Delta x$) and applying the boundary conditions at both end of your domain. This last form of expressing MOL can me more suitable to implement in Matlab, as it can be rearranged in matrix form as a linear system:
$$\vec{U}'(t) = A \vec{U}(t), \quad \vec{U} = (U_0, \ldots, U_N),$$
and then apply, for example, ode23 to solve for every unknown $U_i$. The tridiagonal matrix $A$ is given by the coefficients of the second order discretization:
$$A = \frac{D}{\Delta x^2}<a_{k,j-1}, a_{kk}, a_{k,j+1}>= \frac{D}{\Delta x^2}<1, -2, 1>.$$
However, the $a_{11}, a_{12}$ and $a_{N,N-1}$, $a_{NN}$ terms depend on whether the boundary conditios are either Dirichlet's or Neumann's, as well as the size of $A$ and unknowns $U_i$. I hope this may be useful to you.
Cheers!