Newton-Raphson Method- Same Equations, 2 correct roots, 1 incorrect root

128 Views Asked by At

So I'm trying to calculate the roots for $\phi_1$ to verify this equation to determine the angle of a blade for an Archimedes Screw turbine based upon 'THE TURN OF THE SCREW: OPTIMAL DESIGN OF AN ARCHIMEDES SCREW' (Rorres, 2000):

$$\sin\biggl(\phi_1-\frac{2\pi}{N}\biggr)+\lambda(\phi_1-\phi_0)-\rho \sin(\phi_0)=0$$

$\lambda,\rho,\phi_0$ and N are all known with values of:

  • $\lambda$ = 3/16
  • $\rho$ = 1/2,
  • $\phi_0=\cos^{-1}(-\frac{\lambda}{\rho})$
  • N = 8

With these inputs the value of $\phi_1$ in the paper is listed as 441.96°, with boundary angles $\phi_L$ and $\phi_R$ equal to 162.68° and 343° respectively.

I've written the following MatLab script to determine the roots numerically using the Newton-Raphson method:

%% Newton-Raphson Method

clear all
close all
clc

% Vetruvius' Screw
N = 8;
a = 3/16;
b = 1/2;

% Phi 0: Bucket Begins
phi_0 = acos(-(a/b));

% Phi 1: Bucket Ends
f1 = @(x1) sin(x1-(2*pi/N)) + a*x1 - a*phi_0 - b*sin(phi_0);
df1 = @(x1) cos(x1-(2*pi/N)) + a;

% Phi L: Minimum Angular Boundary L
fl = @(xl) sin(xl) + a*xl - a*phi_0 - b*sin(phi_0);
dfl = @(xl) cos(xl) + a;

% Phi R: Minimum Angular Boundary L
fr = @(xr) sin(xr) + a*xr - a*phi_0 - b*sin(phi_0);
dfr = @(xr) cos(xr) + a;

% Plot functions to estimate initial roots
figure(1)
fplot(f1,[-4*pi 4*pi])
title('Initial Functions')
xlabel('x')
ylabel('y')
xticks([-4*pi -3*pi -2*pi -pi 0 pi 2*pi 3*pi 4*pi])
grid on
grid minor
hold on
fplot(fl)
hold on
legend(' Phi 1 ',' Phi L/R ')

% initial guesses
x1(1) = 7.714;  
xl(1) = 2.8;
xr(1) = 6;

for i = 1:1000
    x1(i+1) = x1(i) - (f1(x1(i))/df1(x1(i)));
    xl(i+1) = xl(i) - (fl(xl(i))/dfl(xl(i)));
    xr(i+1) = xr(i) - (fr(xr(i))/dfr(xr(i)));
end

phi_1 = x1(1000) * (180/pi);
phi_l = xl(1000) * (180/pi);
phi_r = xr(1000) * (180/pi);

figure(2)
plot(x1)
title('Newton-Raphson Method to determine the Blade Angles of an 
Archimedes Screw Turbine')
xlabel('Number of Iterations')
ylabel('Blade Angle (rads)')
hold on
plot(xl)
plot(xr)

I'm expecting a value of about 7.71 rads (441.96°) for $\phi_1$ but my answer is much lower, using the same method the other two angles ($\phi_L$ and $\phi_R$) are pretty much exactly the values I'm expecting at 162.2° and 343.1°.

Can anyone see a reason I'm getting correct values for $\phi_L$ and $\phi_R$ yet incorrect values for $\phi_1$ using the same method?

(Attached are plots of both the functions and the resulting roots from my script.)

EDIT Also attached the results from the paper for comparison.

Functions

Resulting Roots

Resulting Roots with modified loop

Expected Angles

1

There are 1 best solutions below

2
On BEST ANSWER

You probably misread the value for N. Isolating N from the formula for the given value suggests N=3, and indeed the iteration with N=3 has the values

  k           x1                         xl                     xr

  1 : 8.000000000000000          3.000000000000000       6.000000000000000
  2 : 7.730348724105118          2.842377212981511       5.986517565152984
  3 : 7.713755074290903          2.839259277097096       5.986495020153686
  4 : 7.713669321922058          2.839257403087327       5.986495020088732
  5 : 7.713669319598735          2.839257403086645       5.986495020088731
  6 : 7.713669319598736          2.839257403086645       5.986495020088732
  7 : 7.713669319598735          2.839257403086645       5.986495020088731
  8 : 7.713669319598736          2.839257403086645       5.986495020088732