Differential equation concerning echolocation

Question

I am in the process of trying to solve a differential equation using MATLAB.

Given information is the following: $$ c(z) = 4800 - 20.2090 + \frac{17.3368z}{1000}+272.9057e^\frac{-0.7528z}{1000} $$ $$ \frac{d^2z}{dx^2} = -q_0\frac{c'(z)}{c(z)^3} $$ $$ q_0=(\frac{c(z_0)}{cos(B_0)})^2 $$ $$ x = 0, z = z_0, \frac{dz}{dx} = tan(B_0) $$ B0 is the angle between the line z = 2000 and the ray in the start point.

I have solved the equation using MATLAB's ode45 command, but now I'm supposed to use the secant method to find which initial values of $B_0$ also give me an end value of $z = 2500$ at $x = 151900$. I've been told that there's 4, since I'm supposed to use the $\Delta t$ between the arrivals to pinpoint the initial position, but no matter how I try I cannot figure out how to use the secant method in this particular instance.

Currently my program looks something like this:

angle = (-10:0.0035:14); 

z_0 = 2000;



error_tolerance = 1;    % Väljer hur långt ljudvågen
it = 1;

figure
for i = 1:length(angle)

B_0 = tand(angle(i)); 
q_0 = (c_sound(z_0)/cosd(angle(i)))^2;    
[X,Z] = ode45(@soundwave_function, [0 distance], [z_0 B_0]);
if abs(Z(end,1)-2500) <= error_tolerance

    saved_angles(it,1) = angle(i); 

    it = it+1;                           

end

end

angles_str = num2str(saved_angles);

global saved_angles

if numel(saved_angles)>1

for i = 1:numel(saved_angles)
    disp([num2str(saved_angles(i)), ' degrees'])
end

else
    disp([num2str(saved_angles), ' grader'])
end

This does find the angles where $ z(151900) \approx 2500 $ in the interval specified, -10 to 14 degrees (with 0 degrees counting as z = 2000).

I don't know how to approach the problem using the secant method, and since I'm supposed to use that solution for the second bit I don't really know how to continue.

Is there anyone out there that knows how to apply the secant method in this specific case? Would be much appreciated!

Answer 1

You need to use a shooting method. Let's abstract it a bit.

Let's define $F(x_{ic})$ to be a function that takes your initial conditions as its argument and outputs the value of your solution at $x = 151900$.

In other words, $F$ is a solution trajectory of your ODE, and we might think of $F$ as a single run of ode45, executed with initial conditions $x_{ic}$.

If we want to compute $F(x_{ic}) - 2500 = 0$, we can find initial values that satisfy the problem. Therefore, we simply want to write a root finding algorithm.

Newton's method is a good fit, but it is hard to compute $\frac{dF}{dx_{ic}}$. Instead, we might use the secant method, also known as Broyden's method.

It is not so hard: simply implement a function for Broyden's method, and wrap your ODE solver in a function that takes initial conditions as its arguments. For instance:

function z151900 = myODE(IC)
      [x,z] = ode45(myode,[0 151900],IC);
      z151900 = z(end);
end;

function xic = SecantMethod(F,IC1,IC2)
      % your secant method iteration using two initial guesses
end;

IC1 = ic1;
IC2 = ic2;
SecantMethod(z151900,ic1,ic2);

I can post a more concrete example later, or you can search my answer history for "Broyden's method".