数值分析python代码(共12页).doc

《数值分析python代码(共12页).doc》由会员分享，可在线阅读，更多相关《数值分析python代码(共12页).doc（12页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上from math import logfrom math import sin,cosimport random#Q1def bifection(f,a,b,error,n):f stands for the function we want to approximate to, while a and b are theleft and right number of the interval where the f effects in, and n means for thetime we have tried to make a approximationb

2、ifection(f,0,1,10*-6,20)FA = f(a)for i in range(n):p = a + (b-a)/2FP = f(p)if FP = 0 or (b - a)/2 0:a = pelse:b = pdef iteration(f,p0,n,error):iteration(f,0.5,22,10*(-6)p = 0for i in range(n):p = (2.71828)*(-p0)if abs(p - p0) error:return pelse:p0 = pdef fastitertion(f,p0,n,error):fastitertion(f,0.5

3、,20,10*-6)p = 0q0 = p0for i in range(n):p = (2.71828)*(-p0)q = p + g(p)/(1 - g(p)*(p - q0)if abs(q - q0) error:return qelse:p0 = pq0 = qdef newton(f,g,p0,n,error):newton(f,0.5,3,10*-6)for i in range(n):p = p0 - f(p0)/g(p0)if abs(p - p0) error:return pelse:p0 = pdef secant(f,p0,p1,n,error):secont(f,0

4、.5,1,6,10*-6)for i in range(n):p = p0 - f(p0)*(p1 - p0)/(f(p1) - f(p0)if abs(p - p0) error:return pelse:p0 = pdef fastsecant(f,p0,p1,n,error):secont(f,0.5,1,4,10*-6)for i in range(n):p = p0 - f(p0)*(p1 - p0)/(f(p1) - f(p0)if abs(p - p0) group = 2,3,-1,4,1,-5,3,0,4,-1,-5,0,5,-1,1,-4,1,8,-1,3,4,5,-3,4

5、,4,-4,3,7,5,1 gausseli(group)1 0 0 0 0 -10 1 0 0 0 20 0 1 0 0 40 0 0 1 0 -20 0 0 0 1 3for i in range(len(group):first = groupiifor j in range(i,len(group0):groupij = groupij / firstfor temp in range(len(group):if temp != i:second = grouptempifor j in range(i,len(group0):grouptempj -= groupij * secon

6、dreturn groupdef gausseli(group):1.0, 1.5, -0.5, 2.0, 0.5, -2.5-0.0, 1.0, -1.2, 1.5, 1.4, -1.6-0.0, -0.0, 1.0, 0.1, -1.5, -0.9-0.0, -0.0, -0.0, 1.0, -5.0, -17 0.0, 0.0, 0.0, 0.0, 1.0, 3.0for i in range(len(group):first = groupiifor j in range(len(group0):groupij = groupij / firstfor temp in range(i

7、+ 1,len(group):second = grouptempifor j in range(i,len(group0):grouptempj -= groupij * secondreturn groupdef jordan(group):2, 3, -1, 4, 1, -50.0, -4.5, 5.5, -7.0, -6.5, 7.50.0, 0.0, -6.9, -0.8, 10.8, 6.30.0, 0.0, 0.0, -1.0, 5.1, 17.30.0, 0.0, 0.0, 0.0, 83.0 249.0for i in range(len(group):first = gro

8、upiifor temp in range(i + 1,len(group):second = grouptempifor j in range(i,len(group0):grouptempj -= groupij / first * secondreturn group#fittingdef fitting(lst,f,g):f and g is the function which translate the original function into lineari = 0total = 0lst0 = g(item) for item in lst0lst1 = f(item) f

9、or item in lst1x_bar = sum(lst0)/len(lst0)y_bar = sum(lst1)/len(lst1)while i lst1 = 1,2,3,4,5,6,7,8,15.3,20.5,27.4,36.6,49.1,65.6,87.8,117.6a = fitting(lst1, log, lambda x: x)a = (0., 2.8185)#a = e(a)#y = 11.4375e(0.291x)#Q7#1/y = b/t + alst2 = 1,2,3,4,5,6,7,8,4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86

10、b = fitting(lst2,lambda y: 1/y, lambda x: 1/x)b = (0., 0.)#y = t/(0.075t + 0.171)#Q8#y = a + bx3vender = -5,-4,-3,-2,-1,0,1,2,3,4,5beta = -12.9,-5.95,-1.76,0.42,1.20,1.34,1.43,2.25,4.38,8.61,15.6lst = vender,betab = fitting(lst,lambda y: y, lambda x: x*3)(0.11, 1.33)#y = 1.33 + 0.11x3def matrix_tran

11、s(matrix):a,b a,cc,d - b,dreturn rowi for row in matrix for i in range(len(matrix0)def matrix_inver(matrix):for temp in range(len(matrix):matrixtemp.extend(0 for x in range(len(matrix)for i in range(len(matrix):matrixilen(matrix) + i = 1inverse = gaussjordaneli(matrix)for i in range(len(inverse):inv

12、ersei = inverseilen(inverse):return inversedef matrix_mul(A, B):result = 0 * len(B0) for i in range(len(A)for i in range(len(A):for j in range(len(B0):for k in range(len(B):resultij += Aik * Bkjreturn resultdef polyfitting(vender,beta):vender = 1,0,125,1,0,64,1,0,27,1,0,8,1,0,1,1,0,0,1,0,1,1,0,8,1,0

13、,27,1,0,64,1,0,125beta = -12.9,-5.95,-1.76,0.42,1.20,1.34,1.43,2.25,4.38,8.61,15.6new_beta = betai for i in range(len(beta)temp = matrix_mul(matrix_inver(matrix_mul(matrix_trans(vender),vender),matrix_trans(vender)result, i = , 0while i lst =1.0,1.5,2.0,3.5,5.5,6.0,1.0,3.375,8.0,42.875,166.375,216.0

14、lagrangeitp(lst,2.5)15.625lagrangeitp(lst,4.5)91.125return sum(la_multool(lst0,i,x) * lst1i) for i in range(len(lst0)def la_multool(lst,i,x):k = 0total = 1while k j:return 0if i = j :return lst1ielse:return (newton_multool(lst,i+1,j) - newton_multool(lst,i,j - 1)/(lst0j - lst0i)#Q11lst = 0.4,0.5,0.6

15、,0.7,0.8,0.9,-0.,-0.,-0.,-0.,-0.,-0.newtonitp(lst,0.78)- 0.log(0.78)- 0.84996error = 0.0016#Q12lst = 0.785, 0.873, 0.959, 1.047,0.7071,0.7660,0.8192,0.8660newtonitp(lst,0.907)0.7071sin(0.907)0.7876error = 0.1#Q13lst = 0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0

16、.9463q1 = lagrangeitp(lst00:3,lst10:3,0.47)0.sin(0.47)0.error = 0.0008q2 = lagrangeitp(lst04:7,lst14:7,1.6)0.sin(1.6)0.error = 10(-5)class integral(object):def _init_(self,function,a,b):self.f = functionself.a = aself.b = bdef simpson(self,n):h = (self.b - self.a)return h * (self.f(self.a) + 4 * sel

17、f.f(1/2 *(self.a + self.b) + self.f(self.b) / 6def fixed_simpson(self,n):h = (self.b - self.a)/ nresult = self.f(self.a) + self.f(self.b)for k in range(n):result += 4 * self.f(self.a + (k - 1/2)*h) + 2 * self.f(self.a + k * h)return h * result / 6def fixed_cotes(self,n):h = (self.b - self.a)/ nresul

18、t = (self.f(self.a) + self.f(self.b) * 7for k in range(n):result += 32 * self.f(self.a + (k - 1/4)* h) + 12 * self.f(self.a + (k - 1/2) * h) +32 * self.f(self.a + (k - 3/4)* h) + 14 * self.f(self.a + k * h)return h * result / 90def romberg(self,n):R = 0 for x in range(n) for x in range(n)h = (self.b

19、 - self.a)/ nR00 = (self.f(self.a) + self.f(self.b) * h / 2for i in range(n - 1):for j in range(i):Rij = self.T(i,j,n)return Rdef T(self,i,k,n):if i = 0 and k = 0:return (self.f(self.a) + self.f(self.b) * (self.b - self.a)/ (n * 2)if i = 0 and k != 0:return 1/2 * self.T(0,k - 1,n) + (self.b - self.a

20、) / (2 * k) * sum(self.f(self.a + (2*i + 1)*(self.b - self.a)/(2 * k) for i in range(2 * (k - 1) - 1)else:return (4 * i * self.T(i - 1,k+1,n) - self.T(i - 1,k,n)/(4*i - 1)#Q14 E1 = integral(lambda x:(2.71828)*(-x),0,1) E1.simpson(90)0.10494 original = 0.33684#Q15s1 = integral(lambda x: sin(x)/x ,0 +

21、 1/64,5) s1.fixed_simpson(64)1. 70267s1 = integral(lambda x: sin(x)/x ,0 + 1/128,5) s1.fixed_simpson(128)1. 30339s1 = integral(lambda x: sin(x)/x ,0 + 1/256,5) s1.fixed_simpson(256)1. 93797#Q16y1 = integral(lambda x : 1/x , 1, 3)y1.romberg(4)class diff_equation(object):def _init_(self,function,initi

22、al):function = lambda x : lambda y :self.f = functionself.initial = initialdef Euler(self,a,b,h):N = (b - a)/htemp = aresult = self.initialfor i in range(int(N):print(result)result += h * self.f(temp + i * h, result)temp = temp + hreturn resultdef improved_euler(self,a,b,h):N = (b - a)/htemp = aresu

23、lt = self.initialfor i in range(int(N):print(result)temp1 = resultresult += h / 2 * (self.f(temp + i * h, temp1 + h * self.f(temp, temp1) + self.f(temp,temp1)temp = temp + hreturn resultdef Runge_Kutta4(self,a,b,h):N = (b - a)/htemp = aresult = self.initialK1 = self.f(a,result)K2 = self.f(a + 1/2 *

24、h, result + 1 / 2 * K1 * h)K3 = self.f(a + 1/2 * h, result + 1 / 2 * K2 * h)K4 = self.f(a + h, result + K3 * h)for i in range(int(N):print(result)result += 1/6 * h *(K1 + K2 + K3 + K4)temp = temp + hK1 = self.f(temp,result)K2 = self.f(temp + 1/2 * h, result + 1 / 2 * K1 * h)K3 = self.f(temp + 1/2 *

25、h, result + 1 / 2 * K2 * h)K4 = self.f(temp + h, result + K3 * h)return result#Q170.00.20.520.8080.96161. 01. 01. 01. 01. 00.0,0.,0.33684,0.77936,0.59726,0.92415,0.3814,0.43657,0.09424,0.54715#Q18 f = lambda x,y : x + y k = diff_equation(f,1) k.Runge_Kutta4(0,1,0.2)11. 333321. 444431. 837041. 982562

26、. 11494 k.improved_euler(0.1,0.2)11. 221. 55242.3.4.1.0, 1.84079, 1.86987, 2.75187, 2.9733#Q19def diff_sys(b,interval,f0,g0):f = lambda x,y : x * (1 - x * 2 - y * 2) + y * (x*2 - y*2 - b)g = lambda x,y : y * (1 - x * 2 - y * 2) - x * (x*2 - y*2 - b)N = (interval1 - interval0)/0.05temp = interval0res

27、ult = f0,g0for i in range(int(N):print(result)result0 += 0.05 * f(result0,result1)result1 += 0.05 * g(result0,result1)temp = temp + 0.05return result#Q20def partial(h,k):pU/pT = u(n+1,k) - u(n,k)/dtp2U/px2 = u(n,k+1) - 2 * u(n,k) +u(n,k - 1)/dx2u|t = 0 u(0,i) = 4(i*dx)*(1 - i*dx)u|x = 0 = u| x = 1 = 0#Q22 monte_carlo(100,lambda x:2.71828 * x)1. monte_carlo(10000,lambda x:2.71828 * x)1. 00001def monte_carlo(number,f,a,b):count = 0for i in range(1,number + 1):x = random.uniform(0,a)y = random.uniform(0,b)if y f(x):count += 1return count / number * (a * b)专心-专注-专业