乘幂法、复化梯形以及二分法求特征值【相关例题python实现】

写在前面的话

考试结束，《数值计算方法》编程作业发出来给大家看看~

乘幂法

例题

《数值计算方法》P185页例题6.4，误差设定到0.00001时需要10次迭代，所得到的结果与书上的结果接近。当把误差限定得更小则迭代的次数更多。

代码

import numpy as np
import matplotlib.pyplot as plt

k = 0
e = 1
temp_max = []
z0 = np.mat([1, 1, 1])
z0 = z0.T
A = np.mat([[2,0,1],
            [1,-1,2],
            [0,1,5]], dtype=float)
while e > 0.00001:
    y = A * z0
    y1 = 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 + 1
    if k > 1:
        e = abs(y[x] - temp_max[-2])
λ =  y[x]
print('最大特征值为：', λ)
print('相应特征向量为：\n', (A*z0))

结果

二分法求对称三对角矩阵特征值

代码

import numpy as np
import matplotlib.pyplot as plt

data = [
    [2, 2, 0, 0],
    [1, 1, 2, 0],
    [0, 2, 1, 2],
    [0, 0, 2, 1]
]


def P0(x, data):
    return 1


def P1(x, data):
    return data[0][0] - x


def 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 1
    elif a < 0 and b < 0:
        return 1
    else:
        return 0


def P(x, PX, data):
    ret = []

    count = 0
    for 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 += 1
        if ret[i - 1] == 0:
            count += ifTong(ret[i], ret[i - 2])
        else:
            count += ifTong(ret[i], ret[i - 1])

    ret.append(count)
    return ret


if __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 math
import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    y= np.exp(x)*np.cos(x)
    return y
def result(f,a,b,n):
    ti=0.0
    h=(b-a)/n
    ti=f(a)+f(b)
    for k in range(1,int(n)):
        xk=a+k*h
        ti = ti + 2 * f(xk)
    result = ti*h/2
    print("复化梯形公式计算结果为：", result)

if __name__ == '__main__':
    n = 2
    k = 10
    for i in range(1,k):
        print("n= ", n)
        result(f1,0,math.pi,n)
        n = n*2

结果