EGR 103/Spring 2023/Lab 7

7.4.1 Group Gradescope Problems

Be sure to first understand the code in Chapra Figure 4.2!

7.4.2 Individual Lab Report

7.4.2.1 Based on Chapra Problem 4.25

This is mainly meant to illuminate issues with large angle representations for the approximations.

7.4.2.2 Based on Chapra Problem 4.28

The key here is that you need to use a loop to calculate each value of your approximation - you cannot calculate them all at once!

7.4.2.3 Based on Chapra Problem 14.34

This problem shows how Monte Carlo methods can be used to see how variations in parameters might impact output values.

7.4.2.4 Based on Chapra Problem 14.34

This problem shows how using more points in Monte Carlo methods can be used to better see how variations in parameters might impact output values.

7.4.2.5 Finding Roots Using Newton-Raphson

For this one, you will be making changes to either the original or extended version of Chapra 4.2 in order to use a mapping to estimate roots of a function. It is very similar to using the Newton method mapping to estimate values of a square root. Note - though the mapping for the lab assignment comes from a method called Newton-Raphson, its origin is different from the Newton method discussed in class for finding a square root.

Basin Plotter

You will be using the code below to see how various initial guesses evolve into final estimates for the roots of the equation. The sections of code are described below:

• """
  @author: DukeEgr93
  """
  
  # %% Initialize workspace
  import numpy as np
  from poly_root import calc_root
  import matplotlib.pyplot as plt
  

  In addition to importing numpy and matplotlib.pyplot, this will import the calc_root function from your poly_root.py file. Make sure your file is called poly_root.py and make sure your function is called calc_root (not iter_meth).

• # %% Generate guesses and start lists
  xi = np.linspace(0, 5, 1000)
  rootlist = []
  iterlist = []
  

• # %% Run program in a loop and store roots and number of iterations
  for init in xi:
      out = calc_root(init, 1e-12, 1000)
      rootlist += [out[0]]
      iterlist += [out[2]]
  

  We are going to see how 1000 different initial guesses between 0 and 5 evolve to final estimates for the roots. xi will be a set of linearly spaced initial guesses for the root $$x$$. We want to track both where the estimates are as well as how long it takes to get there. To do that, we will be appending the estimates and iteration counts to a list, which means we need to start with two empty lists. The program then runs a loop that looks at each entry in xi, passes it to your calc_root function with a very, very small stopping error and a maximum iteration count of 1000. You function will return the estimate, the final error estimate, the number of iterations it took to get to that estimate, and - if you used the extended version of the code - lists with all the estimates and error estimates. The only things we want to keep track of are the final estimate and final error estimate, so we slice those from the output out and append them to the appropriate lists using the += operator. Note that the values we are appending using += must be in a list for += to work; trying to use += with a list and an int or a float will give an error:

In [1]: a = [1, 2, 3, 4]

In [2]: a+=5
Traceback (most recent call last):

  File "<ipython-input-4-1a160a0ff440>", line 1, in <module>
    a+=5

TypeError: 'int' object is not iterable


In [3]: a+=[5]

In [4]: a
Out[4]: [1, 2, 3, 4, 5]

• # %% Make figure with function and map
  fig0, ax0 = plt.subplots(2, 1, num=0, clear=True)
  fig0.set_size_inches(6, 8, forward=True)
  ax0[0].plot(xi, xi ** 3 - 7 * xi ** 2 + 14 * xi - 8, "k-")
  ax0[0].grid(True)
  ax0[0].set(title="Function", ylabel="$f(x)$", xlabel="$x$")
  

  This will set up figure window 0 to have two rows and one column of subplots (i.e. one subplot over another). It then changes the size of the figure window, plots the function whose roots we are looking to find in the top subplot, adds a grid, and labels the plot.

• ax0[1].plot(
      xi,
      (2 * xi ** 3 - 7 * xi ** 2 + 8) / (3 * xi ** 2 - 14 * xi + 14),
      "r-",
      label="map",
  )
  ax0[1].plot(xi, xi, "k:", label="new=old line")
  ax0[1].set_ylim([-10, 10])
  ax0[1].set(title="Map to Find Roots", ylabel="$x_{k+1}$", xlabel="$x_k$")
  ax0[1].legend(loc=0)
  
  fig0.tight_layout()
  fig0.savefig("RootPlot0.png")
  

  This code plots in the bottom subplot. Specifically, it plots the map we are using for our iterative method. It also plots the line $$y=x$$; from that, you should see that the map maps to the same value when we are at one of the roots of the equation. That is to say, if you look at the value of the map at $$x=1$$, $$x=2$$, and $$x=4$$, the map is equal to 1, 2, and 4, respectively. At all other locations, the map would move us to some other value for the estimate.

• # %% Make figure with roots and iteration counts
  fig1, ax1 = plt.subplots(2, 1, num=1, clear=True)
  fig1.set_size_inches(6, 8, forward=True)
  
  ax1[0].plot(xi, rootlist, "k.")
  ax1[0].set(title="Root Estimate", ylabel="Root", xlabel="Initial Guess")
  
  ax1[1].plot(xi, iterlist, "k.")
  ax1[1].set(title="Iteration Count", ylabel="Iterations", xlabel="Initial Guess")
  
  fig1.tight_layout()
  fig1.savefig("RootPlot1.png")
  

  This figure, which once again has two subplots, will let us look at where the estimate is and how long it took to get there.

• # %% Visualize roots and interation counts differently
  fig2, ax2 = plt.subplots(2, 1, num=2, clear=True)
  fig2.set_size_inches(6, 8, forward=True)
  rli = ax2[0].imshow(np.array([rootlist]), aspect="auto", extent=(xi[0], xi[-1], 0, 1))
  ax2[0].set_yticklabels([])
  fig2.colorbar(rli, ax=ax2[0])
  ax2[0].set(title="Root Estimate", xlabel="Initial Guess")
  

  This figure will also let us look at where the estimate is and how long it took to get there. Instead of making graphs, however, this makes images. The imshow command will take an array and allow us to look at the values by assigning the values in the array to colors on a colormap. The built-in color map (called "viridis" - see more information about colormaps at Choosing Colormaps in Matplotlib at matplotlib.org) goes from dark purple for the lowest values through teal and green to bright yellow for the highest values. Note in this image that the values congregate around the colors representing 1 (dark purple), 2 (medium teal), and 4 (bright yellow).

• tli = ax2[1].imshow(np.array([iterlist]), aspect="auto", extent=(xi[0], xi[-1], 0, 1))
  ax2[1].set_yticklabels([])
  fig2.colorbar(tli, ax=ax2[1])
  ax2[1].set(title="Iteration Count", xlabel="Initial Guess")
   
  fig2.tight_layout()
  fig2.savefig("RootPlot2.png")
  

  This subplot is also an image, but this time representing how long it took to get to a particular estimate. The colorbar is the same, but the values represented by each color are different. In this case, the dark purple represents initial conditions that got to their final estimates relatively quickly while the bright yellow represents initial conditions that needed to evolve for a while before settling.