You do not need to do the next command when using the IPython notebook on SageMathCloud, but it may be needed if you open this notebook on another system where the defaults are set differently. It does "from pylab import *" so numpy and plotting commands are available, and sets it up for plots to appear in the notebook webpage
In [1]:
%pylab inline
In [1]:
def newton(fvals, x0, tol):
xk = x0
kmax = 30 # maximum number of iterations
print " k xk f(xk)"
for k in range(kmax):
fxk, fpxk = fvals(xk) # evaluate f(x) and f'(x)
print "%4i %22.15f %22.15f" % (k, xk, fxk)
if abs(fxk) < tol:
break #leave the loop if tolerance satisfied
xk = xk - fxk / fpxk # update xk using Newton's method
return xk
This next function takes as input some value xk (representing $x_k$ in an iteration of Newton's method) and returns $x_{k+1}$. It also plots the function along with the tangent line at xk that is used to define the next iterate. It plots the function and tangent line over an interval including the two points and adds a legend that contains the two values.
In [2]:
def newton_step_plot(xk, fvals):
fxk, fpxk = fvals(xk)
xkp1 = xk - fxk / fpxk # calculate x_{k+1} from x_k
xdiff = abs(xkp1 - xk)
xlower = min(xk,xkp1) - xdiff
xupper = max(xk,xkp1) + xdiff
x = linspace(xlower,xupper,500)
fx, fpx = fvals(x)
tangent_line = fxk + (x-xk)*fpxk
plot(x, fx, 'b')
plot(x, tangent_line, 'r')
plot([xk], [fxk], 'ro', label="x_k = %g" % xk)
plot([xkp1], [0.], 'go', label="x_{k+1} = %g" % xkp1)
plot([xlower, xupper], [0,0], 'k')
legend(loc='best') # uses the labels specified in the plot commands above
axis('tight')
return xkp1
In [6]:
def fvals(x):
f = (x**2 - 2)*exp(-0.1*x**2) + 0.5
fprime = (2*x -0.2*x*(x**2-2)) * exp(-0.1*x**2)
return f, fprime
# Plot the function over a suitable range to get an idea of where the zeros are:
x = linspace(-10,10,1000)
fx, fpx = fvals(x)
plot(x,fx)
plot(x,zeros(len(x)), 'k') # x-axis
#ylim(-2,0.5)
Out[6]:
In [7]:
xk = newton(fvals, 1., 1e-10)
fxk, fpxk = fvals(xk)
print "At xk = %22.15e, f(xk) = %22.15e " % (xk,fxk)
In [17]:
def fvals(x):
f = 2 + 20*sin(x**2) - exp(x)
fprime = 40*x*cos(x**2) - exp(x)
return f, fprime
# Plot the function over a suitable range to get an idea of where the zeros are:
x = linspace(-6,4,1000)
fx, fpx = fvals(x)
plot(x,fx)
plot(x,zeros(len(x)), 'k') # x-axis
#ylim(-2,0.5)
Out[17]:
Zoom in to look for good initial guesses:
In [19]:
x = linspace(0,4,1000)
fx, fpx = fvals(x)
plot(x,fx)
plot(x,zeros(len(x)), 'k') # x-axis
figure()
plot(x,fx)
plot(x,zeros(len(x)), 'k') # x-axis
xlim([2.5, 3.0])
ylim([-10,10])
Out[19]:
Try various initial guesses to try to find the three positive zeros:
In [12]:
x0 = 1. # note this converges to a an xk < 0
xk = newton(fvals, x0, 1e-10)
fxk, fpxk = fvals(xk)
print "At xk = %22.15e, f(xk) = %22.15e " % (xk,fxk)
In [13]:
x0 = 2.
xk = newton(fvals, x0, 1e-10)
fxk, fpxk = fvals(xk)
print "At xk = %22.15e, f(xk) = %22.15e " % (xk,fxk)
In [14]:
x0 = 2.5
xk = newton(fvals, x0, 1e-10)
fxk, fpxk = fvals(xk)
print "At xk = %22.15e, f(xk) = %22.15e " % (xk,fxk)
In [15]:
x0 = 3.
xk = newton(fvals, x0, 1e-10)
fxk, fpxk = fvals(xk)
print "At xk = %22.15e, f(xk) = %22.15e " % (xk,fxk)
Find the two points where the circle and line intersect in the plot below. The equation of the circle is $(x-1)^2 + y^2 = 4$ and the curve $y = \cos(x)$. Combine these equations by plugging $\cos(x)$ in for $y$ in the equation for the circle. This gives a single nonlinear equation in $x$ that you can solve using Newton's method.
Try to find both intersections and then replot the circle with these two points as red dots.
In [63]:
theta = linspace(0, 2*pi, 1000)
x = 1 + 2*cos(theta)
y = 2*sin(theta)
plot(x,y,'b')
x = linspace(-2,4,1000)
y = cos(x)
plot(x,y,'b')
axis('scaled') # so circle looks circular
Out[63]:
In [ ]: