<jupyter_start><jupyter_text>**Section1-Exercise3**Implement the method of false position with appropriate stopping condition.<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>Take $f(x) = x^3 + 2x^2 - 3x -1$ as example.<jupyter_code>def f(x):
  return (x * x * x + 2 * x * x - 3 * x - 1)<jupyter_output><empty_output><jupyter_text>Algorithm of the method of false position：1. if $f(a)f(b) < 0$ with $a \neq b$, then we can continue2. iterate $p_n = b_n - f(b_n)\frac{b_n-a_n}{f(b_n) - f(a_n)}$3. from the lecture note, we have $|e_n| ≈ \frac{|p_n-p_{n-1}|^2}{|p_n-2p_{n-1}+p_{n-2}|}$, and so terminating the iteration when  $\frac{|p_n-p_{n-1}|^2}{|p_n-2p_{n-1}+p_{n-2}|} < ϵ = tol$ NOTE：Since the method of false position and the bisection method are both linear convergence, the former is sometimes faster and sometmes slower than the latter. In the following code, I set $10 \times⌈ log_2\frac{|b-a|}{tol} ⌉$(10 times than the max iteration number of the bisection method) to be the max iteration number to avoid infinite loops.<jupyter_code># define function of the method of false position
# input : 
#   a, b : the endpoints of the initial interval
#   tol : convergence tolerence
# output :
#   an, bn : the endpoints of the final interval
#   pn : approximated root of the given equation f(x) = 0

def FalsePosition(a, b, tol): # assume that a < b
  fa, fb = f(a), f(b)

  # check if f(a)f(b) < 0 is satisfied or not
  if fa * fb >= 0:
    print('Need to implement the method with condition f(a)f(b) < 0')
    return a, b, -1
  
  p0, p1, p2 = a, b, a # set initial value of p
  max_iter = np.ceil(np.log2(abs(b-a)/tol)).astype(int) * 10 # to avoid infinite loops, give a value as limitation
  print('max iteration:', max_iter)

  for i in range(1, max_iter + 1):
    # use method of false position to iterate p
    p2 = b - fb * (b - a) / (fb - fa)
    fp = f(p2)
    # check if the new p is root or not
    if fp == 0: # case 1: p is a root
      break
    elif fp * fa < 0: # case 2: p is not a root, but a root exists between a and p
      b = p2
    else: # case 3: p is not a root, but a root exists between p and b
      a = p2
    fa, fb = f(a), f(b) # update f(a) and f(b)
    print('i =', '%3d' % i, '  a =', '%.9e' % a, '  b =','%.9e' % b, '  p =', '%.9e' % p2, '  |f(p)| = ', '%.4e' % abs(fp))

    # if there are at least 3 terms of pn, then we can check the stop condition
    if i >= 3:
      approx_e = (p2 - p1) * (p2 - p1) / abs(p2 - 2 * p1 + p0)
      print('          p0=', '%.9e' % p0,'  p1=', '%.9e' % p1,'  p2=', '%.9e' % p2, '  approx_e =', '%.4e' % approx_e)
      p0, p1 = p1, p2
      if approx_e < tol:
        break
  
  print('result:')
  print('i =', '%3d' % (i+1), '  a =', '%.9e' % a, '  b =','%.9e' % b, '  p = ', '%.9e' % p2, '  |f(p)| = ', '%.4e' % abs(fp), '   approx_e = ', '%.4e' % approx_e)
  return a, b, p2
# input: (a, b) = (1, 2)
an, bn, pn = FalsePosition(1, 2, 1e-9)
# input: (a, b) = (-3, -2)
an, bn, pn = FalsePosition(-3, -2, 1e-9)<jupyter_output>max iteration: 300
i =   1   a = -3.000000000e+00   b = -2.833333333e+00   p = -2.833333333e+00   |f(p)| =  8.1019e-01
i =   2   a = -3.000000000e+00   b = -2.907928389e+00   p = -2.907928389e+00   |f(p)| =  4.6300e-02
i =   3   a = -3.000000000e+00   b = -2.912002629e+00   p = -2.912002629e+00   |f(p)| =  2.4451e-03
          p0= -3.000000000e+00   p1= -2.000000000e+00   p2= -2.912002629e+00   approx_e = 4.3501e-01
