I am trying to understand how stiff differential equations are solved.

For instance the equation,

$$\frac{\partial y}{\partial t} = \alpha\frac{\partial ^2 y}{\partial z^2}$$

can be solved using an ode solver by discretizing z-direction at i th node

$$\frac{dy_i}{dt} = \alpha\frac{y_{i+1}-2y_i + y_{i-1}}{\Delta z^2}$$

This ode is often solved using ode15s. From what I understand, stiff solver uses numerical differentiation formals like backward differentiation formula(BDF) to approximate the derivative. Considering a two-step BDF formula,

$$\dot y(t_n)= \frac{\frac{3}{2}\cdot y_n + (-2)\cdot y_{n-1} + \frac{1}{2}y_{n-2}}{\Delta t }$$.

I am not sure how to proceed from here. Instead of using the solver, I am trying to implement the above BDF to understand how the solver works.

I would like to ask for suggestions on how to proceed from here.

EDIT: Equating (2) and (3) equation $$ \frac{\frac{3}{2}\cdot y^n_i + (-2)\cdot y^{n-1}_i + \frac{1}{2}y^{n-2}_i}{\Delta t }= \alpha\frac{y^n_{i+1}-2y^n_i + y^n_{i-1}}{\Delta z^2}$$ Here, n represents time and i respresents space.


3 Answers 3


The way you have set the equation up in your edit, allows you to solve it with thomas algorithm. If you look up an algorithm for TDMA, you can see that this equation will create a tridiagonal matrix A of coefficients of your (n) timestep, and a right hand side b depending on your previous timesteps (n-1, n-2) solution values. And this will allow for an easy solve of your BDF1 timestepping. The wiki page for thomas algorithm is here: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

In case you would like an explanation for how BDF1 works, you can think of it like so, I'll explain it for a steady state euler CFD code. In steady state euler we have our vector of conservative variables, $u$, and our residual operator corresponding to the discretized and integrated PDE, which must equal zero. $$R(u) = 0$$

We begin with $u_n$, for which $R(u_n)$ does not equal zero, and a $R(u_{n+1})$ which we want to be equal to zero. If we take the taylor series expansion of $R(u_{n+1})$ about $u_n$ we get: $$R(u_{n+1}) = R(u_n) + \frac{\partial R}{\partial u_n}\Delta u = 0$$ We can then bring the left term to the RHS and get: $$\frac{\partial R}{\partial u_n}\Delta u = -R(u_n)$$ Now, by solving this linear system, we will not actually solve the nonlinear governing equations in one step, because a linear taylor series expansion is not sufficient to solve the nonlinear equation that quickly, but we can repeat this iteration many times to solve the steady state problem, and you may notice this bears more that a passing resemblance to Newton's Method. Hope this was helpful!


The newest time is supposed to be the unknown, so the implicit method can be reshaped into a root finding equation system by moving the lhs to the rhs.

It seems you discretize y wrong though, the spatial discretization needs to be done separately from the temporal discretization, so the answer is a matrix if you keep the whole history. (Note: you have only shown the expression for the time derivative, not an expression for an integral or the taylor expansion, and shouldn't the lhs for the spatial stencil be the spatial derivative?)

  • $\begingroup$ Sorry if my original post wasn't clear.Please check the updated post. Spatial discretization is done separately and the time derivative is in the second last expression of the original post.Since the expression is in implicit form ,I am not clear about how to go about solving it. Would it be right to rearrange in a matrix form Ay = b and solve for the unknowns? $\endgroup$
    – Natasha
    Commented Jan 10, 2019 at 8:46
  • $\begingroup$ I would need to try at home and come back to you, I am relearning this myself. When I tried to solve an ODE yesterday it was nonlinear so I am not sure what the case is here. I got the implicit method on a form $F(t_0+dx,y(t_0+dx),t_0,y(t_0))=0$ or something $\endgroup$
    – Emil
    Commented Jan 10, 2019 at 16:44
  • $\begingroup$ My initial idea was to solve the spatial part separately and then use those values to step in time, that feels kind of familiar to me, perhaps you need to add boundary/initial(+time derivatives, robin?) conditions to know for sure how one should do it $\endgroup$
    – Emil
    Commented Jan 10, 2019 at 16:59
  • $\begingroup$ I'm using Dirichlet boundary condition at the inlet and Neumann BC at the exit. I am not clearly able to understand whether $F(t_0+dx,y(t_0+dx),t_0,y(t_0))=0$ form has to be solved using a non-linear equation or whether a matrix equation can be solved.I'm going through a reference in this link in the section that explains Crank-Nicolson scheme Ax=b is arrived in the last step. $\endgroup$
    – Natasha
    Commented Jan 11, 2019 at 13:02
  • $\begingroup$ @Natasha: my point was that you get some equation that you want to find the kernel of, in the case of $Ax=b$ it would be something like $F(x)=Ax-b=0$, I think. (Do you have some equation with known solutions to try on?) $\endgroup$
    – Emil
    Commented Jan 11, 2019 at 19:45

Using ode15s

[t,y] = ode15s( @(t,z) fun(t,z,alpha,dz),tspan,y0);


function dy = fun(t,y,alpha,dz)
dy = zeros(size(y));
dy(1) =  % b.c.
dy(2:end-1) = alpha*(y(3:end) -2*y(2:end-1) + y(1:end-2) )/dz^2
dy(end) = % b.c.

where dy(1) and dy(end) depend of the boundary conditions.

If you want to implement your own numeric method, you can search of any method to solve system of linear ODEs.

  • $\begingroup$ I think OP wants to implement the solver by him/herself - not with the help of a black-box function like ode15s $\endgroup$
    – GoHokies
    Commented Jan 11, 2019 at 17:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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