Rootfinding Demonstrations

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.set_printoptions(legacy='1.25')

Examples

A Simple Example

$f(x) = x^3 - 4x^2 - 5$

Find $x$ so that $f(x) = 0$.

def f(x):
    return x**3 - 4*x**2 - 5.

xx = np.linspace(-1, 5, 2000)
plt.figure(figsize=(4,4))
plt.plot(xx, f(xx));
plt.grid('on');
plt.axhline(0, color='k', lw=0.75);

f(4.275)

0.02579687500001171

A Concrete Example

Suppose the grades in a class follow a normal distribution with mean 72% and standard deviation 10%. What grade would beat 90% of the class?

from scipy.special import erf

def gaussian(x, mu, sig):
    g = 1.0 / (np.sqrt(2.0 * np.pi) * sig) * np.exp(-((x - mu) / sig)**2 / 2)
    return g

def myerf(x, mu, sig):
    return 0.5*( 1 + erf((x-mu)/sig/np.sqrt(2)) )

x = np.linspace(0,1,200);
mu, sig = 0.72, 0.1
f = gaussian(x, mu, sig)
plt.figure(figsize=(6,2))
plt.plot(x, f, color='gray', lw=1);

x = np.linspace(0,1,200);
mu, sig = 0.72, 0.1
plt.figure(figsize=(6,2))
plt.plot(x, myerf(x, mu, sig)); plt.grid('on');
plt.plot(x,[0.9]*len(x),'r--');

myerf(0.847, mu, sig)

0.8979576849251809

Fixed Point Example 1

def phi(x):
    # To solve
    #    f(x) = cos(x) - x = 0
    #  ==>  x = cos(x)
    return np.cos(x)

x = np.linspace(0,1.5,200);
plt.figure(figsize=(4,4))
plt.plot(x, phi(x)); plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');

def FixedPointIteration(x0, phi, n_iters=10): 
    x_history = []
    phi_history = []
    for idx in range(n_iters):
        # Plotting code
        x_history.append(x0)
        phi_history.append(phi(x0))

        x0 = phi(x0)  # <- Fixed-point iteration

        # Plotting code
        x_history.append(x0)
        phi_history.append(x0)
        print(x0)
    return x_history, phi_history

xh, phih = FixedPointIteration(0.25, phi)

0.9689124217106447
0.5661963244509128
0.8439474651598078
0.6645181803460105
0.7872140049933128
0.7058216222217144
0.7610788557592939
0.7240923421872705
0.7491010136260781
0.7323013570816245

plt.figure(figsize=(8,8)) #
plt.plot(x, phi(x)); plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go')

[<matplotlib.lines.Line2D at 0x15ce29910>]

Fixed Point Example 2

def phi(x):
    # f(x) = - (x^2 - 2) = 0
    #  ==>  x = x - (x^2 - 2)
    return x - (x**2 - 2.)

x = np.linspace(0,2,200);
plt.figure(figsize=(4,4))
plt.plot(x, phi(x));
plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
#plt.axis([x[0],x[-1],x[0],x[-1]])
plt.show()

xh, phih = FixedPointIteration(0.25, phi, n_iters=200)

