EGR 103/Fall 2022/Lab 8

8.4.1 APT

Note that these are due a week later than the lab!

8.4.2 Group Gradescope Problems

8.4.2.1 Geometric Progressions

Among other things, this problem looks at input validation. You will definitely want to look at the fib_stuff_if.py file in the Lectures/Lec 17 folder on Box. Which is to say,

fd = {}
def fibm(n, f1=1, f2=1):
    if not isinstance(n, int):
        "Error: n must be an integer" 
        return -1
    if n<=0:
        "Error: n must be 1 or greater"
        return -2
    fd[1] = f1
    fd[2] = f2
    if n not in fd:
        fd[n] = fibm(n-1, f1, f2) + fibm(n-2, f1, f2)
    return fd[n]
    
if __name__ == "__main__":
    print(fibm(1.1))
    
    print(fibm(0))
    
    print(fibm(10))

Note that the assignment is different in two ways: you need to check if something is either an int or a float, and you need to check multiple variables. For the former, you can put the variable types that are good in a tuple for the second argument (that is, isinstance(x, (int, float)) will return True if x is either an int or a float). For the latter, well, that's up to you!

You will also want to make sure you think about how you will be checking for valid inputs based on the third requirement and then write and test the code for that. Remember that by the time you are checking to see if the series works, you have already checked and made sure you have positive numerical values for all three entries.

8.4.2.2 Inspired by P&E 2.3.8

The solution to this will require understanding recursion, which requires understanding recursion. For more information beyond what we covered in class, read Chapter 16 (through 16.4) in the Runestone Academy Book. Your program must use recursion to get full credit. Be sure to look at the test case (using 827) and start by getting your program to figure out that the sum of those digits is 17. Then figure out how to get the code to run itself on the number 17 to calculate (and return) 8.

This example page on Fibonacci Numbers also looks at recursion (as well as other methods). This other example page on Python Using Recursion looks at recursion as well as building a list during recursion to reduce the amount of recursion needed. The latter will not be useful for the digit sum problem but is cool nonetheless.

8.4.3 Individual Lab Report

8.4.3.1 Based on Chapra Problem 4.25

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.

8.4.3.2 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.