啥？每一期出作业辅导！不不不，我怂！这件事的起因是三次抽签工程数学都没中，第一次，二志愿落选，第二次，40进17落选，第三次，7进2落选。当时有一丢丢小情绪，所以自己捣鼓一下。

我怂.png

今天打算对这个系列收尾了，不管自己肝到哪，都先告一段落了。一是攒了几天的作业还没做，二是我一个搞传感器的，天天写代码，这属于不务正业。

当我再编写程序时，我其实只是在扮演一个解决问题的执行者，所有公式和方法都是已知的，我只是按照已有的步骤一步步执行下去，这是一个比较low的工作。而进行科学研究，扮演的是问题的解决者，面对一个问题时，很多情况根本没有现有的解决方法，需要科研工作者去探索和试错，这是更难的工作。
圆规正转，前两节主要关于线性方程组的求解，本节打算研究非线性方程求解和插值。同样，能力有限，请多多批评指正。
第一节：https://www.jianshu.com/p/671a94ce586b
第二节：https://www.jianshu.com/p/05ae9194068a

三、非线性方程组求解

1.二分法

《数值方法》（关治等）新书终于到了，今天参考该书的96页到97页。并用例题4.2.1来验证。
这个方法的思路是比较简单的。编写代码的时候要确定运算步数，这里利用误差限来确定。推导如下图所示。为什么要这样推导，主要是为了编写代码方便，因为写一个以2为底的log函数不方便。

运算步数.png

以下是我的代码：

def Half(a,b,Epsilon,isPrint=False):
    '''
    a ：左边端点
    b ：右边端点
    Epsilon ：Epsilon误差
    isPrint ：布尔值，是否打印中间过程
    '''
    mid = (a + b)/2
    k = 1
    n = 1 + int((np.log10(b-a)-np.log10(Epsilon))/np.log10(2)) # 计算需要运算的步数
    print('在这个误差限下需要运算%d次' % n)
    while k <= n:
        if isPrint == True:
            print ("第%d步,左端点a:%f,右端点b:%f,此时x=%f,f(x)=%f" % (k,a,b,mid,f(mid)))
        if f(mid) == 0:                              # f()为所需求解的函数
            return mid
        elif f(a)*f(mid) < 0:
             b = (a+b)/2
        else:
             a = (a+b)/2
        mid = (a+b)/2
        k = k + 1
    return mid

其中，函数里面有一个f(x)函数，这个是认为预先编写好的要求解的函数，以书上例题4.2.1为例。

import numpy as np
def f(x):
    return x**3-15*x**2+319

运行测试一下。

x = Half(5,10,0.000005)
'''
在这个误差限下需要运算20次
x = 5.930750370025635
'''
#其中，当isPrint=True时，可以打印中间过程
x = Half(5,10,0.000005,True)
'''
在这个误差限下需要运算20次
第1步,左端点a:5.000000,右端点b:10.000000,此时x=7.500000,f(x)=-102.875000
第2步,左端点a:5.000000,右端点b:7.500000,此时x=6.250000,f(x)=-22.796875
第3步,左端点a:5.000000,右端点b:6.250000,此时x=5.625000,f(x)=22.369141
第4步,左端点a:5.625000,右端点b:6.250000,此时x=5.937500,f(x)=-0.488525
第5步,左端点a:5.625000,右端点b:5.937500,此时x=5.781250,f(x)=10.883087
第6步,左端点a:5.781250,右端点b:5.937500,此时x=5.859375,f(x)=5.181545
第7步,左端点a:5.859375,右端点b:5.937500,此时x=5.898438,f(x)=2.342397
第8步,左端点a:5.898438,右端点b:5.937500,此时x=5.917969,f(x)=0.925885
第9步,左端点a:5.917969,右端点b:5.937500,此时x=5.927734,f(x)=0.218415
第10步,左端点a:5.927734,右端点b:5.937500,此时x=5.932617,f(x)=-0.135122
第11步,左端点a:5.927734,右端点b:5.932617,此时x=5.930176,f(x)=0.041630
第12步,左端点a:5.930176,右端点b:5.932617,此时x=5.931396,f(x)=-0.046750
第13步,左端点a:5.930176,右端点b:5.931396,此时x=5.930786,f(x)=-0.002561
第14步,左端点a:5.930176,右端点b:5.930786,此时x=5.930481,f(x)=0.019534
第15步,左端点a:5.930481,右端点b:5.930786,此时x=5.930634,f(x)=0.008486
第16步,左端点a:5.930634,右端点b:5.930786,此时x=5.930710,f(x)=0.002962
第17步,左端点a:5.930710,右端点b:5.930786,此时x=5.930748,f(x)=0.000200
第18步,左端点a:5.930748,右端点b:5.930786,此时x=5.930767,f(x)=-0.001181
第19步,左端点a:5.930748,右端点b:5.930767,此时x=5.930758,f(x)=-0.000490
第20步,左端点a:5.930748,右端点b:5.930758,此时x=5.930753,f(x)=-0.000145
'''

结果与教材基本一致，在误差限为0.000005的情况下，需要运算20步。

2.不动点迭代法

这个计算方法也比较简单，但要先将f(x)=0认为地转化为x=g(x)的形式。这种转换不具有唯一性，但有些转换无法收敛，且迭代到全局最优点的速度也不同。这里参考教材《数值方法》（关治等）的4.3节。
以下是我的代码：

def fixedPointIteration(x0, Epsilon, N):
    '''
    x0：初始值
    Epsilon：表示迭代后期相邻两次结果的差值，这个值可以尽量设小一点
    N：最大的迭代次数
    
    下方的g(x)是提前定义好的
    '''
    for k in range(N): 
        print("x(%d)=%.10f"%(k,x0))
        x1 = g(x0)                 # 利用g(x)开始迭代
        if abs(x1-x0) < Epsilon:
            print("x(%d)=%.10f"%(k+1,x1))
            print('成功收敛！')
            return x1,k+1        # 最终的x1为所求的解
        x0 = x1                    # 准备下一次迭代

    print('没有收敛！')
        
    return x1,k+1

