Pages

Bracketing methods


Chapter 5
In [180]:
import matplotlib.pyplot as plt
import math
import numpy as np
%matplotlib inline
##%matplotlib qt ------this will pop up a new window with the graph
from IPython.core.display import Image
In [181]:
Image(filename='example_5.2.PNG')
Out[181]:
In [182]:
x=np.linspace(0,5,1000)
def f(x):
    y=np.sin(10*x)+np.cos(3*x)
    return y
In [183]:
import scipy.optimize as opt
In [184]:
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
In [185]:
fig, ax = plt.subplots(figsize=[16, 8])
ax.plot(x,f(x))
ax.plot(x,np.zeros(len(x)))
ax.set_xlim(0,6)
ax.set_ylim(-3,3)
ax.grid()
ax.set_xticks(np.arange(0,6,.2))
ax.set_yticks(np.arange(-3,3,0.5))
axins = zoomed_inset_axes(ax, 12, loc=1)  # zoom = 6
axins.plot(x,f(x))
axins.plot(x,np.zeros(len(x)))
axins.grid()
axins.set_xlim(4.2,4.3)
axins.set_ylim(-0.05,.05)
mark_inset(ax, axins, loc1=2, loc2=3, fc="none", ec="0.5")
Out[185]:
(<mpl_toolkits.axes_grid1.inset_locator.BboxPatch at 0xa01147d518>,
 <mpl_toolkits.axes_grid1.inset_locator.BboxConnector at 0xa01147df60>,
 <mpl_toolkits.axes_grid1.inset_locator.BboxConnector at 0xa011482898>)
In [186]:
t=opt.bisect(f,4.22,4.24,full_output=True)
t
Out[186]:
(4.229067033678874,       converged: True
            flag: 'converged'
  function_calls: 37
      iterations: 35
            root: 4.229067033678874)
In [187]:
opt.fsolve(f,4.22,full_output=True)
Out[187]:
(array([ 4.22906703]),
 {'fjac': array([[-1.]]),
  'fvec': array([ -4.44089210e-16]),
  'nfev': 8,
  'qtf': array([ -6.68761713e-11]),
  'r': array([ 1.56698199])},
 1,
 'The solution converged.')
In [188]:
import time

start_time = time.clock()
t=opt.fsolve(f,4.22)
print (time.clock() - start_time, "seconds")
print(t)
tt=opt.bisect(f,4.22,4.24)
print (time.clock() - start_time, "seconds")
print (tt)

0.002644758715177886 seconds
[ 4.22906703]
0.027985479347989894 seconds
4.229067033678874
In [189]:
start_time=time.clock()
ttt=opt.newton(f,4.22)
print(time.clock() - start_time,"seconds")
print(ttt)

0.01125393039546907 seconds
4.22906703368
In [190]:
Image(filename='example_5.6.PNG')
Out[190]:
In [191]:
def fun(x):
    return x**10-1
xx=np.linspace(0,1.5,1000)
plt.plot(xx,fun(xx))
plt.grid()

Bisection method

Algorithm

In [192]:
Image(filename='bisection.PNG')
Out[192]:
In [193]:
def new_bisect(xl,xu,error_tolerance,imax,xr,ea):
    iter=0
    fl=fun(xl)
    while(1):
        xrold=xr
        xr=(xl+xu)/2
        fr=fun(xr)
        iter=iter+1
        if xr!=0 :
            ea=abs((xr-xrold)/xr)*100
            #print("ea+++++===",ea)
        test=fl*fr
        ##print("f(xl)=====",fun(xl))
        if test < 0 :
            xu=xr
        elif test > 0 :
            xl=xr
            fl=fr
        else:
            ea=0
        if (ea < error_tolerance or iter >= imax):
            print("xr======",xr)
            print("iter=======",iter)
            print("ea================",ea)
            break
        bisect=xr
In [194]:
new_bisect(0,1.3,.01,100,0,100)

xr====== 0.9999938964843751
iter======= 14
ea================ 0.007934618741565026

False position method(Regula falsi method)

Algorithm

In [195]:
Image(filename='falsepos.PNG')
Out[195]:
In [196]:
  def false_position(xl,xu,error_tolerance,imax,xr,ea):
    iter=0
    fl=fun(xl)
    fu=fun(xu)
    while(1):
        xrold=xr
        xr=xu-fu*(xl-xu)/(fl-fu)
        fr=fun(xr)
        iter=iter+1
        if xr!=0 :
            ea=abs((xr-xrold)/xr)*100
        
        test=fl*fr
        
        if test < 0 :
            xu=xr
            fu=fun(xu)
            iu=0
            il=il+1
                
        elif test > 0 :
            xl=xr
            fl=fun(xl)
                  
        else:
            ea=0
        if (ea < error_tolerance or iter >= imax):
            print("xr======",xr)
            print("iter======",iter)
            print("ea=======",ea)
            break
        false_position=xr
In [197]:
false_position(0,1.3,.01,100,0,100)

xr====== 0.9996886513360608
iter====== 39
ea======= 0.009537956117341062
In [ ]:
In [ ]:

Modified false position method

Algorithm

In [198]:
Image(filename='modfalsepos.PNG')
Out[198]:
In [199]:
  def modified_false_position(xl,xu,error_tolerance,imax,xr,ea):
    iter=0
    il=0
    iu=0
    fl=fun(xl)
    fu=fun(xu)
    while(1):
        xrold=xr
        xr=xu-fu*(xl-xu)/(fl-fu)
        fr=fun(xr)
        iter=iter+1
        if xr!=0 :
            ea=abs((xr-xrold)/xr)*100
        
        test=fl*fr
        
        if test < 0 :
            xu=xr
            fu=fun(xu)
            iu=0
            il=il+1
            
            if il >= 2:
                fl=fl/2
                
        elif test > 0 :
            xl=xr
            fl=fun(xl)
            il=0
            iu=iu+1
            
            if iu >= 2 :
                fu=fu/2
                
                
        else:
            ea=0
        if (ea < error_tolerance or iter >= imax):
            print("xr======",xr)
            print("iter======",iter)
            print("ea=======",ea)
            break
        modified_false_position=xr
In [ ]:
In [200]:
modified_false_position(0,1.3,.01,100,0,100)

xr====== 1.0000057046332558
iter====== 12
ea======= 0.0011642999121982691
In [ ]:

২টি মন্তব্য:

এতে সদস্যতা: মন্তব্যসমূহ (Atom)