工程数学算法python实现

antRain

于 2022-01-10 15:36:05 发布

阅读量1.7k

点赞数 1

分类专栏： python 文章标签： python 算法

本文链接：https://blog.csdn.net/qq_41146650/article/details/122409534

版权

python 专栏收录该内容

15 篇文章 1 订阅

订阅专栏

工程数学算法python实现

非线性方程求根

二分法

def dichotomy(func,a,b,is_print=False,e=1e-8,max_iter=100):
    '''
    func: 需要进行二分法计算零点的函数
    [a,b] 进行二分的区间
    e: 精度，默认1e-8
    max_iter: 最大迭代次数，默认为100
    is_print: 是否打印中间过程，默认不打印（k,ak,bk,xk,f(xk)）  
    '''
    a0=a
    b0=b
    fa = func(a)
    fb = func(b)
    if fa*fb > 0:
        raise Exception('不符合要求，fa*fb>0')
    k = 0
    if is_print:
        print('k\t ak\t bk\t xk\t f(xk)\t')
    while k<max_iter:
        xk = a + (b-a)/2
        fx = func(xk)
        if is_print:
            print('%d\t %f\t %f\t %f\t %f\t' % (k,a,b,xk,fx))
        k+=1
        if abs(xk) < e or (b-a)/2<e:
            if is_print:
                print('| x - '+ chr(945) +'| = ' + str(((b0-a0)/2**k)))
            return xk
        if fa*fx>0:
            a = xk
            fa = fx
        else :
            b = xk
     
    raise Exception('超过最大的迭代次数')

def dichotomy_test_func(x):
    return x**3-x-1
    
print(dichotomy(dichotomy_test_func, 1, 1.5,True,0.005))

简单迭代

def simple_iteration(func,x,is_print=False,e=1e-8,max_iter=100):
    '''
    func: 利用简单迭代法的迭代函数
    x：初始值
    e: 精度，默认1e-8
    max_iter: 最大迭代次数，默认为100
    is_print: 是否打印中间过程，默认不打印（k,ak,bk,xk,f(xk)）  
    '''
    k = 0
    print('k\t xk\t')
    print('%d\t %f\t' % (k,x))
    while k<max_iter:
        k+=1
        xk=func(x)
        if is_print:
            print('%d\t %f\t' % (k,xk))
        if abs(xk-x)<e:
            return xk
        x = xk
        

def iteration_test_func1(x):
    '''
    原函数：f(x) = x^3 - x - 1
    '''
    return (x+1)**(1/3)

def iteration_test_func2(x):
    '''
    原函数：f(x) = x^2 - 3
    '''
    return x - 0.25*(x*x-3)

def iteration_test_func3(x):
    '''
    原函数：f(x) = x^2 - 3
    '''
    return 0.5*(x+3/x)

print(simple_iteration(iteration_test_func1, 1.5,True,1e-5))
print(simple_iteration(iteration_test_func2, 2,True,1e-5))
print(simple_iteration(iteration_test_func3, 2,True,1e-5))

牛顿迭代

def newton_iteration(func,dfunc,x,is_print=False,e=1e-8,max_iter=100):
    '''
    func: 利用牛顿迭代法的原函数 fx
    dfunc: fx的一阶导数
    x：初始值
    e: 精度，默认1e-8
    max_iter: 最大迭代次数，默认为100
    is_print: 是否打印中间过程，默认不打印（k,ak,bk,xk,f(xk)）  
    '''
    k = 0
    print('k\t xk\t')
    print('%d\t %f\t' % (k,x))
    while k<max_iter:
        k+=1
        xk=x - func(x)/dfunc(x)
        if is_print:
            print('%d\t %f\t' % (k,xk))
        if abs(xk-x)<e:
            return xk
        x = xk

def newton_test_func1(x):
    return x*math.exp(x)-1

def newton_test_dfunc1(x):
    return (x+1)*math.exp(x)

print(newton_iteration(newton_test_func1, newton_test_dfunc1, 0.5,True,1e-5))

线性方程组的数值解法

三角形线性代数方程组的直接解法

