Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I've got the following problem that is taken from the numerical analysis book by Kahaner-Moler-Nash (P8-15):

The speed of sound in ocean water depends on pressure, temperature and salinity, all of which vary with depth in fairly complicated ways. Let z denote depth in feet under the ocean surface (so that the positive z axis points down) and let c(z) denote the speed of sound at depth z. We shall ignore the changes in sound speed observed in horizontal directions. It is possible to measure c(z) at discrete values of z; typical results can be found in the table. We need c(z) and also c0(z) between data points. Fit the data in a least-squares sense with the non-linear model function

c(z) = 4800 + p1 + p2*z/1000 + p3*e^(-p4*z/1000)

To obtain startguesses, solve the linear approximation problem with p4 = 1: Make a plot over the data points and the received model curve c(z). Since the sound speed varies with depth, sound rays will travel in curved paths. A fixed underwater point emits rays in all directions. Given a particular point and initial direction we would like to follow the ray path. Thus letting x be the horizontal coordinate we know the initial values: x = 0, z = z0, dz=dx = tan 0, where 0 denotes the angle between the horizontal line z = z0 and the ray in the start point.

z c(z)
0 5050
500 4980
1000 4930
1500 4890
2000 4870
2500 4865
3000 4860
3500 4860
4000 4865
5000 4875
6000 4885
7000 4905
8000 4920
9000 4935
10000 4950
11000 4970
12000 4990

The ray path z(x) is described by the following second order differential equation

d^2z/dx= -q0*c'(z)/c(z)^3

where q0 = (c(z0) / cos b0)^2

Use the Runge-Kutta method (or ode45) to trace the ray beginning at z0 = 2000 feet and b0 = 7.8 degrees. Follow the ray for 25 nautical miles (1 nautical mile is 6076 feet). Plot the curve z(x). You should find that the depth at xf = 25 nautical miles is close to 2500 feet.

Now suppose that a sound source at a depth of 2000 feet transmits to a receiver 25 miles away at a depth of 2500 feet. The above calculation shows that one of the rays from the source to the receiver leaves the source at an angle close to 7.8 degrees. Because of the nonlinearity of the equation there may be other rays leaving at different angles that reach the same receiver. Run your program for b0 in the range from 10 up to 14 degrees, plot the ray paths and print a table of the values z(xf).

We are interested in finding values of 0 for which z(xf) = 2500. Use an efficient algorithm to determine the rays which pass through the receiver. Discuss the accuracy of your results.

So how should I proceed with solving this? I should study the methods, I think I know RK¤ already but not the non-linear least squares. Can you help me?

share|cite|improve this question
up vote 2 down vote accepted

A short (and fairly general) answer.

We have some experimental data points $(x_{i},y_{i})$ and a model $$ y = f(x;\textbf{p})$$ which depends on some parameters. As we change the values of these parameters $\textbf{p}$ the sum of the squares of the differences change, i.e. $$\text{Objective}(\mathbf{p}) = \sum_{i} \left[ y_{i} - f(x_{i};\mathbf{p})\right]^2.$$ We want to find a vector $\mathbf{p}_{0}$ such that $\text{Objective}(\mathbf{p})$ is minimised at $\textbf{p}_{0}$.

MATLAB can deal with these sorts of problems fairly easily. It has a built in solver lsqnonlin which solves nonlinear least square problems (as you may have guessed from the name). If you are actually asking about how that solver works, this answer needs to be a lot longer, but I believe that lsqnonlin implements the Levenberg-Marquadt algorithm.

share|cite|improve this answer
Default for MATLAB is a trust region Gauss-Newton method, though Levenberg-Marquardt can be used too if specifically called for. It's quite related though - one can view the LM method as a trust region method with a weird rule for updating the trust region size. – Nick Alger May 17 '15 at 9:32

Your Answer


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.