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?