Proper use of the Secant method
조회 수: 5 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
답변 (1개)
Roger Stafford
2012년 12월 24일
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
댓글 수: 4
Roger Stafford
2012년 12월 31일
편집: Roger Stafford
2012년 12월 31일
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.
참고 항목
카테고리
Help Center 및 File Exchange에서 Gas Dynamics에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!