Checking for tolerance in Newton's Method

1 view (last 30 days)
Zack
Zack on 11 Oct 2011
I have written code for Newton's Method to solve a system of 10 non linear equations that have been given.
The code works and produces the right roots x1,x2,x3..x10.
The only issue I'm having is that currently, I perform 100 iterations. I want my Newton Method to stop when the difference between the real root and the calculated root is less than some defined tolerance (10e-6).
Can you give me some advice on where to add that to my for loop so I don't iterate when it is not needed.
clc
clear all
format long
%Initial guess
x(1)=1;
x(2)=1;
x(3)=1;
x(4)=1;
x(5)=1;
x(6)=1;
x(7)=1;
x(8)=1;
x(9)=1;
x(10)=1;
x = [x(1);x(2);x(3);x(4);x(5);x(6);x(7);x(8);x(9);x(10)];
dx = [1;1;1;1;1;1;1;1;1;1];
S = x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10);
count=0;
tol=0.000001;
R=4.056734;
n=100;
F=[ x(1)+x(4)-3;
2*x(1)+x(2)+x(4)+x(7)+x(8)+x(9)+2*x(10)-R-10;
2*x(2)+2*x(5)+x(6)+x(7)-8;
2*x(3)+x(5)-4*R;
x(1)*x(5)-0.193*x(2)*x(4);
x(6)*sqrt(x(2))-0.002597*sqrt(abs(x(2)*x(4)*S));
x(7)*sqrt(x(4))-0.003448*sqrt(abs(x(1)*x(4)*S));
x(4)*x(8)-0.00001799*x(2)*S;
x(4)*x(9)-0.0002155*x(1)*sqrt(abs(x(3)*S));
x(4)^2*(x(10)-0.00003846*S)];
a=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10) ));
b=(1/2)*x(6)/sqrt(x(2))-(0.1298500000e-2*(x(4)*(S)+x(2)*x(4)))/sqrt(x(2)*x(4)*(S));
c=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
d=-(0.1298500000e-2*(x(2)*(S)+x(2)*x(4)))/sqrt(x(2)*x(4)*(S));
e=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
f=sqrt(x(2))-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
g=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
h=-0.1298500000e-2*x(2)*x(4)/(sqrt(S));
j=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
k=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*S);
J1=[1 0 0 1 0 0 0 0 0 0];
J2=[2 1 0 1 0 0 1 1 1 2];
J3=[0 2 0 0 2 1 1 0 0 0];
J4=[0 0 2 0 1 0 0 0 0 0];
J5=[x(5) -0.193*x(4) 0 -.193*x(2) x(1) 0 0 0 0 0];
J6=[ a b c d e f g h j k];
J7=[-(0.1724000000e-2*(x(4)*S+x(1)*x(4)))/sqrt(abs(S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) (1/2)*x(7)/sqrt(abs(x(4)))-(0.1724000000e-2*(x(1)*S+x(1)*x(4)))/sqrt(abs(x(1)*x(4)*S))...
-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) sqrt(abs(x(4)))-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S))...
-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S))];
J8=[-0.1799e-4*x(2) -0.1799e-4*x(1)-0.3598e-4*x(2)-0.1799e-4*x(3)-0.1799e-4*x(4)-0.1799e-4*x(5)-0.1799e-4*x(6)-0.1799e-4*x(7)-0.1799e-4*x(8)-0.1799e-4*x(9)-0.1799e-4*x(10) -0.1799e-4*x(2)...
x(8)-0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2) x(4)-0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2)];
J9=[-0.2155e-3*sqrt(abs(x(3)*S))-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*S/sqrt(abs(x(3)*S))...
x(9)-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S))...
-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) x(4)-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S))];
J10=[-0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 2*x(4)*(.99996154*x(10)-0.3846e-4*x(1)-0.3846e-4*x(2)-0.3846e-4*x(3)-0.3846e-4*x(4)-0.3846e-4*x(5)-0.3846e-4*x(6)-0.3846e-4*x(7)-0.3846e-4*x(8)-0.3846e-4*x(9))-0.3846e-4*x(4)^2 ...
-0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 .99996154*x(4)^2];
J=[J1;J2;J3;J4;J5;J6;J7;J8;J9;J10];
for i=1:100
%while (max(dx)>tol)
solution = gaussElim(J,-F);
x=x+solution';
%dx=abs(x-solution')
F=[ x(1)+x(4)-3;
2*x(1)+x(2)+x(4)+x(7)+x(8)+x(9)+2*x(10)-R-10;
2*x(2)+2*x(5)+x(6)+x(7)-8;
2*x(3)+x(5)-4*R;
x(1)*x(5)-0.193*x(2)*x(4);
x(6)*sqrt(x(2))-0.002597*sqrt(abs(x(2)*x(4)*S));
x(7)*sqrt(x(4))-0.003448*sqrt(abs(x(1)*x(4)*S));
x(4)*x(8)-0.00001799*x(2)*S;
x(4)*x(9)-0.0002155*x(1)*sqrt(abs(x(3)*S));
x(4)^2*(x(10)-0.00003846*S)];
a=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10) ));
b=(1/2)*x(6)/sqrt(x(2))-(0.1298500000e-2*(x(4)*(S)+x(2)*x(4)))/sqrt(x(2)*x(4)*(S));
c=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
d=-(0.1298500000e-2*(x(2)*(S)+x(2)*x(4)))/sqrt(x(2)*x(4)*(S));
e=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
f=sqrt(x(2))-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
g=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
h=-0.1298500000e-2*x(2)*x(4)/(sqrt(S));
j=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*(S));
k=-0.1298500000e-2*x(2)*x(4)/sqrt(x(2)*x(4)*S);
J1=[1 0 0 1 0 0 0 0 0 0];
J2=[2 1 0 1 0 0 1 1 1 2];
J3=[0 2 0 0 2 1 1 0 0 0];
J4=[0 0 2 0 1 0 0 0 0 0];
J5=[x(5) -0.193*x(4) 0 -.193*x(2) x(1) 0 0 0 0 0];
J6=[ a b c d e f g h j k];
J7=[-(0.1724000000e-2*(x(4)*S+x(1)*x(4)))/sqrt(abs(S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) (1/2)*x(7)/sqrt(abs(x(4)))-(0.1724000000e-2*(x(1)*S+x(1)*x(4)))/sqrt(abs(x(1)*x(4)*S))...
-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) sqrt(abs(x(4)))-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S))...
-0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S)) -0.1724000000e-2*x(1)*x(4)/sqrt(abs(x(1)*x(4)*S))];
J8=[-0.1799e-4*x(2) -0.1799e-4*x(1)-0.3598e-4*x(2)-0.1799e-4*x(3)-0.1799e-4*x(4)-0.1799e-4*x(5)-0.1799e-4*x(6)-0.1799e-4*x(7)-0.1799e-4*x(8)-0.1799e-4*x(9)-0.1799e-4*x(10) -0.1799e-4*x(2)...
x(8)-0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2) x(4)-0.1799e-4*x(2) -0.1799e-4*x(2) -0.1799e-4*x(2)];
J9=[-0.2155e-3*sqrt(abs(x(3)*S))-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*S/sqrt(abs(x(3)*S))...
x(9)-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S))...
-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) x(4)-0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S)) -0.1077500000e-3*x(1)*x(3)/sqrt(abs(x(3)*S))];
J10=[-0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 2*x(4)*(.99996154*x(10)-0.3846e-4*x(1)-0.3846e-4*x(2)-0.3846e-4*x(3)-0.3846e-4*x(4)-0.3846e-4*x(5)-0.3846e-4*x(6)-0.3846e-4*x(7)-0.3846e-4*x(8)-0.3846e-4*x(9))-0.3846e-4*x(4)^2 ...
-0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 -0.3846e-4*x(4)^2 .99996154*x(4)^2];
J=[J1;J2;J3;J4;J5;J6;J7;J8;J9;J10];
%end
end
fprintf('x(1) %2.5f\n',x(1))
fprintf('x(2) %2.5f\n',x(2))
fprintf('x(3) %2.5f\n',x(3))
fprintf('x(4) %2.5f\n',x(4))
fprintf('x(5) %2.5f\n',x(5))
fprintf('x(6) %2.5f\n',x(6))
fprintf('x(7) %2.5f\n',x(7))
fprintf('x(8) %2.5f\n',x(8))
fprintf('x(9) %2.5f\n',x(9))
fprintf('x(10) %2.5f\n',x(10))
  1 Comment
Walter Roberson
Walter Roberson on 11 Oct 2011
I am unclear on this. If you know the real roots, then why bother going through the calculation? If you do not know the real roots, then how are you going to be able to determine if the difference between the calculated root and the real root is within a specific tolerance?

Sign in to comment.

Answers (0)

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!