2.1875
-0.59765625
1.0451507568359375
1.9528106523212045
0.13934120850203602
2.119925236115228
-0.3741577706029773
1.4858481920944326
1.2781033421441386
1.6445551889441217
0.9399934194610857
2.056405790830941
-0.17239898573208645
1.7978796039864613
0.5655085335559464
2.2457086320313495
-0.7974986279487659
0.5664973104710698
2.2455781077001142
-0.7970429300819117
0.567679637524529
2.2454194666645484
-0.7964891146115565
0.5691159756937423
2.2452229819039022
-0.7958032565655486
0.570893920274119
2.2449740520681667
-0.794934442391197
0.5731447899089996
2.2446498397091683
-0.793803063197227
0.5760736336614722
2.24421280226154
-0.792278299573054
0.5800167964525762
2.243597312285467
-0.790131587409105
0.5855604871692628
2.2426794030353583
-0.7869315017636729
0.5938073097682977
2.2412001886340347
-0.7817780968991981
0.6070449103094699
2.238541387176838
-0.7725261549267639
0.6306771850273056
2.232923473313339
-0.7530237643603677
0.6799314459481738
2.2176246747589996
-0.700234523340959
0.8094370889805009
2.154248687963274
-0.48653872162821266
1.276741350728172
1.646672874068975
0.9351413198743967
2.060652031737968
-0.18563476416784752
1.7799049701645
0.6118432673482102
2.2374910835488766
-0.7688752654118494
0.6399555608260088
2.2304124409938773
-0.744327215946389
0.7016497796551087
2.2093373663650464
-0.6718342320517925
0.8768045325915859
2.1080183442184364
-0.33572299534300143
1.5515670750549215
1.1442066866604372
1.8349977448619812
0.4677810212134246
2.24896193740595
-0.8088678584947746
0.5368649289993028
2.2486409770098765
-0.8077452664780553
0.5398023180042402
2.2484157754814893
-0.8069577239525376
0.5418615078008024
2.2482476141646437
-0.806369720432369
0.5433981535374541
2.2481166002695394
-0.8059116481379327
0.5445947672576681
2.2480113067332343
-0.8055435284672292
0.5455560952773371
2.247924642183082
-0.8052405547508554
0.5463470942336792
2.2478519468560947
-0.8049864281286405
0.5470104224000525
2.247790020185769
-0.8047699546609701
0.5475753654140099
2.2477365846057236
-0.8045831691692795
0.548062754720239
2.247689971608702
-0.8044202368616267
0.5484878456658577
2.247648928822684
-0.8042767784150748
0.5488620852871939
2.247612496621387
-0.8041494383472365
0.5491942424585874
2.2475799265089256
-0.8040355995369417
0.549491155140329
2.2475506255628757
-0.8039331889051979
0.5497582388715214
2.2475241176644047
-0.8038405418187562
0.5499998415097722
2.2475000158489977
-0.803756305392247
0.5502194961499582
2.2474780022064444
-0.8036793681954264
0.5504201049415738
2.2474578130176806
-0.8036088082765351
0.5506040749838319
2.247439227595031
-0.8035438541379181
0.550773420339262
2.247422059787053
-0.8034838550304269
0.5509298396750171
2.247406151430677
-0.8034282580577705
0.551074776096486
2.247391367246694
-0.8033765903282699
0.5512094637842533
2.247377590818929
-0.8033284448961648
0.5513349647245445
2.24736472139673
-0.8032834695818716
0.5514521979146387
2.247352671329753
-0.8032413580032234
0.551561962789914
2.2473413639932516
-0.8032018423217973
0.5516649581690736
2.247330732097388
-0.8031646873319938
0.551761797690907
2.2473207162998055
-0.8031296856104655
0.5518530224807694
2.2473112640596087
-0.8030966535095878
0.5519391116121133
2.2473023286849445
-0.8030654278278302
0.5520204907998738
2.24729386853694
-0.8030358630267855
0.5520975396660401
2.2472858463607452
-0.8030078288925862
0.552170597844629
2.2472782287205337
-0.8029812085607659
0.5522399701375258
2.2472709855200304
-0.802955896840138
0.5523059308895117
2.2472640895937817
-0.8029317987839866
0.5523687277175251
2.247257516357248
-0.8029088284668982
0.5524285847030148
2.2472512435060388
-0.8028869079333987
0.5524857051357475
2.2472452507564036
-0.802865966290808
0.5525402738811191
2.2472395196204973
-0.8028459389236664
0.5525924594301102
2.2472340332110923
-0.8028267668111004

plt.figure(figsize=(4,4))
plt.plot(x, phi(x));
plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
#plt.axis([x[0],x[-1],x[0],x[-1]])
plt.plot(xh, phih,'g', lw=0.5)
plt.plot(xh[0],phih[0],'go')
plt.show()

def phi(x):
    # f(x) = - (x^2 - 2) = 0
    #  ==>  x = x - (x^2 - 2) / 2
    return x - (x**2 - 2.)/2.

xh, phih = FixedPointIteration(0.25, phi)

plt.figure(figsize=(4,4))
plt.plot(x, phi(x));
plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go')
plt.show()

1.21875
1.47607421875
1.3866766691207886
1.425240576778826
1.4095852259304056
1.4161199713497692
1.4134220847219336
1.4145410899320852
1.4140778423789595
1.4142697702253928

Newton's Method

# Let's illustrate what Newton's method does.

def f(x):
    return 2 - x**2
def fp(x):
    return -2.*x

x = np.linspace(0,2,200)
plt.figure(figsize=(8,8))

plt.plot(x,np.zeros(x.shape),'k', x, f(x)); plt.grid('on');
a = plt.axis()
x0 = 1.; f0 = f(x0);
m = fp(x0)
f_left = f0 + m*(a[0]-x0); f_right = f0 + m*(a[1]-x0);
plt.plot(a[:2], [f_left, f_right], 'r');
plt.axis(a);
plt.plot([x0], [f0], 'ro');
x1 = x0 - f0/m
plt.plot([x1], [0], 'ro');

# Now let's try it on phi(x) = cos(x)

def phi_Newton(x):
    return x - (-np.cos(x))/np.sin(x)

xh, phih = FixedPointIteration(0.8, phi_Newton)

1.7712146006504745
1.5680690668947184
1.5707963335566546
1.5707963267948966
1.5707963267948966
1.5707963267948966
1.5707963267948966
1.5707963267948966
1.5707963267948966
1.5707963267948966

plt.figure()
plt.plot(x,x,'r', x, x+1./np.tan(x))
plt.axis([0.5, 2, 0.5, 2])
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go');

/var/folders/97/ppgk101x3214hn7bzlb85jk00000gn/T/ipykernel_54126/2142578811.py:2: RuntimeWarning: divide by zero encountered in divide
  plt.plot(x,x,'r', x, x+1./np.tan(x))