Description
1 Introduction
Last time we discussed Newton’s method for nonlinear equations in one real or complex variable. In this
lab, we will extend the discussion to two or more dimensions. We will consider a common application of
Newton’s method to nonlinear least squares fitting of a function to given data.
2 A vector function
In many cases, the nonlinear problem we are facing involves vectors. For example, X and F may represent
vectors of length 2, and we are asked to solve the system:
F0(X) = 1 − x0
F1(X) = 10(x1 − x
2
0
)
When we are working on a vector problem like this, the derivative information that we need is known as the
Jacobian, which is often symbolized by J. For this problem, we will have
J =
∂F0
∂x0
∂F0
∂x1
∂F1
∂x0
∂F1
∂x1
!
=
−1 0
−20×0 10
It would be nice to think that our old Newton code would work with this problem. That’s almost true, but
we need to make several modifications in order to adjust to the fact that vectors are involved.
After we get an initial vector version of Newton’s method working, we will explore some of the new issues
that arise when working with vector problems.
3 A Newton code for vector problems
Suppose x is a vector in R
n, x = (x1, x2, . . . , xn), and F(x) a differentiable vector function of x:
F = (Fi(x1, x2, . . . , xN )), for i = 1, 2, . . . , n,
with Jacobian matrix J,
Jij =
∂Fi
∂xj
.
Recall that Newton’s method can be written
x
(k+1) − x
(k) = −J(x
(k)
)
−1F(x
(k)
). (1)
The superscript indicates iteration number, rather than vector index.
As an implementation note, the inverse of J should not be computed. Instead, the system
J(x
(k)
)(x
(k+1) − x
(k)
) = −F(x
(k)
).
1
should be solved. As you will see later in this course, this will save about a factor of three in computing
time over computing the inverse. numpy.linalg provides the x=np.linalg.solve(A,b) function for this
solution. If you wish to solve the system
J dx = −F
then the solution can be written as
import numpy a s np
dx = np . l i n a l g . s o l v e ( J , −F )
4 Exercise 1
1. Make a copy of the file newton3.py from Lab 4, changing its name to vnewton1.py
2. Change the signature to
def vnewton1 ( func , d func , x0 , maxIts = 100 ) :
”””x , numIts = vnewton1 ( func , d func , x0 , maxI ts = 100 )
func : e v a l u a t e s t h e f u n c t i o n .
d func : e v a l u a t e s t h e Jacob ian .
x0 : t h e i n i t i a l s o l u t i o n e s t im a t e .
maxI ts : maximum number o f i t e r a t i o n s , d e f a u l t i n g t o 1 0 0 .
”””
. . . Old t e x t o f newton3 ( ) al g o ri t hm . . .
return x , numIts
and modify it so that it:
• sets the function vector F using func();
• sets the Jacobian matrix J using dfunc();
• uses the numpy.linalg function solve() instead of scalar division;
• uses the numpy.linalg function norm() instead of abs();
• has no print statements;
3. Create a new problem file v1.py, in which v1F() defines the function vector F and v1J() defines the
jacobian matrix J, for the problem whose function is:
f0(x) = 1 − x0
f1(x) = 10(x1 − x
2
0
)
Both F and J are numpy arrays, so you need to declare them first as arrays of zeros before using them.
For instance, inside the v1.py file, the definition of v1F() might look like this:
def v1F ( x ) :
import numpy a s np
F = np . z e r o s ( 2 )
F [ 0 ] = 1. 0 − x [ 0 ]
F [ 1 ] = 1 0. 0 ∗ ( x [ 1 ] − x [ 0 ] ∗ ∗ 2 )
return F
4. Test vnewton1() on problem v1, using the starting point x0=np.array([-1.2,2.0]). In order to
access the function and jacobian information from v1, you will need a statement like
from v1 import v1F , v1J
5. Check your result by printing the norm of the function, which should be very small.
2
5 A complex function revisited
It is possible to rewrite a complex function of a complex variable as twice as many real vector functions of
twice as many real vector variables. This is done by equating the real part of a complex number with the
first component of a vector and the imaginary part of a complex number with the second component of a
vector.
Previously, we used Newton’s method to find complex roots of f(z) = z
2 + 9. Let us replace the complex
argument z by the real vector x = [x0, x1] where z = x0 + x1j, and similarly f(z) = f0 + f1j. Here Python
understands an expression like 2j to be the imaginary quantity we would mathematically represent as 2i.
Plugging these variables in yields
f0 + f1j = (x
2
0 − x
2
1 + 9) + (2×0 x1)j.
This can be written in an equivalent real vector form as
f0(x0, x1)
f1(x0, x1)
=
x
2
0 − x
2
1 + 9
2 x0 x1
(2)
6 Exercise 2
1. Create a file v2.py containing v2F and v2J, defining the values of the function F and Jacobian J of this
problem as real-valued vector and matrix quantities, respectively. The definitions might look like this;
def v2F ( x ) :
. . . set up and c a l c u l a t e F . . .
return F
def v2J ( x ) :
. . . set up and c a l c u l a t e J . . .
return J
2. Test vnewton1() on the vector problem defined in v2.py starting from [1,1]. Compare this result with
what you obtained in the previous lab, where newton3() used complex arithmetic on f8.py, starting
from 1+1j.
3. Fill in the table below. One approach to these four calculations is simply to repeat the full set of
commands four times. But you may be interested to try using the for loop to simplify things like this:
x 0 a r r a y = np . a r r a y ( [ [ 1 . 0 , 1 . 0 ] , [ 1 . 0 , −1. 0] , [ 1 0 . 0 , 5 . 0 ] , [ 1 0 . 0 , 1. 0E−25] ] )
fo r x0 in x 0 a r r a y :
x , numits = vnewton1 ( v2F , v2J , x0 )
. . . compute e r r o r . . .
. . . print numIts and e r r o r norm . . .
Initial guess numIts norm(error)
[1,1] _______ _________
[1,-1] _______ _________
[10,5] _______ _________
[10,1.e-25] _______ _________
7 Exercise 3
At this point, you should be confident that vnewton1() is correctly programmed and yields the same results
as newton3() from the previous lab. In this exercise, you will apply vnewton1() to a simple nonlinear
problem.
3
Use Newton’s method to find the two intersection points of the parabola y = x
2 − x and the ellipse
x
2/16 + y
2 = 1.
1. Plot the parabola and ellipse together over the range −4 ≤ x ≤ 4. To do this, you have to rewrite the
ellipse as a formula for y; because you are taking the square root to do this, there will be two formulas,
one for the upper half of the ellipse, and one for the lower half. Thus you will need to make three
plt.plot() statements to show the whole picture. Please include this plot when you send me your
work.
2. In order to work with Newton’s method, it is necessary to think of x and y as x[0] and x[1], and to
write both equations as functions. Thus, you should now think of the first equation as
F [ 0 ] = x [ 0 ] ∗ ∗ 2 − x [ 0 ] − x [ 1 ]
with a similar reformulation of the second equation.
3. Write a file v3.py to compute the vector function F and its jacobian J.
def v3F ( x ) :
. . . compute F . . .
return F
def v3J ( x ) :
. . . compute J . . .
return J
If we can solve F(x) = 0, then we have found a pair of values corresponding to an intersection of the
parabola and the ellipse.
4. Starting from the initial vector guess [2,1], what is the intersection point that vnewton1() finds, and
how many iterations does it take?
5. Find the other intersection point by choosing a different initial guess. What initial guess did you use,
what is the intersection point, and how many iterations did it take?
8 Exercise 4
In Exercise 2, you saw that the guess [10,1.e-25] resulted in a relatively large number of iterations. While
not a failure, a large number of iterations is unusual.
In this exercise, we will look more carefully at the iterates in the poorly-converging case from Exercise 2.
Because we are interested in the sequence of iterates, we will generate a special version of vnewton1.py that
returns the full sequence of iterates. It will return the iterates as a matrix whose rows are the iterates. (One
would not normally return so much data from a function call, but we have a special need for this exercise.)
1. Make a copy of vnewton1.py and call it vnewton2.py. Change its signature line to the following
def vnewton2 ( func , d func , x0 , maxIts ) :
”””x , numIts , i t e r a t e s = vnewton2 ( func , d func , x0 , maxI ts )
. . . comments
”””
x = ? ? ?
numIts = ? ? ?
i t e r a t e s = ? ? ?
return x , numIts , i t e r a t e s
2. Just before the loop, add the following statement
n = len ( x )
i t e r a t e s = np . z e r o s ( [ 1 , n ] )
i t e r a t e s [ 0 , : ] = x
4
to initialize the variable iterates.
3. Add the following just after the new value of x has been computed.
x = x + dx
i t e r a t e s = np . i n s e r t ( i t e r a t e s , numIts , x , a x i s = 0 )
This statement appends x as a new row to the array iterates.
4. Test your modified function vnewton2() using the same v2.py from above, starting from the vector
[1,1]. You should get the same results as before for solution and number of iterations. Check that
the size of the matrix of iterates is (iterations+1)×2, that is, it has (numIts+1) rows and 2 columns.
9 Exercise 5
The study of the poorly-converging case continues with this exercise. The objective is to try to understand
what is happening.
1. Apply vnewton2() to the function defined in v2.py starting from the vector [10,1.e-25]. This is the
poorly-converging case.
2. Plot the sequence of iterates on the plane using the command
p l t . pl o t ( i t e r a t e s [ : , 0 ] , i t e r a t e s [ : , 1 ] , ’∗− ’ )
p l t . show ( )
It is hard to interpret this plot, but it is a place to start. You do not have to send me a copy of your
plot. Note that most of the iterates lie along the x-axis, with some quite large.
3. If it’s available to you, use the zoom feature to focus on the portion of the plot with horizontal extent
around [-20,20]. It is clear that most of the iterates appear to be on the x-axis. Please send me a copy
of this plot.
4. Look at the formula you used for the Jacobian in v2.py. Explain why, if the initial guess starts with
x1 = 0 exactly, then all subsequent iterates also have x1 = 0. You have used computation to help
formulate a hypothesis that can be proved.
5. You should be able to see what is happening. We used an initial guess with x1 ≈ 0, so for a long time,
the subsequent iterates stayed near the horizontal-axis. Since the root is at x1 = 3, it takes many
iterates before the x1 component of the iterate pulls away from the horizontal axis and approaches the
region of quadratic convergence near the root at (0, 3).
6. Use a semilog plot to see how the second component of the iterates grows in the vertical direction:
p l t . s emil o g y ( np . abs ( i t e r a t e s [ : , 1 ] ) )
p l t . show ( )
The plot shows the x1 component seems to grow exponentially (linearly on the semilog plot) with
seemingly random jumps superimposed. Please send me a copy of this plot.
You now have a rough idea of how this problem is behaving. It converges slowly because the initial guess
is not close enough to the root for quadratic convergence to be exhibited. The second component grows
exponentially, however, and eventually the iterates become close enough to the root to exhibit quadratic
convergence.
10 Nonlinear least squares
A common application of Newton’s method for vector functions is fitting a nonlinear curve by the method
of least squares. This application falls into the general topic of optimization where the extremum of some
function D : R
n → R is sought. If D is differentiable, then its extrema are given by the solution of the system
5
of equations ∂D/∂xk = 0 for k = 0, 1, . . . , n − 1, which we can regard as the vector equation F(x) = 0. The
value x which zeroes out the vector equation is the same value which will minimize the scalar function D.
The differential equation describing the motion of a weight attached to a damped spring without forcing
is
m
d
2v
dt2
+ c
dv
dt + k v = 0,
where v is the displacement of the weight from equilibrium, m is the mass of the weight, c is a constant
related to the damping of the spring, and k is the spring stiffness constant. The physics of the situation
indicate that m, c and k should be positive. Solutions to this differential equation are of the form
v(t) = e
−x0t
(x2 sin x1t + x3 cos x1t),
where x0 = c/(2m), x1 =
p
k − c
2/(4m2), and x2 and x3 are values depending on the position and velocity
of the weight at t = 0. One common practical problem (called the parameter identification problem) is to
estimate the values x0 . . . x3 by observing the motion of the spring at many instants of time.
After the observations have progressed for some time, you have a large number of pairs of values (tn, vn) for
n = 0, . . . , N−1. The question to be answered is, “What values of xk would best reproduce the observations?”
In other words, find the values of x = [x0, x1, x2, x3] that minimize the norm of the differences between the
formula and observations. Define
D(x) =
N
X−1
n=0
(vn − e
−x0tn (x2 sin x1tn + x3 cos x1tn))2
(3)
and we seek the minimum of D.
In order to solve this problem, it is best to note that when the minimum is achieved, the gradient vector
F = ∇D must be zero. The components Fk of the gradient can be written as
F0 =
∂D
∂x0
= 2
K
X−1
k=0
(vk − e
−x0tk
(x2 sin x1tk + x3 cos x1tk)) tke
−x0tk
(x2 sin x1tk + x3 cos x1tk)
F1 =
∂D
∂x1
= −2
K
X−1
k=0
(vk − e
−x0tk
(x2 sin x1tk + x3 cos x1tk)) e
−x0tk
(x2 cos x1tktk − x3 sin x1tktk) (4)
F2 =
∂D
∂x2
= −2
K
X−1
k=0
(vk − e
−x0tk
(x2 sin x1tk + x3 cos x1tk)) e
−x0tk
sin(x1tk)
F3 =
∂D
∂x3
= −2
K
X−1
k=0
(vk − e
−x0tk
(x2 sin x1tk + x3 cos x1tk)) e
−x0tk cos(x1tk)
To minimize D, we want the components of F to be zero. So we want to apply Newton’s method to
F (and not to D!). Thus, we will also need the derivatives of F. The function D is a scalar function of a
vector. Its gradient, F = ∇D, is a vector function. The gradient of F is the Jacobian matrix, and is actually
the second derivative of D. In the least squares context, J is called the “Hessian” matrix of D, but for our
Newton application, we think of J as the Jacobian of F.
11 Exercise 6
1. Download the file objective.py which defines a function objectiveF() and jacobian objectiveJ() for
the least squares problem.
6
2. For this exercise, set x to the vector that actually minimizes the objective function:
x = np . a r r a y ( [ 0. 1 5 , 2 . 0 , 1 . 0 , 3 . 0 ] )
3. At the minimizer, F(x) should be zero. Verify this by calling objectiveF(x).
4. Compute the Jacobian by calling objectiveJ(x).
5. Compute the determinant of J(x) to see that it is nonsingular there:
v al u e = np . l i n a l g . de t ( J )
6. Since the matrix J is the second derivative of the underlying scalar function D, it must be positive
definite and symmetric. Check that J(x) is symmetric or very nearly so, using a command like
e r r = np . l i n a l g . norm ( J − np . t r a n s p o s e ( J ) )
7. Check that J is positive definite by computing its eigenvalues, which should all be positive real values,
using the command
val , vec = np . l i n a l g . e i g ( J )
which returns the eigenvalues (a vector) and eigenvectors (in a matrix).
12 Exercise 7
In this exercise, you will see that Newton’s method applied to this system can require the convergence
neighborhood to be quite small.
The point of this exercise is to see how sensitive Newton’s method can be when the initial guess is
changed. Fill in the following table with either the number of iterations required or the word “failed,” using
vnewton1() and the indicated initial guess. Note that the first line is the correct solution. For this exercise,
restrict the number of iterations to be no more than 100.
Initial guess Number of iterations
[0.15, 2.0, 1.0, 3] ______________
[0.15, 2.0, 0.9, 3] ______________
[0.15, 2.0, 0.0, 3] ______________
[0.15, 2.0,-0.1, 3] ______________
[0.15, 2.0,-0.3, 3] ______________
[0.15, 2.0,-0.5, 3] ______________
[0.15, 2.0, 1.0, 4] ______________
[0.15, 2.0, 1.0, 5] ______________
[0.15, 2.0, 1.0, 6] ______________
[0.15, 2.0, 1.0, 7] ______________
[0.15, 1.99, 1.0, 3] ______________
[0.15, 1.97, 1.0, 3] ______________
[0.15, 1.95, 1.0, 3] ______________
[0.15, 1.93, 1.0, 3] ______________
[0.15, 1.91, 1.0, 3] ______________
[0.17, 2.0, 1.0, 3] ______________
[0.19, 2.0, 1.0, 3] ______________
[0.20, 2.0, 1.0, 3] ______________
7
You can see from the previous exercise that Newton can require a precise guess before it will converge.
Sometimes some iterate is not far from the ball of convergence, but the Newton step is so large that the next
iterate is ridiculous. In cases where the Newton step is too large, reducing the size of the step might make
it possible to get inside the ball of convergence, even with initial guesses far from the exact solution. This is
the strategy examined in the following section.
13 Damping
The exercise in the previous section suggests that Newton gets in trouble when its increment is too large.
One way to mitigate this problem is to “dampen” the iteration by putting a fractional factor on the iterate.
x
(n+1) = x
(k) − αJ(x
(k)
)
−1F(x
(k)
) (5)
where α is a number smaller than one. It should be clear from the convergence proofs you have seen for
Newton’s method that introducing the damping factor α destroys the quadratic convergence of the method.
This raises the question of stopping. In the current version of vnewton1(), you stop when norm(increment)
gets small enough, but if increment has been multiplied by alpha, then convergence could happen immediately if alpha is very small. It is important to make sure norm(increment) is not multiplied by alpha
before the test is done.
14 Exercise 8
In this exercise, we experiment with damping.
1. Starting from your vnewton1.py file, copy it to a new file named vnewton3.py, change its signature to
def vnewton3 ( func , d func , x0 , maxIts ) :
”””x , numIts = vnewton3 ( func , d func , x0 , maxI ts )
. . . comments
”””
return x , numIts
and change it to conform with Equation (5) with the fixed value α = 1/2. Don’t forget to change the
comments and the convergence criterion.
2. Apply vnewton3() to the objectiveF function. First solve the problem for the exact solution [ 0.15,
2.0, 1.0, 3.0 ]. Then increase the first component of x0 by 0.01 and solve again. Keep increasing x0
by 0.01 and solving, until the iteration diverges. What is the largest value of x0 for which vnewton3()
still converges? Your code might look something like this:
x0 = np . a r r a y ( [ 0. 1 5 , 2 . 0 , 1 . 0 , 3. 0 ] )
while ( True ) :
try :
x , numits = vnewton3 ( o b j e c ti v eF , o b j e c ti v e J , x0 )
print ( ” vnewton3 : x0 =” , x0 , ” numits =” , numits )
x0 [ 0 ] = x0 [ 0 ] + 0. 0 1
except Excep ti on :
print ( ” vnewton3 : x0 =” , x0 , ” ( Did not c on v e r ge ! ) ” )
break
Damping by a constant factor can improve the initial behavior of the iterates, but it destroys the quadratic
convergence of the method. Further, it is hard to guess what the damping factor should be. There are tricks
to dampen the iteration in such a way that when it starts to converge, the damping goes away (α → 1).
8
One such trick is to compute
∆x = −J(x
(k)
)
−1
f(x
(k)
)
α =
1
1 + βk∆xk
(6)
x
(k+1) = x
(k) + α∆x
where the value β = 10 is a conveniently chosen constant. You should be able to see how you might prove
that this damping strategy does not destroy the quadratic convergence rate, or, at least, allows a superlinear
rate.
Note that the expression in (6) is designed to keep the largest step below one tenth of the Newton step.
This is a very conservative strategy. Note also that another quantity could be placed in the denominator,
such as kfk, so long as it becomes zero at the solution.
15 Exercise 9
1. Make a copy of your vnewton3.py file as the new file vnewton4.py, and change it so that α is not the
constant 1/2, but is determined from Equation (6). Don’t forget to change the comments.
2. Apply vnewton4() to the objectiveF function. for the initial guess x0 set to [ 0.20, 2.0, 1.0,
3.0 ]. You should see that vnewton4() has a slightly larger ball of convergence than vnewton3(),
and converges much faster.
3. Using vnewton4(), repeatedly increase the first component of the starting point x0 by increments of
0.01, until the iteration diverges. What is the largest such value for which vnewton4() is still able to
solve the problem?
Another strategy is to choose α so that the objective function almost always decreases on each iteration.
This strategy can be useful when the the previous approach is too slow. One way to implement this strategy
is to add an additional loop inside the Newton iteration loop. Begin with a full Newton step, but check that
this step reduces the objective function. If it does not, halve the step size and try again. Keep halving the
step size until either the objective function is reduced or some fixed maximum number of halvings (10, for
example) have been tried. This can be written as an algorithm:
1. On Newton iteration k + 1, start with α = 1.
2. Compute a trial update:
x˜
(`) = x
(k) − αJ(x
(k)
)
−1
f(x
(k)
)
3. If |f(x˜
(`)
)| ≤ |f(x
(k)
)| or if α < 1/1024, then set x
(k+1) = x˜
(`) and continue with the next Newton
iteration,
otherwise replace α with α/2 and return to Step 2.
The limit α < 1/1024 represents the (arbitrary) choice of a maximum of 10 halvings.
16 Exercise 10
Here is the outline of a revised Newton code that we will need for the next exercise.
def vnewton5 ( func , d func , x0 , maxIts ) :
”””x , numIts = vnewton5 ( func , d func , x0 , maxI ts )
. . . comments
”””
import numpy a s np
9
EPSILON = 5. 0 e−05
x = x0 . copy ( )
fo r numIts in range ( 1 , maxIts + 1 ) :
f = func ( x )
fnrm = np . l i n a l g . norm ( f ) # Need t h i s now f o r compar ison
J = d func ( x )
dx = − np . l i n a l g . s o l v e ( J , f )
alph a = 1. 0 # i n i t i a l i z e al p h a
while ( True ) : # Loop s e a r c h i n g f o r good al p h a
x2 = ? ? ?
f 2 = ? ? ?
f2nrm = ? ? ?
i f ( f2nrm < fnrm ) :
x = np . copy ( x2 ) # x = x2 doesn ’ t work f o r py thon v e c t o r s !
break # s u c c e s s !
alph a = alph a / 2. 0 # maybe t r y sm a l l e r al p h a ?
i f ( alph a <= 1. 0 / 1 0 2 4. 0 ) :
ra is e Excep ti on ( ” vnewton5 : damping f a i l e d . ” )
e r r = np . l i n a l g . norm ( dx )
i f ( e r r <= EPSILON ) :
return x , numIts
ra is e Excep ti on ( ” vnewton5 : i t e r a t i o n did not c on ve r ge . ” )
Of course, it is necessary for you to replace the ”???” text by the appropriate formulas.
1. Using vnewton5(), to the nearest 0.01, how large can the first component of the initial guess get before
the iteration diverges? (Leave the other three values at their correct values.)
Warning: You may find this strategy is not better than the previous one!
17 Continuation or homotopy methods
As can be seen in the previous few exercises, there are ways to improve the radius of convergence of Newton’s
method. For some problems, such as the curve-fitting problem above, they just don’t help enough. There
is a class of methods called continuation or homotopy methods that can be used to find good initial guesses
for Newton’s method.
In the previous section, we were concerned with finding the value x which minimizes a scalar objective
D1(x). Suppose there is another, much simpler, objective D2(x), whose minimizer z is easy to find using
Newton’s method. One possible choice would be D2(x) = kx − zk
2
, where z might be arbitrary, or better,
chosen as any approximate guess for the minimizer of D1(x) . For 0 ≤ p ≤ 1, consider the new objective
function
G(x, p) = pD1(x) + (1 − p)D2(x). (7)
When p = 0, the function G reduces to D2 and the corresponding vector equation F2(x) = ∂D2
∂x is easy
for Newton’s method to solve with any initial guess. (Naturally, our starting guess here will be z!). When
p = 1, G is equal to D1, and we are back to solving F(x) = ∂D1
∂x . The only new thing is that our starting
guess for this iteration is likely to be much better than z, or any other guess we might make in advance.
That’s because we plan to solve a sequence of problems, turning up the value of p gradually, so that Newton
starts on the easy problem, and works its way towards the hard problem, each time using the solution of the
previous step as its starting point.
10
18 Exercise 11
1. We suppose that z is a fixed vector of size 4, which is a rough approximate minimizer of D1(x). Define
D2(x) = X
4
k=1
(xk − zk)
2
.
Its gradient is
F2k =
∂D2
∂xk
and the Jacobian matrix is given by
J2k,` =
∂F2
∂x`
The following code, easy objective.py, computes these quantities: similar to objective.py.
def e a s y o b j e c ti v e D ( x , z ) :
import numpy a s np
D2 = np .sum ( ( x − z ) ∗∗2 )
return D2
def e a s y o b j e c t i v e F ( x , z ) :
import numpy a s np
n = len ( x )
F2 = np . z e r o s ( n )
F2 = ? ? ?
return F2
def e a s y o b j e c t i v e J ( x , z ) :
import numpy a s np
n = len ( x )
J2 = np . z e r o s ( ( n , n ) )
fo r k in range ( 0 , n ) :
J2 [ k , k ] = 2. 0
return J2
Make a copy of this file, and fill in the missing calculations for F2.
The chosen value for z is problem-related and amounts to a vague approximation to the final solution.
2. Copy the following text to a file named homotopy.py, and fill in the missing details indicated by
question marks:
def homotopyD ( x , p , z ) :
from e a s y o b j e c t i v e import e a s y o b j e c ti v e D
from o b j e c t i v e import o b j e c ti v eD
D1 = o b j e c ti v eD ( x )
D2 = e a s y o b j e c ti v e D ( x , z )
D = p ∗ D1 + ( 1. 0 − p ) ∗ D2
return D
def homotopyF ( x , p , z ) :
from e a s y o b j e c t i v e import e a s y o b j e c t i v e F
from o b j e c t i v e import o b j e c ti v e F
F1 = ?
F2 = ?
F = ?
return F
def homotopyJ ( x , p , z ) :
?
return J
11
3. Create the following code named homotopy run.py.
def homotopy run ( z )
from homotopy import homotopyF , homotopyJ
from vnewton1 import vnewton1
import numpy a s np
x0 = z
STEPS = 1000
MAX ITERS = 100
p = 0. 0
print ( ” p n x [ 0 ] x [ 1 ] x [ 2 ] x [ 3 ] ” )
fo r k in range ( 0 , STEPS + 1 ) :
p = k / STEPS
x , n = vnewton1 ( lambda X: homotopyF ( X, p , z ) , \
lambda X: homotopyJ ( X, p , z ) , x0 , MAX ITERS )
i f ( k == STEPS or ( k % 50 ) == 0 ) :
print ( p , n , x )
x0 = np . copy ( x ) # Here i s our s t a r t i n g p o i n t f o r nex t Newton i t e r a t i o n .
4. Explain why the expression
lambda X: homotopyF ( X, p , z )
is necessary for this problem, in order to successfully call vnewton1(). You can look up the term
python lambda for information. Why must this construct be used instead of simply using homotopyF?
5. Can you start from z=[0.24, 2, 1, 3] and reach the solution at p=1? Do you get the same solution
as in Exercises 5, 6, 8 and 9?
6. Can you start from z=[0,2,1,3] and reach the solution?
7. Can you start from z=[-0.5,2,1,3] and reach the solution?
8. Start from z=[.25,2,1,3] and solve the problem. Repeatedly increase the first component of z by
0.05 and solve the problem. At what value of the first component of z does this process fail? You
could automate this experiment as follows:
z = np . a r r a y ( [ 0. 2 5 , 2 . 0 , 1 . 0 , 3. 0 ] )
try :
while ( True ) :
homotopy run ( z )
z [ 0 ] = z [ 0 ] + 0. 0 5
except :
print ( ” homotopy run f a i l e d f o r z [ 0 ] = ” , z [ 0 ] )
12