محاسبات عددی با پایتون

خلاصه: کدهای پایتون درس محاسبات عددی: دوبخشی، نیوتن، نقطه ثابت، درون‌یابی لاگرانژ، و انتگرال‌گیری ذوزنقه‌ای و سیمپسون.

کدهایی که برای این درس تهیه کرده بودم رو در ادامه میزارم.

Bisection

import math
import decimal
from prettytable import PrettyTable
 
def get_length(interval):
    return decimal.Decimal((interval[1] - interval[0]))
 
 
def split_interval(interval):
    left = decimal.Decimal(interval[0])
    right = decimal.Decimal(interval[1])
    interval_length = decimal.Decimal(get_length(interval))
    if not (fun(left) * fun(right)) < 0:
        print ("bAze EntekhAb Shode Eshtbh Ast.")
        return False
 
    split_point = decimal.Decimal(( left + (interval_length / 2) ))
 
    if (fun(left) * fun(split_point)) < 0:
        return [left,split_point]
    elif (fun(split_point) * fun(right)) < 0:
        return [split_point, right]
    else:
        print ('bAze Sahih nist.')
        return False
 
 
def run_bisection(interval, target_length):
    interval_length = decimal.Decimal(get_length(interval))
    while interval_length > target_length:
        interval = split_interval(interval)
        if not interval:
            print ("Function returned False! (spliting the interval)")
            return False
        interval_length = get_length(interval)
        t.add_row([interval[0], interval[1],interval_length])
    print ("Bisection Method\n\n|_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_--| \n\n  Result: [%f : %f]" % (interval[0], interval[1]))
    print (t)
 

 

def fun(x):
    return pow(x,2)- 2

interval = [0, 2] 
length = decimal.Decimal('0.00000000001') 

t = PrettyTable(['a', 'b','interval size'])

run_bisection(interval, length)
k=input(":")

Newton

import time
import numpy as np
from prettytable import PrettyTable
tab = PrettyTable(['m=(f(x+h) - f(x))/h', 'b=f(x)-m','-b/m'])
def rootNR(f, start_x):
    step = 1e-10 
    x = start_x
    eps = 1e-15 
    
    while True:
        m = (f(x+step) - f(x))/step
        b = f(x)-m*x
        x_0 = -b/m
        tab.add_row([m, b,x_0])
        if (np.abs(x-x_0) <= eps):
            break
        else:
            x = x_0
            
    return x, np.abs(x-x_0)
 
if __name__ == "__main__":

    
  

    t = time.time()   
    
    
    def funct(x):
        return (x*x)-2
    
    root, err = rootNR(funct, 1)
    print 'Root:', root
    print 'Elapsed time: %s ms' % ((time.time() - t)*1000)
    print 'Error:', err
    print tab
    k=input("enterkey to close ")

Fixed Point

import math
def fixedpoint( f, x,  tol, maxiter ):
    x1 = x + 8.0 * tol
    x0 = x + 4.0 * tol
    a= 2 * tol / ( f( x - tol ) - f( x + tol ) )
    k = 0

    while k <= maxiter and abs( x - x0 ) >= tol:
        x2, x1, x0 = ( x1, x0, x )
        x = x + a * f( x )
        print ("%2d %18.11e %18.11e") % ( k, x, abs( x - x0 ) )
        k = k + 1

    if k > maxiter:
        print ("Error: exceeded %d iterations") % maxiter

    rate = math.log( abs((x - x0) / (x0 - x1)) ) / \
           math.log( abs((x0 - x1) / (x1 - x2)) )

    return ( x, rate )

#------------------------------------------------------------------------------

if __name__ == "__main__":

     
    maxiter = 100

    
    tol = 1e-10
    
    def  f(x):
        return x*x-2

        
    def df(x):
        return x-2
     

        

    a, b = ( 0.0, 2.0 )
    first_guess = 1
    
   
    dfr=0
    print ("\n-------- Fixed-Point Method -------------------------\n")
    x,  rate = fixedpoint( f, first_guess,  tol, maxiter )
    print ("root = "), x
    print ("Estimated Convergence Rate = %5.2f") % rate


    k=input(":")

Lagrange Polynomial

import numpy as np import matplotlib.pyplot as plt import sys def main(): if len(sys.argv) == 1 or "-h" in sys.argv or "--help" in sys.argv: print ("python lagrange.py .. ") print ("Example:") print ("python lagrange.py 0.1 2.4 4.5 3.2") exit() points = [] for i in xrange(len(sys.argv)): if i != 0: points.append((int(sys.argv[i].split(".")[0]),int(sys.argv[i].split(".")[1]))) #points =[(0,0),(25,30),(50,10), (57,0)] P = lagrange(points) nr = 2 print (" + str(points[nr][0]) + ", " + str(points[nr][1]) +") P(" + str(points[nr][0]) +")= " +str(P(points[nr][0])) plot(P, points) def plot(f, points): x = range(-10, 100) y = map(f, x) print y plt.plot( x, y, linewidth=2.0) x_list = [] y_list = [] for x_p, y_p in points: x_list.append(x_p) y_list.append(y_p) print x_list print y_list plt.plot(x_list, y_list, 'ro') plt.show() def lagrange(points): def P(x): total = 0 n = len(points) for i in xrange(n): xi, yi = points[i] def g(i, n): tot_mul = 1 for j in xrange(n): if i == j: continue xj, yj = points[j] tot_mul *= (x - xj) / float(xi - xj) return tot_mul total += yi * g(i, n) return total return P if __name__ == "__main__": main()

Trapezoidal Integral

import time def trapezoidal(f, a, b, n): h = float(b - a) / n s = 0.0 s += f(a)/2.0 for i in range(1, n): s += f(a + i*h) s += f(b)/2.0 return s * h a=input("Enter a:" ) a1=eval(a) b=input("Enter b:" ) b1=eval(b) n=input("Enter n: ") n1=eval(n) t = time.time() trap=trapezoidal(lambda x:x*x, a1,b1, n1) #mitunid tabe ro ba'd az ":" avaz konid masalan lambda x: 3*x print(trap) print ("Elapsed time in ms:") print( (time.time() - t)*1000 )

Simpson Integral

import time def simpson(f, a, b, n): h=(b-a)/n k=0.0 x=a + h for i in range(1,n//2 + 1): k += 4*f(x) x += 2*h x = a + 2*h for i in range(1,n//2): k += 2*f(x) x += 2*h return (h/3)*(f(a)+f(b)+k) t = time.time() def function(x): return x*x l =simpson(function, 3.0, 6.0, 100) print (l) print ("Elapsed time in ms:") print( (time.time() - t)*1000 )