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;

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;                           



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'])

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

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!


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);

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

IC1 = ic1;
IC2 = ic2;

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.