目录
写在前面的话乘幂法例题代码结果二分法求对称三对角矩阵特征值代码结果复化梯形例题代码结果写在前面的话
考试结束，《数值计算方法》编程作业发出来给大家看看~
乘幂法
例题
《数值计算方法》P185页例题6.4，误差设定到0.00001时需要10次迭代，所得到的结果与书上的结果接近。当把误差限定得更小则迭代的次数更多。
代码
import numpy as npimport matplotlib.pyplot as pltk = 0e = 1temp_max = []z0 = np.mat([1, 1, 1])z0 = z0.TA = np.mat([[2,0,1],[1,-1,2],[0,1,5]], dtype=float)while e > 0.00001:y = A * z0y1 = y.copy()y1 = abs(y1)x = y1.argmax()z0 = y/y[x]print('迭代次数：',k,'y：',y.T,'m：',y[x],'z：',z0.T)print(' ')temp_max.append(y[x])k = k + 1if k > 1:e = abs(y[x] - temp_max[-2])λ = y[x]print('最大特征值为：', λ)print('相应特征向量为：\n', (A*z0))
结果
二分法求对称三对角矩阵特征值
代码
import numpy as npimport matplotlib.pyplot as pltdata = [[2, 2, 0, 0],[1, 1, 2, 0],[0, 2, 1, 2],[0, 0, 2, 1]]def P0(x, data):return 1def P1(x, data):return data[0][0] - xdef P2(x, data):return (data[1][1] - x) * P1(x, data) - data[1][0] ** 2 * P0(x, data)def P3(x, data):return (data[2][2] - x) * P2(x, data) - data[2][1] ** 2 * P1(x, data)def P4(x, data):return (data[3][3] - x) * P3(x, data) - data[3][2] ** 2 * P2(x, data)def ifTong(a, b):if a > 0 and b > 0:return 1elif a < 0 and b < 0:return 1else:return 0def P(x, PX, data):ret = []count = 0for i in range(PX):if i == 0:ret.append(1)elif i == 1:ret.append(data[i - 1][i - 1] - x)else:ret.append((data[i - 1][i - 1] - x) * ret[-1] - data[i - 1][i - 2] ** 2 * ret[-2])for i in range(1, len(ret)):if ret[i] == 0:a = ret[i - 1]count += 1if ret[i - 1] == 0:count += ifTong(ret[i], ret[i - 2])else:count += ifTong(ret[i], ret[i - 1])ret.append(count)return retif __name__ == '__main__':# x1,x2,x3,x4,x5=P0(x,data),P1(x,data),P2(x,data),P3(x,data),P4(x,data)# print(x1,x2,x3,x4,x5)x = np.arange(-3, 5.5, 0.000025)ret = []for i in x:ret.append(P(i, 5, data))ret = np.transpose(ret)ret = np.insert(ret, 0, x, axis=0)print(ret)ans = []for i in range(len(ret[0]) - 1):ans.append(ret[-1][i] - ret[-1][i + 1])# print(ans)for i in range(len(ans)):if ans[i] != 0:print('有一个值在', ret[0][i], '和', ret[0][i + 1], '之间')pass
结果
复化梯形
例题
《数值计算方法》P103页例题4.1中表4.3的结果和代码运行结果一致。
代码
import mathimport numpy as npimport matplotlib.pyplot as pltdef f1(x):y= np.exp(x)*np.cos(x)return ydef result(f,a,b,n):ti=0.0h=(b-a)/nti=f(a)+f(b)for k in range(1,int(n)):xk=a+k*hti = ti + 2 * f(xk)result = ti*h/2print("复化梯形公式计算结果为：", result)if __name__ == '__main__':n = 2k = 10for i in range(1,k):print("n= ", n)result(f1,0,math.pi,n)n = n*2
