EGR 103/Spring 2022/Lab 5

Errors / Additions

None yet!

5.1 Introduction

The main purpose of this lab will be to look at debugging, dictionaries, and using Chapra Figure 4.2.

5.2 Resources

See main EGR 103 page for links to these

5.3 Getting Started

Beginning of Lab

5.4 Debugging Introduction

Learn more about debugging by reading Python:Debugging page.

5.5 Group Connect Problems

The second Connect assignment is active. For Chapra 3.4, note that you will be using your function as part of the lab assignment and also note that the table in the Connect assignment may have different values from the book. Use the values from the Connect assignment in all cases. This Connect assignment can be done in small groups, but remember, Connect will give different values to different people so be sure to have your values in the code you write since the autograder will be using your values.

5.6 Individual Lab Assignment

5.6.1

Note that the skeleton of the code is already given. The chapra_0304.py file should not do any of the graphing - it should only contain a function that takes two arguments and then returns an array of temperatures.

5.6.2

For this one, you will be making changes to an extended version of the Chapra 4.2 code in order to use accumulation to estimate values of cosine. It is very similar to using an accumulation to estimate values of the exponential.

5.6.3

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.