I have problem with understanding why I get the complex convergence when using Newton-Raphson's method (own implemented) to the convergence tables down below:
%These equations represent the different bows that belong to the
%two ellipses
f1 = @(x) sqrt(36 - 36/16.*(x-4).^2) + 2;
fp1 = @(x) -3.*(x-4)./(2.*sqrt(-(x-8).*x)); %Derivative to f1
f2 = @(x) -sqrt(36 - 36/16.*(x-4).^2) + 2;
fp2 = @(x) 3.*(x-4)./(2.*sqrt(-(x-8).*x)); %Derivative to f2
f3 = @(x) 0.5.*x + sqrt(10 - 0.15.*x.^2);
fp3 = @(x) 0.5 + - 0.15.*x/sqrt(10-0.15.*x^2); %Derivative to f3
f4 = @(x) 0.5.*x - sqrt(10 - 0.15.*x.^2);
fp4 = @(x) 0.5 + + 0.15.*x/sqrt(10-0.15.*x^2); %Derivative to f4
n = 10; %Number of iterations
%Table over convergence to the first intersection point between the vertical
%ellipse's upper bow and the crooked ellipse's upper bow
xStart = 4;
f = @(x) f1(x)-f3(x);
fp = @(x) fp1(x)-fp3(x);
x = xStart;
disp('iteration x dx')
for iter = 1:n
dx = -f(x)/fp(x);
x = x + dx;
disp([iter x dx])
end
%Table over convergence to the second intersection point between the
%vertical ellipse's upper bow and the crooked ellipse's upper bow
xStart = 4;
f = @(x) f1(x) - f3(x);
fp = @(x) fp1(x) - fp3(x);
x = xStart;
disp('iteration x dx')
for iter = 1:n
dx = -f(x)/fp(x);
x = x + dx;
disp([iter x dx])
end
%Table over convergence to the intersection point between the vertical
%ellipse's upper bow and the crooked ellipse's lower bow
xStart = 8;
f = @(x) f1(x) - f4(x);
fp = @(x) fp1(x) - fp4(x);
x = xStart;
disp('iteration x dx')
for iter = 1:n
dx = -f(x)/fp(x);
x = x + dx;
disp([iter x dx])
end
%Table over convergence to the intersection point between the vertical
%ellipse's lower bow and the crooked ellipse's lower bow
xStart = 4;
f = @(x) f2(x) - f4(x);
fp = @(x) fp2(x) - fp4(x);
x = xStart;
disp('iteration x dx')
for iter = 1:n
dx = -f(x)/fp(x);
x = x + dx;
disp([iter x dx])
end
The output I get when running these codes is as following:
iteration x dx
1.0000 15.4861 11.4861
2.0000 + 0.0000i 5.5336 - 0.6701i -9.9525 - 0.6701i
3.0000 + 0.0000i 8.4377 + 0.9203i 2.9041 + 1.5904i
4.0000 + 0.0000i 6.8854 + 0.2759i -1.5523 - 0.6444i
5.0000 + 0.0000i 7.5608 - 0.1090i 0.6755 - 0.3849i
6.0000 + 0.0000i 7.4558 - 0.0122i -0.1050 + 0.0967i
7.0000 + 0.0000i 7.4578 + 0.0000i 0.0019 + 0.0122i
8.0000 + 0.0000i 7.4578 - 0.0000i 0.0001 - 0.0000i
9.0000 + 0.0000i 7.4578 - 0.0000i -0.0000 + 0.0000i
10.0000 + 0.0000i 7.4578 + 0.0000i -0.0000 + 0.0000i
iteration x dx
1.0000 15.4861 11.4861
2.0000 + 0.0000i 5.5336 - 0.6701i -9.9525 - 0.6701i
3.0000 + 0.0000i 8.4377 + 0.9203i 2.9041 + 1.5904i
4.0000 + 0.0000i 6.8854 + 0.2759i -1.5523 - 0.6444i
5.0000 + 0.0000i 7.5608 - 0.1090i 0.6755 - 0.3849i
6.0000 + 0.0000i 7.4558 - 0.0122i -0.1050 + 0.0967i
7.0000 + 0.0000i 7.4578 + 0.0000i 0.0019 + 0.0122i
8.0000 + 0.0000i 7.4578 - 0.0000i 0.0001 - 0.0000i
9.0000 + 0.0000i 7.4578 - 0.0000i -0.0000 + 0.0000i
10.0000 + 0.0000i 7.4578 + 0.0000i -0.0000 + 0.0000i
iteration x dx
1 8 0
2 8 0
3 8 0
4 8 0
5 8 0
6 8 0
7 8 0
8 8 0
9 8 0
10 8 0
iteration x dx
1.0000 -0.5192 -4.5192
2.0000 + 0.0000i 0.6812 + 1.5032i 1.2005 + 1.5032i
3.0000 + 0.0000i 1.7776 + 0.2145i 1.0963 - 1.2887i
4.0000 + 0.0000i 1.2888 - 0.0404i -0.4888 - 0.2549i
5.0000 + 0.0000i 1.3235 - 0.0006i 0.0347 + 0.0398i
6.0000 + 0.0000i 1.3234 + 0.0000i -0.0001 + 0.0006i
7.0000 + 0.0000i 1.3234 - 0.0000i -0.0000 - 0.0000i
8.0000 + 0.0000i 1.3234 - 0.0000i 0.0000 + 0.0000i
9.0000 + 0.0000i 1.3234 - 0.0000i 0.0000 + 0.0000i
10.0000 + 0.0000i 1.3234 - 0.0000i 0.0000 + 0.0000i
To me the result is not making any sense, since I have drawn the equations to respective ellipse in the graph and here is the picture that I got:
So again, my question is why do I get complex convergence when the solutions to the intersection points between the ellipses when (according to the graph) there should be only real solutions? Have I implemented something wrong or could it be something theoretically that I am not aware of?
Thanks for all responses in advance!
You have found two phenomena that Newton's method succumbs to.
Firstly observe that the function you are using it on is complex-valued outside of a certain domain, but the tangent to the curve at many points in that domain will cut the x-axis outside. Hence you get complex values along the way.
Secondly observe that in one case you start from a point where the derivative is infinite, and hence the tangent cuts the x-axis at the same x-coordinate and the iteration is stuck there.
To avoid these you should know exactly why and when Newton's method works. $\def\wi{\subseteq}$
Define $[m] = \{ z : |z| \le m \}$ for any complex $m$.
Given complex function $f$ such that $f(r) = 0$ and $f'(r) \ne 0$ and $f''$ exists and $\dfrac{f''}{2f'(r)} \in [m]$:
Let $a = f'(r)$.
Given $d \in [\frac{1}{5m}]$:
$f(r+d) \in ad+[\frac{1}{2}f'']d^2 \wi ad+a[m]d^2$ [by Taylor's theorem].
$f'(r+d) \in ad+[f'']d^2 \wi ad+2a[m]d^2$ [by Taylor's theorem].
$d - \dfrac{f(r+d)}{f'(r+d)} \in d - \dfrac{ad+a[m]d^2}{a+2a[m]d} = d - \dfrac{d+[m]d^2}{1+2[m]d}$
$ \wi d - (d+[m]d^2)(1+\frac{10}{3}[m]d)$ [because $\dfrac{1}{1+t} \in 1+[\frac{5}{3}]t$ for any $t \in [\frac{2}{5}]$]
$ \wi \frac{13}{3}[m]d^2+\frac{10}{3}[m]^2d^3$
$ \wi 5[m]d^2 \wi [1]d$.
Therefore Newton-Raphson converges quadratically if starting within $[\frac{1}{5m}]$ of $r$.
I leave you to check that both phenomena are due to either starting too far away or having unbounded $\dfrac{f''}{f'}$.