I'm solving an assignment in numerical analysis where I use this model function for soundwaves under the water after fitting the model function in a least-squares sens and finding the coefficients.
$c(z)=4800+20.2090+(17.3368)\frac{z}{1000}+(272.9057)exp(-\frac{-0.7528z}{1000})$ and
$c'(z)=(-0.7528*272.9067/1000)exp(-\frac{-0.7528z}{1000}) +17.3368/1000$
Is the above correct? I get the expected answer ("You should find that the depth at xf = 25 nautical miles is close to 2500 feet.") but the graph doesn't look as expected. It seems to indicate that the sound wave is moving down, then up and then down again. Could it be correct or is there mistakenly done?
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.
I got some help at scicomp https://scicomp.stackexchange.com/questions/8346/ode45-usage-in-this-case but I'm not sure about the solution.
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.
function dZ=sys(x,Z)
c=@(z)4800 + 20.2090 + (17.3368)*z/1000+ (272.9057)*exp(-z*0.7528/1000); % c(z)
c=c(2000);
% Z(1):=z
% Z(2):=u
dZ=zeros(2,1); % a column vector
dZ(1)=Z(2);
dZ(2)=-(c/cosd(7.8))^2*(((-272.9057*0.7528/1000)*exp(-Z(1)*0.7528/1000)) + 17.3368/1000)/...
(4800 + 20.2090 + (17.3368)*Z(1)/1000+ (272.9057)*exp(-Z(1)*0.7528/1000))^3;
end
I ran the above code with the script
x=0:0.5:6076*25;
[X,Z]=ode45(@sys,x,[2000 tand(7.8)]);
plot(X,Z(:,1),'r')
Then I get this graph which actually has the correct value (2500) at 25 nautical miles. Is the solution correct? Then I iterate over initial values to find which go to the receiver. Could that be feasible?
hold on
x=0:5:6076*25;
for i=-10:14
[X,Z]=ode45(@sys,x,[2000 tand(i)]);
plot(X,Z(:,1),'r') %Z(:,1) is z(x) and Z(:,2) is z'(x).
i=i+1;
end
hold off
This is not an answer but it is too long for a comment.
First, your double minus sign making a problem to me, I started from your data and fitted, just as you did, the model $$c(z)=4800+a+\frac{b z}{1000}+c e^{\frac{d z}{1000}}$$ and obtained $a=-20.209$ (not $+20.209$), $b=17.3368$, $c=272.906$, $d=-0.752778$ and the fit is effectively very good.
Concerning the derivative, it write $$c'(z)=\frac{1}{1000}\Big(b+c\, d \,e^{\frac{d z}{1000}}\Big)$$ and it only cancels once at $$z_{*}=\frac{1000 }{d}\log \left(-\frac{b}{c d}\right)$$ ($3284$ using the numbers); so $c(z)$ decreases for $0<z<z_{*}$ and increases for $z_{*}<z$; it never decreases again.