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.


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 Jan 10 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 Jan 10 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 Jan 10 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 Jan 11 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 Jan 11 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 Jan 11 at 17:12

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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