def back_substitution(A,b):
    '''
    回代法求解方程组Ax=b的解, 其中A为下三角行列式
    '''
    n = len(b)
    x = [0 for i in range(n)]
    for i in range(n-1,-1,-1):
        sum = 0
        for j in range(n-1,i,-1):
            sum+=A[i][j]*x[j]
        x[i] = (b[i]-sum)/A[i][i]
    return x

# A = [[1,1,1],[0,1,3],[0,0,2]]
# b = [1,-2,4]
# print(back_substitution(A, b))

def front_substitution(A,b):
    '''
    前代法求解方程组Ax=b的解, 其中A为上三角行列式
    '''
    x = []
    n = len(b)
    for i in range(n):
        sum = 0
        for j in range(i):
            sum+=A[i][j]*x[j]
        x.append((b[i]-sum)/A[i][i])
    return x

# A = [[2,0,0],[3,1,0],[1,1,1]]
# b = [4,-2,1]
# print(front_substitution(A, b))

Gauss消元法

def gauss(A,b,is_print=False):
    '''
    高斯消元法求解方程组Ax=b的解
    '''
    n = len(b)
    ga = []
    gb = [i for i in b]
    # 深度复制
    for i in range(n):
        tmp = []
        for j in range(n):
            tmp.append(A[i][j])
        ga.append(tmp)
    for k in range(n-1):
        for i in range(k+1,n):
            l = ga[i][k] / ga[k][k]
            ga[i][k] = 0
            for j in range(k+1,n):
                ga[i][j] -= l*ga[k][j]
            gb[i] -= l * gb[k]
        if is_print:
            print('k= '+str(k)+':')
            print(ga)
            print(gb)
    return back_substitution(ga, gb)

# A = [[1,1,1],[1,2,4],[1,3,9]]
# b = [1,-1,1]
# print(gauss(A, b,True))

列选主元Gauss消元法

def column_pivot_element_gauss(A,b,is_print=False):
    '''
    利用列选主元高斯消元法求解方程组Ax=b的解
    '''
    n = len(b)
    # 深度复制
    ga = [ [ j for j in i] for i in A]
    gb = [i for i in b]
    for k in range(n-1):
        # 先进行列选主元
        max = ga[k][k]
        index = k
        for i in range(k+1,n):
            if ga[i][k] > max:
                max,index = ga[i][k],i
        for i in range(k,n):
            ga[k][i],ga[index][i] = ga[index][i],ga[k][i]
        gb[k],gb[index] = gb[index],gb[k]
        if is_print:
            if index == k:
                print('----------')
                print(ga)
                print(gb)
            print('k= '+str(k)+':')
            # sys.stdout.flush() # import sys 用来强制刷新缓存，输出
            # todo: 两次输出一样
        for i in range(k+1,n):
            l = ga[i][k] / ga[k][k]
            ga[i][k] = 0
            for j in range(k+1,n):
                ga[i][j] -= l*ga[k][j]
            gb[i] -= l * gb[k]
        if is_print:
            print(ga)
            print(gb)
    return back_substitution(ga, gb)

# A = [[1,-2,3],[-1,1,2],[2,4,1]]
# b = [-20,-8,9]
# print(column_pivot_element_gauss(A, b,True))

LU分解

def Doolittle(A,b,is_print=False):
    '''
    杜立特尔(Doolittle)分解法
    '''
    n = len(b)
    # 初始化
    L = [ [ 0 for j in range(n)] for i in range(n)]
    U = [ [ 0 for j in range(n)] for i in range(n)]
    # 求解U的第一行，L的第一列
    for i in range(n):
        U[0][i] = A[0][i]
        L[i][0] = A[i][0]/U[0][i]
        L[i][i] = 1
    for k in range(1,n):
        for j in range(k,n):
            sum = 0
            for q in range(k):
                sum += L[k][q]*U[q][j]
            U[k][j] = A[k][j]-sum
        for i in range(k+1,n):
            sum = 0
            for q in range(k):
                sum += L[i][q]*U[q][k]
            L[i][k] = (A[i][k]-sum)/U[k][k]
    y = front_substitution(L, b)
    x = back_substitution(U, y)
    if is_print:
        print("L =",L)
        print("U =",U)
        print("y =",y)
    return x

