Webster equation is a popular generalization of 1D wave equation used for ducts of variable cross-section $S \equiv S(x)$. Assuming harmonicity in time, the spatial equation for propagation of variable $\Phi(x)$ is:

$$ \frac{\partial^2 \Phi}{\partial x^2} + \frac{\partial (\ln S(x))}{\partial x}\frac{\partial \Phi}{\partial x} + k^2\Phi^2 = 0 $$

What FD scheme would be appropriate assuming that we have boundary conditions at $\Phi(0)$ and $\Phi(L)$, i.e. no $\partial_x\Phi(0)$ etc.?

  • $\begingroup$ I'm not sure I understand your point regarding boundary conditions. Do you mean you have $\phi(0)$ and $\phi(L)$ prescribed, i.e. dirichlet boundary conditions? $\endgroup$
    – Steve
    Feb 28, 2016 at 16:37
  • $\begingroup$ That's it. Sorry for confusing assignment. $\endgroup$ Feb 28, 2016 at 16:49

1 Answer 1


This is not something I've worked on, so I do not know if there is a recommended approach to this, only how I would approach this problem.

I'm assuming you can evaluate $C(x) = (\ln(S(x)))_x$ wherever you want.

Since $\phi$ is prescribed at $0$ and $L$, I would discretise on $x_i = L\frac{i}{N+1}$ for $i$ from $1$ to $N$. Then the finite difference operators could be combined into a matrix $A$ and boundary condition correction term vector $\mathbf{b}$ with $$\begin{gathered} A = \frac{1}{\Delta x^2}\delta^2_x + \frac1{2\Delta x}\mathrm{diag}(C(x_i))\delta_{2x}\,,\\ \mathbf{b}_i = \begin{cases} \left(\frac1{\Delta x^2} - \frac1{2\Delta x}c(x_1)\right)\phi(0)\quad&\text{if}\quad i=1\,,\\ \left(\frac1{\Delta x^2} + \frac1{2\Delta x}c(x_N)\right)\phi(L) &\text{if}\quad i=N\,,\\ 0&\text{otherwise}\,, \end{cases} \end{gathered} $$ so that $$\left[\phi_{xx} + C(x)\phi_x\right]_{x_i} \approx (A\underline{\phi} + \mathbf{b})_i\,.$$

To find the solution, I would use Newton's method, looking for the root of $$F(\underline{\phi}) = A\underline{\phi} + \mathbf{b} + k^2 \mathrm{diag}(\phi_i^2)$$ which has tridiagonal Jacobian $$\mathrm{J}(\underline{\phi}) = A + 2k^2 \mathrm{diag}(\phi_i)\,.$$

(For completion) If there were a Neumann boundary condition, for example $$\phi_x(0) = \alpha\,,$$ I would include $\phi(0)$ in the solution vector $\underline{\phi}$, then I would construct $\mathrm{b}$ by considering a ghost point at $x = -\Delta x$ satisfying $$\frac{\phi(\Delta x) - \phi(-\Delta x)}{2\Delta x} = \alpha\,.$$ I believe mixed boundary conditions wouldn't be much more complicated than non-zero Neumann.


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.