This thesis examines the nonlinear behaviour of thermoacoustic systems by using approaches from the field of nonlinear dynamics. The underlying behaviour of a nonlinear system is determined by two things: first, by the type and form of the attractors in phase space, and second, by the mechanism that the system transitions from one attractor to another. For a thermoacoustic system, both of these things must be understood in order to define a safe operating region in parameter space, where no high-amplitude oscillations exist. Triggering in thermoacoustics is examined in a simple model of a horizontal Rijke tube. A triggering mechanism is presented whereby the system transitions from a stable fixed point to a stable limit cycle, via an unstable limit cycle. The practical stability of the Rijke tube was investigated when the system is forced by stochastic noise. Low levels of noise result in triggering much before the linear stability limit. Stochastic stability maps are introduced to visualise the practical stability of a thermoacoustic system. The triggering mechanism and stochastic dependence of the Rijke tube match extremely well with results from an experimental combustor. The most common attractors in thermoacoustic systems are fixed points and limit cycles. In order to define the nonlinear behaviour of a thermoacoustic system, it is therefore important to find the regions of parameter space where limit cycles exist. Two methods of finding limit cycles in large thermoacoustic sytems are presented: matrix-free continuation methods and gradient methods. Continuation methods find limit cycles numerically in the time domain, with no additional assumptions other than those used to form the governing equations. Once the limit cycles are found, these continuation methods track them as the operating condition of the system changes. Most continuation methods are impractical for finding limit cycles in large thermoacoustic systems because the methods require too much computational time and memory. In the literature, there are therefore only a few applications of continuation methods to thermoacoustics, all with low-order models. Matrix-free shooting methods efficiently calculate the limit cycles of dissipative systems and have been demonstrated recently in fluid dynamics, but are as yet unused in thermoacoustics. These matrix-free methods are shown to converge quickly to limit cycles by implicitly using a 'reduced order model' property. This is because the methods preferentially use the influential bulk motions of the system, whilst ignoring the features that are quickly dissipated in time. The matrix-free methods are demonstrated on a model of a ducted 2D diffusion flame, and the safe operating region is calculated as a function of the Peclet number and the heat release parameter. Both subcritical and supercritical Hopf bifurcations are found. Physical information about the flame-acoustic interaction is found from the limit cycles and Floquet modes. Invariant subspace preconditioning, higher order prediction techniques, and multiple shooting techniques are all shown to reduce the time required to generate bifurcation surfaces. Two types of shooting are compared, and two types of matrix-free evaluation are compared. The matrix-free methods are also demonstrated on a model of a ducted axisymmetric premixed flame, using a kinematic G-equation solver. The methods find limit cycles, period-2 limit cycles, fold bifurcations, period-doubling bifurcations and Neimark-Sacker bifurcations as a function of two parameters: the location of the flame in the duct, and the aspect ratio of the steady flame. The model is seen to display rich nonlinear behaviour and regions of multistability are found. Gradient methods can also efficiently calculate the limit cycles of large systems. A scalar cost function is defined that describes the proximity of a state to a limit cycle. The gradient of the cost function is used in an optimisation routine to iteratively converge to a limit cycle (or fixed point). The gradient of the cost function is found with a forwards-backwards process: first, the direct equations are marched forwards in time, second, the adjoint equations are marched backwards in time. The adjoint equations are derived by partially differentiating the direct governing equations. The gradient method is demonstrated on a model of a horizontal Rijke tube. This thesis describes novel nonlinear analysis techniques that can be applied to coupled systems with both advanced acoustic models and advanced flame models. The techniques can characterise the rich nonlinear behaviour of thermoacoustic models with a level of detail that was not previously possible.