Solution of the Helmholtz equation

Directly to the input form

• The Helmholtz equation is the classical example of an elliptic partial differential equation in two space dimensions. It reads
  Δ u + λ × u = f, Δ=(d/dx)²+(d/dy)²
  λ=0 includes the Poisson equation. We consider here only an axis parallel rectangle as region for (x,y):
  (x,y) in [a,b]x[c,d] and λ <= 0. Under these restrictions the equation has an unique solution, given properly defined boundary conditions. For continuous f and continuous boundary conditions u is continuous on the closed rectangle and two times continuously differentiable on the open one. Observe that the second derivatives need not be continuous on the closed rectangle.
• Elliptic partial differential equations describe a physical system in steady state. You might imagine of u as the deformation of a thin membrane under external load f and under the influence of its own mass.
• On the boundary of [a,b]x[c,d] where must be given boundary conditions in order to get a well posed problem. These might be prescriptions of the function value (so called Dirichlet data), prescriptions of the directional derivative in the direction normal to the boundary (von Neumann conditions) or a proper linear combination of these (Robin conditions). Robin conditions are not implemented here.
• The problem is solved here using the standard finite difference method with the standard discretization of the second order partial derivatives resulting in the so called 5 point star:
  (1/δ_x**2)*(u_i-1,j-2*u_i,j+u_i+1,j)+ (1/δ_y**2)*(u_i,j-1-2*u_i,j+u_i,j+1)+λ*u_i,j =f_i,j
  
• The rectangle [a,b]x[c,d] becomes discretized by a finite grid (x_i,y_j), i=0,....,n, j=0,...,m with equidistant values x_i and y_j, the values 0 and n resp. m corresponding to the boundaries. u_i,j in the equation above is an approximate value for the true solution at (x_i,y_j). The difference equation is solved for 1 <=i <= n-1, 1 <= j < m-1.
• The discretization error of this method is second order in δ_x, δ_y (measured in the max norm) provided the solution u is four times differentiable on the closed rectangle. Weaker smoothness allows weaker error estimates only. The discrete system is a linear system of equations which is solved here by a special routine.
• the computational routine used is hwscrt from netlib/fishpack.

Input

• There are two predefined cases with known exact solution (on the unit square) and hence known discretization error and you also have the possibility to define a problem of your own.
• In the first predefined case the solution is C⁴ in the closed rectangle and hence the convergence is of order O(h²) with h=max(δ_x,δ_y). Here λ =0 and we have Dirichlet data on all of the boundary. (We first choose a function as a solution and derive the data from this)
• In the second case u is C³ in the interior of the rectangle only, and only C² on the boundary. Hence convergence is remarkably weaker here. Again we have λ =0 but von Neumann data for x=0 and y=0 and Dirichlet data for x=b and y=d.
• If you want to specify a problem of your own, then you must give
  1. The function f(x,y). You define it as a formula in the language FORTRAN.
  2. The specification of the type of boundary value for each of the four boundary parts of the rectangle: either Dirichlet or von Neumann. This can be done individually for each of the four edges.
  3. The formulas for the functions defining the boundary values. attention: you specify the partial derivative on the boundary for von Neumann data, hence (d/dx)u at x=a, not the outer normal! and similar for y.
  4. The parameter λ. Important : we allow only λ <= 0 !.
  5. The parameters a,b,c,d defining the rectangle.
• The parameters n and m defining the gridsizes.
• Since you will get a 3D plot of the computed u you must specify your view point via the rotation angle for the x-axis and the u-axis.

Output

• You get a text output, either with the true solution and the maximal error or simply ''Helmholtz: u= user defined problem''.
• A 3D-Plot showing the piecewise linear interpolation of the discretized solution u over the rectangle, with some isolines.

Questions ?!

• Can you verify the theoretical convergence properties?
• What takes place in case of reduced smoothness properties? For example, play with discontinuous boundary data.
• In order to get four times differentiability on the closed rectangle, which conditions must be satisfied by the boundary functions?

To the input form