i =   4   a = -3.000000000e+00   b = -2.912217267e+00   p = -2.912217267e+00   |f(p)| =  1.2858e-04
          p0= -2.000000000e+00   p1= -2.912002629e+00   p2= -2.912217267e+00   approx_e = 5.0526e-08
i =   5   a = -3.000000000e+00   b = -2.912228552e+00   p = -2.912228552e+00   |f(p)| =  6.7599e-06
          p0= -2.912002629e+00   p1= -2.912217267e+00   p2= -2.912228552e+00   approx_e = 6.2632e-07
i =   6   a = -3.000000000e+00   b = -2.912229146e+00   p = -2.912229146e+00   |f(p)| =  3.5539e-07
          p0= -2.912217267e+00   p1= -2.912228552e+00   p2= -2.912229146e+00  [...]<jupyter_text>Compare with the bisection method：<jupyter_code># define function of the bisection method
# input : 
#   a, b : the endpoints of the initial interval
#   tol : convergence tolerence
# output :
#   an, bn : the endpoints of the final interval
#   pn : approximated root of the given equation f(x) = 0
def bisection(a, b, tol):
  fa, fb = f(a), f(b)
  if fa * fb >= 0:
    print('Need to implement the method with condition f(a)f(b) < 0')
    return a, b, -1
  
  max_iter = np.ceil(np.log2(abs(b-a)/tol)).astype(int)
  print('max iteration:', max_iter)

  for i in range(1, max_iter + 1):
    p = 0.5 * (a + b)
    fp = f(p)
    print('i =', '%3d' % i, '  p = ', '%.16e' % p, '  |f(p)| = ', '%.4e' % abs(fp))

    if fp == 0:
      break
    elif fa * fp < 0:  
      b = p
      fb = fp
    else:
      a = p
      fa = fp
  
  print('result:')
  print('i =', '%3d' % (i+1), '  a =', '%.9e' % a, '  b =','%.9e' % b, '  p = ', '%.9e' % p, '  |f(p)| = ', '%.4e' % abs(fp))
  return a, b, p
# input: (a, b) = (1, 2)
an, bn, pn = bisection(1, 2, 1e-9)
# input: (a, b) = (-3, -2)
an, bn, pn = bisection(-3, -2, 1e-9)<jupyter_output>max iteration: 30
i =   1   p =  -2.5000000000000000e+00   |f(p)| =  3.3750e+00
i =   2   p =  -2.7500000000000000e+00   |f(p)| =  1.5781e+00
i =   3   p =  -2.8750000000000000e+00   |f(p)| =  3.9258e-01
i =   4   p =  -2.9375000000000000e+00   |f(p)| =  2.7710e-01
i =   5   p =  -2.9062500000000000e+00   |f(p)| =  6.4301e-02
i =   6   p =  -2.9218750000000000e+00   |f(p)| =  1.0475e-01
i =   7   p =  -2.9140625000000000e+00   |f(p)| =  1.9812e-02
i =   8   p =  -2.9101562500000000e+00   |f(p)| =  2.2347e-02
i =   9   p =  -2.9121093750000000e+00   |f(p)| =  1.2931e-03
i =  10   p =  -2.9130859375000000e+00   |f(p)| =  9.2531e-03
i =  11   p =  -2.9125976562500000e+00   |f(p)| =  3.9784e-03
i =  12   p =  -2.9123535156250000e+00   |f(p)| =  1.3422e-03
i =  13   p =  -2.9122314453125000e+00   |f(p)| =  2.4469e-05
i =  14   p =  -2.9121704101562500e+00   |f(p)| =  6.3434e-04
i =  15   p =  -2.9122009277343750e+00   |f(p)| =  3.0494e-04
i =  16   p =  -2.9122161865234375e+00   |f(p)| =  1[...]
{'fraction_non_alphanumeric': 0.18521514073115497, 'fraction_numerical': 0.164024587512132, 'mean_word_length': 2.483380281690141, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 20, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '58745', 'license': 'MIT', 'n_tokens_mistral': 3399, 'n_tokens_neox': 2845, 'n_words': 766, 'path': 'sec1exercise03.ipynb', 'repo_name': 'teshenglin/computational_mathematics', 'size': '6182'}