# Transforming an implicit difference scheme into an explicit one

Given is the following ordinary differential equation:

$\epsilon y''+y' = 0$,

$y(0) = 0$,

$y(1) = 1$

How can one transform the implicit difference scheme $$\epsilon \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^2} + \frac{y_{i+1}-y_{i-1}}{2 h} = 0, \quad i = 1, \dots, n-1\\ y_0 = 0, \qquad y_n = 1$$ into an explicit one?

My idea:

Transforming the difference scheme above into the form $A y = 0$ with a matrix $A$.

Thus, we would get

$A = \pmatrix{c_0 & 0 & \dots & 0 \\ \frac{\epsilon}{h^2} - \frac{1}{2 h} & \frac{2 \epsilon}{h^2} & \frac{\epsilon}{h^2} + \frac{1}{2 h} \\ \dots & \dots & \dots & \dots \\ 0 & 0 & -a_n & c_n}$

Unfortunately, I' ve got no idea how to continue and what to do with $c_0$ and $-a_n$, $c_n$.

Could anyone of you help me, please?

Thank you in advance,

Peter123.

• Your difference problem is missing boundary conditions. It should be apparent now what are the actual values of $c_0, c_n$ and $a_n$. – uranix Nov 23 '16 at 16:00

## 3 Answers

You just need to write the equation for $i=1$ and $i = n$ explicitly. Because I am not sure whether $\epsilon$ is constant, I will write your first equation as

$$a_i y_{i+1} + b_{i}y_i + c_i y_{i-1} = 0 \quad\mbox{for}\quad i = 1, \cdots, n - 1 \tag{1}$$

with

$$a_i = \frac{\epsilon(x_i)}{h^2} + \frac{1}{2h}, \quad c_i = \frac{\epsilon(x_i)}{h^2} - \frac{1}{2h}, \quad\mbox{and}\quad b_i = -\frac{2\epsilon(x_i)}{h^2} \tag{2}$$

From Eq. (1) we can evaluate

• $i = 1$:

$$a_1 y_2 + b_1y_1 + c_1 \underbrace{y_0}_{=0} = 0$$

• $i = n-1$:

$$a_{n-1} \underbrace{y_n}_{=1} + b_{n-1}y_{n-1} + c_{n-1}y_{n-2} = 0$$

The system of equations is then reduced by two $y_1, \cdots, y_{n-1}$

$$\left( \begin{array}{8} b_1 & a_1 & & & & & & \\ c_2 & b_2 & a_2 & & & & & \\ & c_3 & b_3 & a_3 & \cdots & & & \\ & & \vdots & \ddots & & & & \\ & & & & & c_{n-2} & b_{n-2} & a_{n-2} \\ & & & & & & c_{n-1} & b_{n-1} \\ \end{array} \right) \left( \begin{array}{c} y_2 \\ y_3 \\ \vdots \\ y_{n-2} \\ y_{n-1}\end{array}\right) = \left( \begin{array}{c} 0 \\ 0 \\ \vdots \\ 0 \\ -a_{n-1}\end{array}\right) \tag{3}$$

There are many efficient algorithms to solve this problem, here is an example using $\epsilon = 0.5$. Blue is the actual solution, red is solving Eq. (3) with 10 nodes. • Thank you very much, too! You also helped me a lot. Do you know what problems there are if we choose ϵ=hϵ=h for the step size regarding the convergence of yiyi for h ---> 0? – Peter123 Nov 23 '16 at 16:53
• If you make $\epsilon = h \to 0$ you run into another problem: the equation becomes $y' = 0$ with $y(0) = 0$ and $y(1) = 0$ which has no solution – caverac Nov 23 '16 at 17:00

Well, you can apply tridiagonal matrix algorithm to your matrix. The result may be viewed as an explicit scheme (forward pass and backward pass).

