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?