利用书上例题4.3.1来做测试，x的初始值选为1.50

def g(x):
    return (10/(4 + x))**0.5
x = fixedPointIteration(1.5, 0.00000001, 15)
'''
x(0)=1.5000000000
x(1)=1.3483997249
x(2)=1.3673763720
x(3)=1.3649570154
x(4)=1.3652647481
x(5)=1.3652255942
x(6)=1.3652305757
x(7)=1.3652299419
x(8)=1.3652300225
x(9)=1.3652300123
x(10)=1.3652300136
成功收敛！
'''

可以看到，和书上的结果是一致的。
3.牛顿法
此处参考教材4.5节。考虑牛顿迭代法的一般情况,也就是不考虑重根的情况，那就非常简单了。
此处要求给定的f(x)的导函数，不是一个数值，二是导函数。这在MATLAB上是可以比较方便实现的，Python也有一个库“sympy”也可以实现，但限定用numpy，所以直接认为定义导函数。实际上，用sympy库写也没有简单多少。
以下是我的代码：

def newton_iteration(x0, Epsilon, N):
    '''
    x0:初始值
    Epsilon:误差容限，跟前面一样
    N:最大迭代数
    Y:所求解的函数，提前定义
    diff_Y:所求解函数的导数，提前定义
    '''
    for k in range(N): 
        print("x(%d)=%.10f"% (k,x0) )
        x1 = x0 - Y(x0)/diff_Y(x0)    # 牛顿法的迭代公式
        if abs(x1 - x0) < Epsilon:
            print("x(%d)=%.10f"%(k+1,x1))
            print("成功收敛！")
            return x1  # x1为所求的解
        x0=x1
    return x1

用教材例题4.5.1来做测试。初始值为1.50

def Y(x):
    return x**3+4*x**2-10
def diff_Y(x):
    return 3*x**2+8*x
x = newton_iteration(1.5,0.00001,10)
'''
x(0)=1.5000000000
x(1)=1.3733333333
x(2)=1.3652620149
x(3)=1.3652300139
x(4)=1.3652300134
成功收敛！
'''

可以看到与书上结果一致，并且确实比二分法迭代收敛速度快。

四、插值法

1.拉格朗日插值

插值法的代码编写有一个小困难，因为插值之后得到的是一个函数，而不是一个数值。这个函数是n项的和，每一项是一个系数乘上（x-xi），其中n表示n次拉格朗日插值。具体可以参考《数值方法》（关治等）的6.1章节。
因此，这里我编写代码时，输出是某一个点在这个插值函数的值。如果想要得到插值函数的系数，可以使用scipy的相关库。
以下是我的代码：

def Lagrange(X, Y, Point):
    '''
    X：X(0),X(1),...,X(n)，list或array形式，表示选择用于插值的点
    Y：Y(0),Y(1),...,Y(n)，list或array形式
    Point：对于n次拉格朗日插值，代入某点。目的是求这一点的拉格朗日插值结果
    '''
    result=0.0
    for i in range(len(Y)):
        temp = Y[i]
        for j in range(len(Y)):
            if i != j:
                temp *= (Point-X[j])/(X[i]-X[j])  # 插值基函数的值乘于Y(i)的值
        result += temp  # 这n项的值累加起来
    return result

同样，用教材例题6.1.1的三个小问来验证这个代码。因为上述代码输出并非一个函数而是一个数值，但是教材上的结果是一个函数。因此，这里选择做图展示，来验证代码。
先生成数据

X1 = [0.2,1.0]
Y1 = np.cos(X1)
LaX1 = np.arange(0,1.2,0.01)
LaY1 = [Lagrange(X1,Y1,point) for point in LaX1]
X2 = [0.0,0.6,1.2]
Y2 = np.cos(X2)
LaX2 = np.arange(0,1.2,0.01)
LaY2 = [Lagrange(X2,Y2,point) for point in LaX2]
X3 = [0.0,0.4,0.8,1.2]
Y3 = np.cos(X3)
LaX3 = np.arange(0,1.2,0.01)
LaY3 = [Lagrange(X3,Y3,point) for point in LaX3]

之后开始做图。左边这一列是例题三个小问中选取的节点，右边这一列是插值之后的效果图。这里是对cosx函数进行插值，选取节点越多，得到的插值函数图像应该越接近cosx函数图像。

plt.figure()
plt.subplot(321)
plt.scatter(X1,Y1)
plt.subplot(322)
plt.plot(LaX1,LaY1,'-')
plt.plot(LaX1,np.cos(LaX1),'-')
plt.subplot(323)
plt.scatter(X2,Y2)
plt.subplot(324)
plt.plot(LaX2,LaY2,'-')
plt.subplot(325)
plt.scatter(X3,Y3)
plt.subplot(326)
plt.plot(LaX3,LaY3,'-')
plt.tight_layout()
plt.show()

结果如下图所示：

结果展示.png

可以看到选取的插值节点越多，得到的插值函数越接近我们想要的情况。其中，第一行右边的图黄色的线表示cosx的函数图像，蓝色的线表示插值函数图像，可见，在选取两个节点的情况下，插值误差大。

2.Hermite插值

下午15:,30，累了。长时间写代码容易掉头发，而且久坐工位容易腰痛和大腹便便，年轻人就该多运动。暂告一段落了。

今日小结

人生的执行者还是探索者？也许，人生的一大乐趣就是你有机会做一些无用的事。