Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

This is a somewhat basic question, I guess. Take the ODE boundary value problem $$ \frac1\lambda y''-y'=0, \qquad y(0) = 0, \quad y(1) = 1, $$ with the solution $$ y(x) = \frac{e^{x\lambda}-1}{e^\lambda -1}. $$ Discretize it on a uniform grid $x_k = kh$, $h=1/n$, $k=0,\ldots,n$, as follows: $$ \frac1\lambda \frac{y_{k-1}-2y_k+y_{k+1}}{h^2}-\frac{y_{k+1}-y_{k-1}}{2h} = 0, $$ and solve the resulting tridiagonal system, for, e.g., $n=100, \lambda=10^4$:

enter image description here

With increasing $n$, the scheme nevertheless converges.

a. Can you suggest a good reference that describes this kind of thing? I think that stiffness of initial-value problems is not quite the same thing, which is why I'm asking this question, but maybe I'm wrong.

b. Short of choosing a better method, is there a "filter" I can apply to the output $\mathbf{y}_{0:n}$ to recover a better approximation of the true solution (in particular, make the transformed output positive, as it should be).

share|improve this question

1 Answer 1

up vote 5 down vote accepted

This is a convection dominated problem. Since $1/\lambda$ is small, your equation is approximately $y'=0$. But the solution of this cannot satisfy your boundary conditions. The solution will rapidly change from 0 to 1 as you approach $x=1$, and a boundary layer is created at $x=1$. The term $y''/\lambda$ cannot be neglected in the boundary layer.

If you write

$$ y_k = \frac{1}{2}(1 + \lambda h/2) y_{k-1} + \frac{1}{2}(1 - \lambda h/2) y_{k+1} $$

the solution will be monotone if

$$ \lambda h < 2 $$

which is called the cell Peclet condition. So your mesh must be finer than

$$ h < \frac{2}{\lambda} $$

You can see that increasing $\lambda$ requires very fine mesh. If your mesh satisfies above condition you should get a good solution.

The other alternative is to use an upwind scheme. In your case it is (please check stability as above)

$$ \frac{1}{\lambda} \frac{y_{k-1} - 2 y_k + y_{k+1}}{h^2} - \frac{y_k - y_{k-1}}{h} = 0 $$

but this is only first order accurate.

Another way to analyze this is in terms of an M-matrix.

For more, you can refer to some texts, e.g.,

  1. Strikwerda: Finite Difference Schemes and Partial Differential Equations
  2. Wesseling: Principles of Computational fluid dynamics
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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