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


کدهایی که برای این درس تهیه کرده بودم رو توی ادامه مطلب میزارم
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 )