How does a stiff equation solver work?

Question

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.

Answer 1

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!

Answer 2

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?)

Answer 3

Using ode15s

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

with

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.