Second approach is to introduce new unknown $v_n$, defined as $$v_i = \frac{y_{i+1} - y_i}{h}, \quad i = 0, \dots, n-1$$ Then your equation can be rewritten as an explicit scheme $$\epsilon\frac{v_i - v_{i-1}}{h} + \frac{v_i + v_{i-1}}{2} = 0, \quad i = 1, \dots, n-1,$$ but at cost of changing the boundary conditions to the integral-like condition $$v_{n-1} + v_{n-2} + \dots + v_0 = \frac{y_n - y_0}{h} = \frac{1}{h}$$ Despite this fact the problem can be easily solved: $$v_i = \frac{2\epsilon - h}{2\epsilon + h} v_{i-1}\\ \frac{1}{h} = \sum_{i=0}^{n-1} v_i = v_0 \sum_{i=0}^{n-1} q^i = v_0 \frac{1 - q^n}{1 - q}, \quad q = \frac{2\epsilon - h}{2\epsilon + h}$$ Solving the condition for the $v_0$ we obtain $$v_0 = \frac{1-q}{h(1-q^n)}, \quad v_i = q^i\frac{1-q}{h(1-q^n)}$$ and $$y_{i+1} = y_i + h v_i, \quad i = 0, \dots, n-1\\ y_0 = 0$$

Third approach views the scheme a Cauchy problem for a second-order equation with unknown second condition. $$\epsilon\frac{y_{i+1} - 2y_i + y_{i-1}}{h^2} + \frac{y_{i+1} - y_{i-1}}{2h} = 0\\ y_0 = 0\\ y_1 = \beta.$$ The scheme is now explicit, but the correct value for $\beta$ is unknown yet. Due to linearity of the problem, every $y_i(\beta)$ is a linear function in $\beta$, so we could solve the problem for, say, $\beta = 0$ and $\beta = 1$. Let's denote the solution with $\beta = 0$ as $y_i^{(0)}$ and the solution with $\beta = 1$ as $y_i^{(1)}$. Then the solution for arbitrary $\beta$ is given by $$y_i^{(\beta)} = (1-\beta) y_i^{(0)} + \beta y_i^{(1)}.$$ Using the right boundary condition $y(1) = 1$ and solving it for $\beta$ $$y_n^{(\beta)} = (1-\beta) y_n^{(0)} + \beta y_n^{(1)} = 1$$ gives the correct value of $\beta$ along with the correct solution $y_i^{(\beta)}$.

• Thank you very much! That helps me a lot. Do you also know what problems there are if we choose $\epsilon = h$ for the step size regarding the convergence of $y_i$ for h ---> 0? – Peter123 Nov 23 '16 at 16:51
• Plug it in the second method and observe that $q = 1/3, y_1 \to 2/3, v_0 \to \infty$ so the limiting solution is discontinuous. – uranix Nov 23 '16 at 18:36

You get the difference equation $$(h+2ϵ)y_{i+1}−4ϵy_i-(h-2ϵ)y_{i−1}=0.$$ Its characteristic polynomial has the solutions $$q_{1,2}=\frac{4ϵ\pm\sqrt{16ϵ^2+4(h+2ϵ)(h-2ϵ)}}{2(h+2ϵ)}=\frac{2ϵ\pm h}{h+2ϵ}=\frac{2\pm h/ϵ}{2+h/ϵ}$$ which only will give sensible non-oscillating solutions for $h< ϵ$.

With $y_0=0$ you get the solution formula $$y_k=c\left(1-\frac{(1-\frac{1}{2nϵ})^k}{(1+\frac{1}{2nϵ})^k} \right)$$ The constant is obtained from $$1=c\left(1-\frac{(1-\frac{1}{2nϵ})^n}{(1+\frac{1}{2nϵ})^n}\right)$$ With $n\gg\frac1ϵ$, or $h\llϵ$, this is close to and converges towards the exact solution $$y(t)=\frac{1-e^{-t}}{1-e^{-1}}.$$