k-Wave Toolbox |

## Comparison Of Modelling Functions Example

On this page…
Overview About the modelling functions Running the simulation Numerical comparison

### Overview

This example provides a short comparison between the simulation functions `kspaceFirstOrder2D` and `kspaceSecondOrder`. It builds on the Homogeneous Propagation Medium and Using A Binary Sensor Mask examples.

- open the file in the MATLAB Editor
- run the file in MATLAB

For a more detailed discussion of the second order model used in k-Wave, see Treeby, B. E. and Cox, B. T., "A k-space Green's function solution for acoustic initial value problems in homogeneous media with power law absorption," J. Acoust. Soc. Am., vol. 129, no. 6, pp. 3652-3660, 2011.

Back to Top

### About the modelling functions

In the previous examples, the acoustic simulations are all performed using `kspaceFirstOrder2D`. This function is based on the sequential calculation of the acoustic pressure, particle velocity, and acoustic density using three coupled first-order equations (the conservation of mass, conservation of momentum, and an equation of state). For a homogeneous medium, these equations can also be combined into a single second-order acoustic wave equation and solved using Green's function methods. The function `kspaceSecondOrder` is an efficient numerical implementation of this solution specifically for initial value problems. For homogeneous media, the two approaches give the same results.

In comparison, the function `kspaceSecondOrder` has greater computational efficiency and the time steps can be arbitrarily large (as the solution is exact). It also allows both an initial pressure and an initial pressure gradient to be defined (see the Setting An Initial Pressure Gradient Example). Power law absorption is also more accurately encoded (see the Modelling Power Law Absorption Example). However, `kspaceSecondOrder` has slightly less functionality than its first-order counterparts (`kspaceFirstOrder1D`, `kspaceFirstOrder2D`, and `kspaceFirstOrder3D`). For example, it only supports homogeneous media and binary sensor masks, and does not support time varying sources or the full range of visualisation options.

The second-order code does not implement an absorbing boundary layer (see Controlling The Absorbing Boundary Layer Example for more details on the perfectly matched layer used in `kspaceFirstOrder2D`). Instead, to delay the appearance of wave wrapping, it extends the size of the computational grid. Although this reduces the computational efficiency, it allows the exactness of the solution to be maintained. The grid extension is controlled by setting the optional input parameter `'ExpandGrid'` to `true`. Note, the wrapped waves will still appear if the simulation time is longer than the time it takes for the waves to propagate through the grid extension.

Back to Top

### Running the simulation

The function `kspaceSecondOrder` (this is used for all dimensions) takes the same input structures as the first-order simulation functions. For example, to run the Homogeneous Propagation Medium Example, the Cartesian sensor mask is first converted to a binary sensor mask using `cart2grid`, and then `kspaceSecondOrder` is called with `'ExpandGrid'` set to `true`.

    % convert the cartesian sensor mask to a binary sensor mask
    sensor.mask = cart2grid(kgrid, sensor.mask);

    % run the simulation using the second order code
    sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'ExpandGrid', true);

Back to Top

### Numerical comparison

When the medium is not absorbing and no compensation for wave wrapping is used, the two simulation functions give identical results to machine precision. A plot of the signals recorded at the first sensor position and the difference between them is given below (set `example_number = 1` within the example m-file).

When the medium is absorbing and wave wrapping is mitigated using a perfectly matched layer, there is a small error in the numerical result computed using `kspaceFirstOrder2D`. A plot of the signals recorded at the first sensor position and the difference between them is given below (set `example_number = 4` within the example m-file)

Back to Top

© 2009-2012 Bradley Treeby and Ben Cox.