How to proceed solving this problem?

Question

I'm working on the problem "soundwaves under the water" (page 16 in the document is in English) from a numerical analysis book.

I've got the following problem that is taken from the numerical analysis book by Kahaner-Moler-Nash (P8-15). I've made the plot of the data point and the approximation model function:

z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 0 0 -1];
opt = optimset('MaxFunEvals',1000);
p=lsqnonlin(fun,x0,[],[],opt);
fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t);
tt=linspace(0,12000,1000);
plot(z,data,'r-',tt,fitf(tt),'b-');

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.

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

$\displaystyle\frac{d^2z}{dx^2}= \frac{-q_0c'(z)}{c(z)^3}$

where $\displaystyle q_0 = \left(\frac{c(z_0)}{\cos b_0}\right)^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 RK4 already in theory but I dont' know how to do it in matlab and the non-linear least squares is done above. Can you help me?

Answer 1

This looks like quite an easy system of ODEs to implement in Matlab and a good way to learn how to use the ODE solvers. I'd start by reading through this explanatory article from The MathWorks and by looking at the documentation for ode45. The examples may be particularly helpful. Finally this blog post may be helpful, particularly if you have difficulty transforming your second order equation to a system of two first order ODES.

Within Matlab itself you'll find additional examples that illustrate the capabilities of the various solvers (though you won't need event detection or dynamic plotting like many of these use). Run odeexamples. The orbitode and ballode examples are nice; rigidode is quite simple. You can view the code for any of these: type edit orbitode, for example.