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

I share my academic life!

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

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

 

 

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 )

 

 

Sympson 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 )

 

Tags: , , , , , , ,