Proper use of the Secant method

I'm having difficulty utilizing the secant method is this code I'm developing. The program takes in inclination angles of a curve that is approximated has straight lines. A subprogram reads these angles and calculates what the Mach angle is to be. The nature of this equation requires that the Mach angle (beta) be solved by numerical methods...I choose the Secant method to avoid using the derivative in the Newton method. Two problems I'm having.

1) The Secant method requires that I have two initial guesses. betaguess(1) and betaguess(2). I choose betaguess(1) = incl_ang(i), where incl_ang(i) is the angle I'm solving beta for. The reason I choose this is because for oblique shock waves, beta will always be greater than the inclination angle. As for betaguess(2), I choose it to be 0.001 + betaguess(1). For the iterations I start at i = 3, but if I do this I won't solve for betaguess(1) and betaguess(2).

2) The while loop keeps running. For the values that I have supplied in the code, it hangs up at i = 3. Even if I adjust the tolerance value.

Is there something that I am messing up on.

The red dots are points that are used to calculate the inclination angle by taking the tangent of the slope. The number of data points are specified by the user as well as the Mach number and the attack angle. For debugging purposes, the attack angle is 0 degrees.

    %MACH ANGLE - Calculates the angle (beta) of the shockwave given the   inclination angle.

%function [beta_ang] = Mach_angle(Mach, incl_ang, gamma)
incl_ang = [1.22255465965392,1.40187754402271,1.56039222519828,1.38106934082949,1.20174645646069,1.02242357209190,0.843100687723107,0.663777803354314,0.484454918985502,0.372860505867722,0.328994564000931,0.285128622134138,0.241262680267350,0.197396738400557,0.153530796533770,0.109664854666978,0.0657989128001851,0.0219329709333962,0.0219329709333962,0.0219329709333962,0.0657989128001851,0.109664854666978,0.153530796533770,0.197396738400557,0.241262680267350,0.285128622134138,0.328994564000931,0.372860505867723,0,0.484454918985502,0.663777803354314,0.843100687723107,1.02242357209190,1.20174645646069,1.38106934082949,1.56039222519828,1.40187754402271,1.22255465965392];

incl_ang = pi/2 - (incl_ang);       %inclination angle relative to blunt body axis of symmetry
incl_ang_deg = incl_ang*180/pi;
Mach = 5;
gamma = 1.4;
format long e

beta_ang = zeros(1,length(incl_ang));

tol = 0.001;

f = @(beta_guess)2*cot(beta_guess).*((Mach^2*sin(beta_guess).^2-1)./(Mach.^2.*(gamma+cos(2*beta_guess))+2))-tan(incl_ang(i));

for i = 3:length(incl_ang)
  error = 1;

    beta_guess(1) = abs(incl_ang(i-2));
    beta_guess(2) = beta_guess(1) + 0.001;

    while (error > tol)

        beta_guess(i) = beta_guess(i-1) - (f(beta_guess(i-1)))*((beta_guess(i-1) - beta_guess(i-2))/(f(beta_guess(i-1)) - f(beta_guess(i-2))));
i
        if f(beta_guess(i)) > tol
            beta_guess(i) = beta_guess(i) + tol;
        elseif f(beta_guess(i)) < tol
            beta_ang(i) = beta_guess(i);
        end

    end

 end

 beta_deg = beta_ang*180/pi;

 end

Answer 1

[Roger Stafford]

Your code is not carrying out the secant method properly, Harold, but rather than my engaging in a prolonged discourse on how to correct it, I prefer to show you a shortcut which avoids having to use the secant method altogether.

The equality you are trying to solve is this:

 2*cot(b)*(M^2*sin(b)^2-1)/(M.^2.*(g+cos(2*b))+2) = tan(a),

using the substitutions b = beta_ang, M = Mach, g = gamma, and a = incl_ang for ease of notation. It can be shown that this equality is equivalent to the following cubic equation in tan(b):

 tan(a)*(M^2*(g-1)+2)*tan(b)^3-2*(M^2-1)*tan(b)^2+tan(a)*(M^2*(g+1)+2)*tan(b)+2 = 0

To solve this second equation there is no need to use the secant method or any of the other high-powered methods in matlab which require an initial estimate. You can simply use matlab's 'roots' function to solve for t = tan(b) in the cubic equation and then find the arctangent of t using matlab's 'atan' function to find b.

You will face three obstacles with this method. For some of your incl_ang values there will be one real root and two complex roots, so you must select the real root.

For other values of incl_ang there will be three real roots and you must choose which of these you want to use - they will each satisfy the above equation. You would have had this problem even with the secant method.

The third obstacle is that 'atan' will only give you values between -pi/2 and pi/2. If your desired value for b lies outside that range, you will have to correct it by adding (or subtracting) pi.

Presumably you will want to make the above selections in such a way that you get a continuous curve as incl_ang varies throughout the vector you have defined.

Roger Stafford

[Harold]

Thank you for the response. I will attempt to implement this during the next few days.

[Harold]

I've given some time to this problem and I tried using your method but none of the values for the roots match up to what's the online calculator finds. Here is the calculator that I am using to check my work. It's the Oblique Shock Relations calculator with the Mach set to 5 and the turn angle at weak. http://www.dept.aoe.vt.edu/~devenpor/aoe3114/calc.html.

As a test I used 1.22255465965392 rads (19.9527777777776 deg). The online calculator gives 29.7468 deg. The roots function gives 4.356445770919350e+00, -4.800000000000000e+01, 2.250830314974997e+01, 2.000000000000000e+00

I modified my code and it works pretty decently until it gets to Mach and theta combinations that give a detached shock, as tested with the Oblique Shock Relations calculator. I have yet to code in a check for shock detachment.

[Roger Stafford]

Something is amiss with your calculations, Harold. The angle 1.22255465965392 radians is not equal to 19.9527777777776 degrees. It is 70.04722218 degrees. Also the method I proposed involved the solution by 'roots' of a cubic equation which will produce three roots, not four. Don't forget that you must find the arctan of the root that you choose. Finally, remember that you may have to add or subtract some multiple of pi (180 degrees) to get whatever angle you have in mind.

[Harold]

Yes, I see what I did wrong. Fixing the code right now.