Newton’s method is one of the most powerful and well-known numerical methods for

Question

Newton’s method is one of the most powerful and well-known numerical methods for solving a root-finding problem f(x) = 0. The algorithm of Newton’s method is as following:

float Newton(float p0, float Tolerance) //initial approximation p0

{

   float p;

   while(1)

   {

            p=p0-f(p0)/f’(p0);

            if (fabs(p-p0)<Tolerance)

                        return p;

          p0=p;                               //update p0

   }

}             

In this project, let f (x)=Z4-1 and Tolerance = 0.001. Let z=x+iy in which x is the real and y the imaginary parts be a point in the range W=[-2 £ x £ 2, -2 £ y £ 2]. Now consider your screen in the VGA mode 640 x 480 to map to W.   In this case, z = x+iy becomes a pixel in the range 640 x 480. The Z4-1=0 has four root: -1, 1, i, -i. We can color each pixel in following way:

If z=x+iy converge the root –1 by Newton’s method, we color the corresponding pixel by Red.

If z=x+iy converge the root 1 by Newton’s method, we color the corresponding pixel by Green.

If z=x+iy converge the root i by Newton’s method, we color the corresponding pixel by Blue.

If z=x+iy converge the root -i by Newton’s method, we color the corresponding pixel by Black.

Of course, you can choose another four colors. Therefore, this project is that

you make a picture which is in 640x480 VGA mode such that each pixel corresponding to a point in W is colored by converging which root by Newton’s method.   

As the previous project, use clock() function to measure your running time. Try to speed up your program as fast as possible without concerning memory space. Output should include the running time of your program.

Explanation / Answer

float Newton(float p0, float Tolerance) //initial approximation p0

{

%tic for displaying running time of program

   float p;
% purpose compute destination roots for Newton's method at points in complex plane

x = linspace(xrange(1), xrange(2), N);
y = linspace(yrange(1), yrange(2), N);

[X,Y] = meshgrid(x,y);

z = X+i*Y;

% set roots
z1 = -1; z2 = 1; z3 = i; z4 = -i;

f (x)=Z4-1

Tolerance =0.001;

% initial value

p0=0.7

while(1)

   {

            p=p0-f(p0)/f’(p0);

            if (fabs(p-p0)<Tolerance)

                        return p;

          p0=p;                               //update p0

   }

}             

% color channels for plot (as an RGB image)
r = zeros(N, N); g = zeros(N, N); b = zeros(N, N); black = zeros(N, N);

% find roots (to a given tolerance tol)


% color starting values that got close to each root differently
ids = find(abs(z-z1)< Tolerance); r(ids) = 255; g(ids) = 255;
ids = find(abs(z-z2)< Tolerance); r(ids) = 255;
ids = find(abs(z-z3)< Tolerance); g(ids) = 255;
ids = find(abs(z-z4)< Tolerance); b(ids) = 255;

% build rgb image and draw it
c(:,:,1) = r; c(:,:,2) = g; c(:,:,3) = b; c(:,:,4) = black;
image(x, y, uint8(c)); axis equal; axis off; drawnow; pause(.05);
[time] = CLOCK();

return

function fz = f(z, z1, z2, z3, z4)

fz = (z-z1).*(z-z2).*(z-z3).*(z-z4);

return

function fprime = dfdz(z, z1, z2, z3, z4)

fprime = ...
      (z-z2).*(z-z3).*(z-z4) + ...
      (z-z1).*(z-z3).*(z-z4) + ...
      (z-z1).*(z-z2).*(z-z4) + ...
      (z-z1).*(z-z2).*(z-z3) ;

return


function[time] = CLOCK()

tic

time=toc;

Related Questions

Navigate

• Browse (All)
• Browse N
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.