1 Lab2B: Newton’s Method, Explorations and Error
In this lab you will be implementing Newton’s method for solving nonlinear equations of the form:
f(x) = 0 . (1)
As you work through the lab you should use some examples which you know will work with each of the methods.
Copyright By PowCoder代写 加微信 powcoder
Task 1: Implement Newton’s method
Newton’s method is a fixed point iteration (FPI) method, so here you use your work from Lab2A on FPI to help implement Newton’s method:
∆xk = − f(xk) (2) f′(xk)
xk+1 = xk + ∆xk (3) Note that in this Lab you will be storing values in vectors with components xk (or in Matlab
notation x(k)). The header below outlines what your code should take as input and output:
Note that you should implement this incrementally, not trying to add all the conditions at once, since this is much harder to debug. For example just start with the main loop and the maxiters stopping condition.
Testing your code
To ensure that your function mynewton is working correctly, you should test the function on some problems which you know the solution for. For example, the function
f(x) = (x + 1)(x − 1/2) 1
has roots at x = −1 and x = 1/2, and the derivative is f′(x) = 2x+1/2.
Thus, the Matlab code
[x1,fl]=mynewton(f,df,-1.2,0.001,10); [x2,fl]=mynewton(f,df,0.6,0.001,10);
should give output vectors x1 and x2 of the iterates that yield x1 ≈ −1 and x2 ≈ 0.5. You will also need to look at the order of convergence α, as discussed in lectures, we have
|xk+1 − x∗| = |ek+1| ≈ C (4) |xk − x∗|α |ek|α
multiplying by the denominator and taking logarithms gives,
loge |ek+1| ≈ α loge |ek| + loge (C) (5)
which is a linear model in loge |ek| with slope α and intercept logC. For your Newton method do the following, on your example:
Plot the logarithms of the errors against each other.
By hand draw what you think is the line which best fits the data.
Estimate the slope of the line (the order of convergence) and compare it to the expected value. Do the different schemes meet expectations for their relative order of convergence? Why or why not?
Fractals from Newton’s method (Task 2) View this lecture on YouTube
The three cube roots of unity are the solutions of z3 = 1. We can solve this equation using Eu- ler’s formula, exp (iθ) = cos θ + i sin θ. We write z3 = exp (i2πn), with n an integer, and take the cube root to find the solutions z = exp ( i 2πn/3). The three unique complex roots correspond to n = and 2, and using the trig values of special angles, we find
r1 =1, r2 =−12+ 23i, r3 =−12− 23i.
These roots can be located on the unit circle in the complex plane, and are shown below.
We now define the root-finding problem to be f (z) = 0, where f (z) = z3 − 1, with derivative f ′ (z) = 3z2. Newton’s method for determining the three complex roots of unity is given by the iteration
zn+1=zn− f(z). f′(z)
Here, we want to determine which initial values of z0 in the complex plane converge to which of the three cube roots of unity. So if for a given z0, the iteration converges to r1, say, we will color red the point z0 in the complex plane; if it converges to r2, we will color z0 green; and if to r3, blue.
In the next lecture, we show you how to construct a Matlab program that grids up a rectangle in the complex plane, and determines the convergence of each grid point z0. By coloring these grid points, we will compute a beautiful fractal.
From the Text Numerical Methods for Engineers, by . Task 2: Find the solutions of z^4 = 1 as complex numbers and illustrate them in terms of exp and the graph above.
Coding the Newton fractal (Tasks 3) View this lecture on YouTube
Our code first defines the function f (z) and its derivative to be used in Newton’s method, and defines
the three cube roots of unity:
f = @(z) z.^3-1; fp = @(z) 3*z.^2;
root1 = 1; root2 = -1/2 + 1i*sqrt(3)/2; root3 = -1/2 – 1i*sqrt(3)/2;
We let z = x + iy. The complex plane is represented as a two-dimensional grid and we define nx and ny to be the number of grid points in the x- and y-directions, respectively. We further define xmin and xmax to be the minimum and maximum values of x, and similarly for ymin and ymax. Appropriate values were determined by numerical experimentation.
nx=2000; ny=2000;
xmin=-2; xmax=2; ymin=-2; ymax=2;
The grid in x and y are defined using linspace. x=linspace(xmin,xmax,nx); y=linspace(ymin,ymax,ny);
We then use meshgrid to construct the two-dimensional grid. Suppose that the x-y grid is defined by x=[x1 x2 x3] and y=[y1 y2 y3]. Then [X, Y]=meshgrid(x,y) results in the matrices
x1 x2 x3 X = x1 x2 x3 ,
and we can grid the complex plane using
[X,Y]=meshgrid(x,y);
y1 y1 y1 Y = y2 y2 y2 ,
We now iterate Newton’s method nit times. These lines constitute the computational engine of the code.
for n=1:nit
Z = Z – f(Z) ./ fp(Z);
We next test to see which roots have converged. We use the logical variables Z1, Z2 and Z3 to mark the grid points that converge to one of the three roots. Grid points that have not converged are marked by Z4, and our convergence criteria is set by the variable eps. The function abs returns the modulus of a complex number.
eps=0.001;
Z1 = abs(Z-root1) < eps; Z2 = abs(Z-root2) < eps;
Z3 = abs(Z-root3) < eps; Z4 = ~(Z1+Z2+Z3);
Finally, we draw the fractal. The appropriate graphing function to use here is image, which can color pixels directly. We first open a figure and set our desired colormap. Here, our map will be a four- by-three matrix, where row one of the matrix corresponds with the RGB triplet [1 0 0] that specifies red; row two of the matrix [0 1 0] specifies green; row three of the matrix [0 0 1] specifies blue; and row four of the matrix [0 0 0] specifies black. The numbers 1-2-3-4 in our image file will then be colored red-green-blue-black.
map = [1 0 0; 0 1 0; 0 0 1; 0 0 0]; colormap(map);
We construct the image file by combining the values of our four Z matrices into a single Z matrix containing elements 1, 2, 3, or 4.
Z=(Z1+2*Z2+3*Z3+4*Z4);
The graph is created in our final lines of code. We use the limits of x and y to specify the location of our image in the complex plane. One needs to be aware that the function image assumes that the first row of pixels is located at the top of the image. So by default, image inverts the y-axis direction by setting the ’YDir’ property to ’reverse.’ We need to undo this when plotting data from a computation because the usual convention is for the first row of pixels to be at the bottom of the image. We therefore set the ‘YDir’ property to ‘normal.’
image([xmin xmax], [ymin ymax], Z); set(gca,'YDir','normal');
xlabel('$x$', 'Interpreter', 'latex', 'FontSize',14);
ylabel('$y$', 'Interpreter', 'latex', 'FontSize',14);
title('Fractal from $f(z)=z^3-1$', 'Interpreter', 'latex','FontSize', 16)
The resulting plot looks like
From the Text Numerical Methods for Engineers, by .
Task 3: Implement the Matlab above for z^4 = 1. Discuss aspects of the implementation. Watching the 2 excellent videos by Chasnov will also be helpful.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com