# A = [[1,1,1,1],[1,2,1,3],[4,3,2,1],[6,11,3,7]]
# b = [6,12,14,45]
# print(Doolittle(A,b,True))

迭代法

def jacobi(A,b,x0,is_print=False,e=1e-8,N=100,norm=None):
    '''
    利用Jacobi迭代法求解方程Ax=b的解
    x0为初始值
    e：代表迭代精度
    N: 最大的迭代次数
    norm: 代表向量范数，默认为2范数
    '''
    # 2范数
    def two_norm(x1,x2):
        n = len(x1)
        sum = 0
        for i in range(n):
            d = x1[i]-x2[i]
            sum += d*d
        return math.sqrt(sum)
    if not norm:
        norm = two_norm
    # 初始化
    n = len(b)
    k = 0
    x = [0 for i in range(n)]
    while k<N:
        k += 1
        for i in range(n):
            # 注释部分误差大
            # sum = A[i][i]*x0[i] # 省略掉下面的if判断
            # for j in range(n):
            #     sum += A[i][j]*x0[j]
            sum = 0
            for j in range(n):
                if i != j:
                    # sum += A[i][j]*x0[j] # 表示使用雅可比迭代
                    sum += A[i][j]*x[j] # 表示使用高斯赛德尔迭代
            x[i] = (b[i]-sum)/A[i][i]

        if norm(x,x0)<e:
            return x
        for i in range(n):
            x0[i] = x[i]
        if is_print:
            print('x',k,' = ',x0,sep='')
    print('已超过最大迭代次数')

A = [[10,-1,-2],[-1,10,-2],[-1,-1,5]]
b = [72,83,42]
print(jacobi(A, b,[0,0,0],True))

数据拟合

线性回归

def linearRegression(x,y):
    '''
    线性回归
    x: 横坐标数据，List
    y: 纵坐标数据，List
    return: (a,b) 即y=a+bx
    '''
    n = len(x)
    avg_x = sum(x)/n
    avg_y = sum(y)/n
    sum1 = 0
    sum2 = 0
    for i in range(n):
        sum1 += (x[i]-avg_x)*(y[i]-avg_y)
        sum2 += (x[i]-avg_x)*(x[i]-avg_x)
    b = sum1/sum2
    a = avg_y-b*avg_x
    return (a,b)

# print(linearRegression([0,1,2],[1,3,2]))

一般最小二乘拟合

def normal_equation(x,y,funcs,w=None,is_print=False):
    '''
    利用最小二乘拟合来拟合x,y数据
    funcs 为需要拟合的函数组
    w：为每点的权，默认为1
    '''
    n = len(funcs)
    m = len(x)
    # 初始化
    if not w:
        w = [1 for i in range(m)]
    A = [[0 for j in range(n)] for i in range(n)]
    b = [0 for i in range(n)]
    for i in range(n):
        for j in range(i+1):
            sum = 0
            for t in range(m):
                sum += w[t]*funcs[i](x[t])*funcs[j](x[t])
            A[i][j]=A[j][i] = sum
        sum = 0
        for t in range(m):
            sum += w[t]*funcs[i](x[t])*y[t]
        b[i] = sum
    if is_print:
        print('A =',A)
        print('b =',b)
    return column_pivot_element_gauss(A, b)

# x = [-3,-2,-1,0,1,2,3]
# y = [4,2,3,0,-1,-2,-5]
# print(normal_equation(x, y, [lambda x:1,lambda x:x,lambda x:x*x],None,True))

# x = [-3,-1.5,0,1.5,3]
# y = [-0.1385,-2.1587,0.8330,2.2774,-0.5110]
# print(normal_equation(x, y, [lambda x:math.sin(x),lambda x:math.cos(x)],None,True))

数值微积分

三点微分公式

def three_point_differential(y0,y1,y2,h):
    '''
    三点微分公式求解，其中x2-x1=x1-x0=h
    '''
    coefficient = 0.5/h
    df1 = coefficient * (-3*y0+4*y1-y2)
    df2 = coefficient * (-y0+y2)
    df3 = coefficient * (y0-4*y1+3*y2)
    return (df1,df2,df3)

复化求积公式

def complexification_trapezoid(func,a,b,n):
    '''
    复化梯形求积公式
    func: 被积函数
    求积区间[a,b] 分为n等分
    '''
    sum = 0.5 * (func(a)+func(b))
    h = (b-a)/n
    for i in range(1,n):
        sum+= func(a+h*i)
    return h*sum


def complexification_simpson(func,a,b,n):
    '''
    利用复化Simpson求积公式
    func: 被积函数
    把区间[a,b] 分为2n等分
    '''
    h = (b-a)/(2*n)
    sum = func(a)+func(b)
    for i in range(1,n):
        sum += 2*func( a+h*2*i)
    for i in range(n):
        sum += 4*func(a+h*(2*i+1))    
    return h*sum/3

def test_func(x):
    return 4/(1+x*x)

print(complexification_trapezoid(test_func, 0, 1,8))
print(complexification_simpson(test_func, 0, 1,4))

变长求积公式

def elongate_simpson(func,a,b,e=1e-6,is_print=False,N=100):
    '''
    变步长梯形求积公式
    func被积函数，求积区间[a,b],e代表精度，N最大迭代次数
    '''
    def recursion(func,a,h,n,T1):
        '''
        求解T的递推公式 T(h) = 1/2 * T(2h) + h*sum(f(2k-1),k=1,n/2)
        func 被积函数，a积分下限，
        '''
        sum = 0
        for i in range(1,n,2):
            sum += func(a+h*i)
        return 0.5*T1 + h*sum
    T1 = complexification_trapezoid(func, a, b, 1) # 求初始值
    h = (b-a)/2
    k = 1
    c = 1/3
    T2 = recursion(func, a, h, 2, T1)
    S1 = T2 +  c*(T2-T1)
    if is_print:
        print('S',k,' = ',S1,sep='')
    T1 = T2 
    while k<N:
        h /= 2
        k += 1
        T2 = recursion(func, a, h, 2**k, T1)
        d = c*(T2-T1)
        S2 = T2 + d
        d = abs(S2-S1)/15
        if is_print:
            print('S',k,' = ',S2,sep='')
            print('误差估计：',d)
        if d < e:
            return S2
        T1 = T2
        S1 = S2    
    print('已超过最大迭代次数')

def elongate_test_func(x):
    if x == 0:
        return 1
    return math.sin(x)/x

def elongate_trapezoid(func,a,b,e=1e-6,is_print=False,N=100):
    '''
    变步长梯形求积公式
    func被积函数，求积区间[a,b],e代表精度，N最大迭代次数
    '''
    T1 = complexification_trapezoid(func, a, b, 1) # 求初始值
    h = b-a
    if is_print:
        # print('T(',h,') = 0.5*[ f('a')
        print('T(',h,') = ',T1,sep='')
    k = 0
    c = 1/3
    while k<N:
        sum = 0
        k += 1
        h /= 2
        for i in range(1,2**k,2):
            sum += func(a+h*i)
        T2 = 0.5*T1 + h*sum
        d = abs(c*(T2-T1))
        if is_print:
            print('T(',h,') = ',T2,sep='')
            print('误差估计：',d)
        if d < e:
            return T2
        T1 = T2    
    print('已超过最大迭代次数')
 

# print(elongate_trapezoid(elongate_test_func, 0, 1,1e-3,True))
# print(elongate_simpson(elongate_test_func, 0, 1,1e-6,True))

龙贝格(Romberg)求积法

def test_func(x):
    if x == 0:
        return 1
    return math.sin(x)/x

def romberg(func,a,b,e=1e-3,is_print=False,N=100):
    '''
    龙贝格(Romberg)数值积分方法
    '''
    r = []
    k = 0
    r.append([complexification_trapezoid(func, a, b, 1)])
    while k<N:
        k+=1 
        tmp = []
        tmp.append(complexification_trapezoid(func,a,b,2**k))
        for i in range(1,k+1):
            tmp.append(tmp[i-1]+(tmp[i-1]-r[k-1][i-1])/(4**i-1))
        r.append(tmp)
        if abs(tmp[k]-r[k-1][k-2])<e:
            if is_print:
                for item in r:
                    print(item)
            return tmp[k]
    print('已超过最大迭代次数')

print(romberg(test_func, 0, 1,0.5*1e-5,True))