基于Python的数值分析方法应用_数值计算方法应用实例

引言

本次数值分析编程内容包括利用共轭梯度法求解大规模稀疏方程组、最小二乘拟合问题的求解、非线性方程组的迭代解法、工程应用（牛顿分段二次插值多项式进行工程用表插值计算、共轭梯度法求解超静定问题方程组），编程语言选择为Python。

*注：完整PDF文档及程序请参考：https://gitee.com/zhang-maojie/programming-and-operation.git，如有侵权，联系删除。

1. 利用共轭梯度法求解大规模稀疏方程组

共轭梯度法是把求解线性方程组Ax=b（A是对称正定矩阵）的问题转化为求解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩阵A的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代n次（n为矩阵A的阶数），就可求得二次函数的极小点，也就求得了线性方程组Ax=b的解。共轭梯度法求解的公式如下所示。

1.1程序使用方法

在程序设计时主要考虑了系数矩阵A，向量b，初始迭代值x，精度e，计算结果最大小数点位数pot，如下所示：

1. #共轭梯度法
2. class CGM():
3.     def __init__(self, A: np.ndarray, b: np.ndarray, x=[], e=0.5*10**(-2), pot=16):
4.         """
5.         初始化共轭梯度法！
6.         :param A: 创建实例时需传入Ax=b的矩阵A与向量b,numpy类型；
7.         :param b: 创建实例时需传入Ax=b的矩阵A与向量b,numpy类型；
8.         :param x: 初始向量x可以被修改，默认x(0)=[0,0,0……],列表类型；
9.         :param e: 精度e默认为0.5*10^(-2),可以修改；
10.         :param pot:小数点位数，默认16位，最大16位。
11.         """

在附件的ConjugateGradientMethod.py文件中，主要包含了一个用于创建共轭梯度法求解方程组的类CGM()，其输入参数如上述内容所示，主要包括了共轭梯度法计算中的迭代过程的类方法forward()，以及对已经创建实例的相关参数做修改的类方法params_change()。其次文件中还包含了两个函数，分别是用于储存计算结果的函数save_ans()和用于生成课后计算实习算例3.2相关矩阵、向量的函数Ab()。文件使用方法如下：

1. #计算实习3.2
2.     print("A的阶数为100时的结果:")
3.     A, b = Ab(100)
4.     cgm3_2 = CGM(A=A, b=b)
5.     ans3_2_100 = cgm3_2.forward()
6.     save_ans("./data/共轭梯度法3_2_100", ans3_2_100)

（1）利用 Ab()函数生成系数矩阵 A 和向量 b，需要参数 n 来确定系数矩阵的阶数，也可根据需要自己定义 A，b，只需要保证 A，b 数据类型是 np.ndarrary、数据为真实有效数据就行；

1. #3.2 题的矩阵生成函数
2. def Ab(n): 3. """ 4. 生成系数矩阵 A 和值 b 的方法。
5. :param n: 系数矩阵 A 的阶数。
6. :return:A,b
7. """

（2）利用 CGM()类生成计算实例，生成实例时必须要将参数系数矩阵 A，向量 b 传入，若不做特殊修改，初始向量默认为 x=[0,0,0……]T、精度默认为 0.5×10-2，小数点位数默认为 pot=16；

（3）因为在类定义时将计算和初始化分开考虑，所以需要求得计算结果还需要调用类方法 forward()，该方法会返回一个 np.ndarray 类型的 x 向量，在调用时格式如下：接收参数名称 = 实例名称.forward()，运行之后调试窗口可以打印出计算结果；

（4）如果还需要将结果保存，可以调用 save_ans()函数，其包括文件保存位置及名称参数 path 和计算结果参数 ans（即第三步所用的接收参数名称），此函数可以自动将求解结果以 csv 文件格式保存在对应位置。

1.def save_ans(path, ans):
2.    """
3.    保存结果的函数。
4.    :param path: 文件名及路径,str类型
5.    :param ans: 传入结果。
6.    :return: 无
7.    """

其次程序还设计了相关错误文字提示，使用者可以根据相关情况描述解决一些程序问题。

1.2程序框图

共轭梯度法结合程序后的程序框图如下：

图1 共轭梯度法的程序框图

1.3计算实例

计算实习3.2：用共轭梯度法求解线性方程组Ax=b，其中，

矩阵A的阶数分别取为100,200,400，指出计算结果是否可靠。

各阶计算调用程序如下：

1.#计算实习3.2
2.    print("A的阶数为100时的结果:")
3.    A, b = Ab(100)
4.    cgm3_2 = CGM(A=A, b=b)
5.    ans3_2_100 = cgm3_2.forward()
6.    save_ans("./data/共轭梯度法3_2_100", ans3_2_100)
7.
8.    print("A的阶数为200时的结果:")
9.    A, b = Ab(200)
10.    cgm3_2.params_change(A=A, b=b)
11.    ans3_2_200 = cgm3_2.forward()
12.    save_ans("./data/共轭梯度法3_2_200", ans3_2_200)
13.
14.    print("A的阶数为400时的结果:")
15.    A, b = Ab(400)
16.    cgm3_2.params_change(A=A, b=b)
17.    ans3_2_400 = cgm3_2.forward()
18.    save_ans("./data/共轭梯度法3_2_400", ans3_2_400)
19.
20.    print("A的阶数为400时的结果（加大精度后）:")
21.    A, b = Ab(400)
22.    cgm3_2.params_change(A=A, b=b, e=0.5*10**(-3))
23.    ans3_2_400 = cgm3_2.forward()
24.    save_ans("./data/共轭梯度法3_2_400(高精度)", ans3_2_400)

由于求解结果在程序调试窗口中较长，只做截取展示部分内容并用相关保存文件的数据进行比较。

下面图1展示了几种情况csv文件中的计算结果前17项（其中(a)是n=100、默认精度的结果，(b)是是n=200、默认精度的结果，(c)是n=400、默认精度的结果，(d)是n=400、精度修改为默认值的十分之一的结果。）：

图2 共轭梯度法计算实习3.2结果对比

通过理论分析可以知道本题的标准解为x=[1,1,1……]T，图一中的解基本上趋近这一标准解，对于n=100和n=200的(a)、(b)其结果是基本上可以认为是可靠的；对于n=400分别使用了两种不同的精度进行计算，通过计算发现低精度计算结果(c)中存在诸如0.925这种偏差较大的结果，而通过提高精度的(d)中则改善了这种情况，因此(c)的结果不如(d)可靠。

通过上述计算，可以知道，随着矩阵A阶数的增大在保持计算精度不变的情况下会使得计算结果偏差变大。并且通过试验发现，在精度保持不变时，随着阶数n的增加，其计算结果会变得极不可靠，在(c)后续结果中甚至存在0的情况。

通过分析可以知道，迭代精度不变，随着矩阵阶数n的增大，矩阵的条件数也在增大，条件数大的矩阵为病态矩阵，相应方程组为病态方程组，病态问题在数据很小变化时会令解产生巨大变化，因此在迭代过程中，由于小数的舍入，会导致病态方程组的结果变化剧烈而产生很大误差。

为了用实例证明上述内容，下面分别将算例矩阵阶数调大到500、1000、2000，并用默认精度计算，部分计算结果如下（185~202项）：

图4 阶数为500、1000、2000的计算结果

通过上述计算结果的对比分析可以知道，在精度不变的前提下，随着矩阵阶数的提升，其相关计算的内容会出现“死点”，即无论怎样增大矩阵阶数，其计算结果相对来说不再改变，对于高阶的矩阵的结果中间内容全部为0，并且呈现首尾对称，中间为0的现象。

2. 最小二乘拟合问题的求解

设f(x)是定义在点集X={x1,x2,…,xm}上的列表函数或者f(x)是定义在区间[a,b]上表达式复杂的连续函数，通过构造广义多项式

使得

达到最小，其中c0,c1,…,cn为待定参数，φi(x)(i=0,1,…,n)为已知的一组线性无关的函数（基函数），取p(x)作为f(x)的近似表达式就是最优平方逼近。

当f(x)为列表函数时，设在点xi处的函数值yi，则：

误差表达式为：

当f(x)为[a,b]的连续函数时：

误差表达式为：

ω等于1时的正规方程组如下：

由(8)式可以求得系数向量c，进而完成最小二乘拟合，其中G是由基函数确定的矩阵。

2.1程序的使用方法

1.#最小二乘拟合
2.class LSF():
3.    def __init__(self, dispersed=True, polynomial=True, x=[], y=[], step=0.001, pot=11, **kwargs):
4.        """
5.        这是最小二乘拟合。
6.        :param dispersed:列表函数还是连续函数，默认列表函数，若为连续函数则需在kwargs中输入相关表达式。
7.        :param polynomial: 是否为多项式拟合，选False时激活非多项式拟合的基函数的使用。
8.        :param x: 传入拟合数据x，列表
9.        :param y: 传入拟合数据y，列表
10.        :param step: 拟合曲线的步长，一般情况无需调整，可以调大观察区别，一般拟合首位不适配就需要调整此参数
11.        :param kwargs: 用于传递各种参数，多项式拟合的最高次幂：highest_power=？
12.                                         非多项式拟合的基函数：fai=[?, ?, ...]
13.                   fai内只能是基本初等函数，且线性无关，可以是(c为已知数字,x为变量,为自己所定义的变量符号)：
14.                                         输入的x必须为np.ndarray类型，且shape=（1，n）
15.                                        x**c,c**x,np.log(x)或np.log10(x)及其组合的其他底数的对数函数,
16.                                         np.sin(x),np.cos(x)……
17.                                         np.arcsin(x) ……
18.                                         需要常数项c时输入np.ones((1, len(x)))
19.                          若想在结果中得到表达式、拟合图像、误差等需要输入自己的字符串型的基本初等函数，如：
20.                                         expression = ['1', 'x**2', 'np.log(x)'……]
21.                                         例子：
22.        LSF(polynomial=False, x=X, y=Y,
23.        fai=[np.ones((1, len(X))), np.array([X])**2, np.log(np.array([X])), np.sin(np.array([X])),
24.                    np.exp(np.array([X]))],
25.               expression=['1', 'x**2', 'np.log(x)', 'np.sin(x)', 'np.exp(x)'])
26.                                persed_expression=“y=np.exp(x) + 1”，连续函数表达式,未知数统一用x表达。
27.                                         persed_ab = [a,b],拟合及绘制范围
28.                                         persed_step = 1000,拟合区间插值节点数。默认1000
29.        """

在附件的LeastSquaresFitting.py文件中，定义了用来进行最小二乘拟合的类LSF()，其包含的功能有：列表函数的多项式最小二乘拟合及任意基函数拟合、连续函数的多项式最小二乘拟合及任意基函数拟合。

各参数的使用及含义如下表所示：

表1 LSF()类的参数及含义

其次LSF()类主要包含了五个类方法：Polynomial(self, n)、painter(self, my_step)、caculate(self, x)、Non_polynomial(self, fai)、save_data(self, path)，一个文件读取函数：red_csv_data(path)。

表2 LSF()的类方法说明及文件函数说明

因此利用本文件内容实现最小二乘拟合的使用步骤如下：

1. 自己定义列表函数数据或者读取csv文件数据形成列表函数数据，或者定义一个连续性函数；

2. 利用LSF()传入相关参数生成实例计算得到拟合结果；

3. 根据需要调用save_data(self, path)函数保存计算结果。

详细操作可以参考附件或者下一小节的计算实例。

2.2程序框图

图5 最优平方逼近的程序框图

2.3计算实例

5个相关实例程序如下：

1.if __name__ == '__main__':
2.    # 5.1数据,多项式拟合
3.    X = [i*0.1 for i in range(1, 10, 1)]
4.    Y = [5.1234, 5.3057, 5.5687, 5.9375, 6.4370, 7.0978, 7.9493, 9.0253, 10.3627]
5.
6.    print("多项式拟合：")
7.    lsf1 = LSF(polynomial=True, x=X, y=Y, highest_power=4, pot=5)
8.    lsf1.save_data("./data/最优平方逼近5_1")
9.
10.    print("**"*20)
11.    # #5.1数据,非多项式拟合
12.    print("非多项式拟合：")
13.    lsf2 = LSF(polynomial=False, x=X, y=Y,
14.               fai=[np.ones((1, len(X))), np.array([X])**2, np.log(np.array([X])), np.sin(np.array([X])),
15.                    np.exp(np.array([X]))],
16.               expression=['1', 'x**2', 'np.log(x)', 'np.sin(x)', 'np.exp(x)'],
17.               pot=5)
18.    lsf2.save_data("./data/最优平方逼近5_1_非多项式拟合")
19.
20.    # 自定义数据计算
21.    print("**"*20)
22.    print("自定义数据：y=100sin(x)/x")
23.    print("多项式拟合")
24.    X_a = np.linspace(-10, 10, num=20)
25.    Y_a = 100 * np.sin(X_a) / X_a
26.    lsf3 = LSF(polynomial=True, x=X_a, y=Y_a, highest_power=10, pot=11)
27.    print("**" * 20)
28.
29.    print("非多项式拟合")
30.    lsf4 = LSF(polynomial=False, x=X_a, y=Y_a,
31.               fai=[np.ones((1, len(X_a))), np.array([X_a])**(2), np.sin(np.array([X_a]))],
32.               expression=['1', 'x**(-1)', 'np.sin(x)'],
33.               pot=11)
34.
35.    print("**" * 20)
36.    xy = red_csv_data("./data/zmj.csv")
37.    lsf5 = LSF(polynomial=True, x=xy[0], y=xy[1], highest_power=10, pot=11)
38.    print("x={}时结果是{}".format(1.88, lsf5.caculate(1.88)))
39.    lsf5.save_data("./data/正态分布数据zmj拟合")
40.
41.    #连续函数的拟合
42.    print("**" * 20)
43.    lsf6 = LSF(dispersed=False, polynomial=True, highest_power=10, pot=6, persed_expression="y=np.sin(x)",
44.               persed_ab=[-2*np.pi, 2*np.pi], persed_step=100)
45.    lsf6.save_data("./data/连续函数sin(x)的最小二乘拟合")

1. 计算实习5.1（多项式拟合）

给定数据如下表所示，求其最小二乘拟合四次多项式及其误差：

根据题目要求获得的四次多项式拟合结果如下：

图6 计算实习5.1的4次多项式拟合曲线

图7 计算实习5.1的4次多项式拟合结果

可见其拟合曲线拟合效果良好，估计误差仅为0.00058，相关表达式已在图中展示。

2. 计算实习5.1（非多项式拟合）

基函数选择为：1，x2，ln(x)，sin(x)，ex。

图8 计算实习5.1的非多项式拟合曲线

图9 计算实习5.1的非多项式拟合结果

通过对比计算实习5.1的4次多项式拟合和这里选择的基函数实现的非多项式拟合，可以发现其拟合效果都比较精确，只是这里所选取的基函数所拟合的曲线的误差为0.00114，大于4次多项式的拟合估计误差。

3. y=100sin(x)/x数据形成的列表函数的多项式及非多项式拟合图像对比

图10 y=100sin(x)/x数据形成的列表函数的10次多项式拟合

图11 y=100sin(x)/x数据形成的列表函数的非多项式拟合

在拟合中，多项式拟合选择最高项次数为10，非多项式拟合的基函数为：1，x2，ex。通过对比可以看出10次多项式的拟合效果较好，而所选基函数的非多项式拟合效果较差。

4. 标准正态分布表的数据拟合

通过读取相关标准正态分布表的部分对应数据，对标准正态分布表进行了拟合，拟合效果如下：

图12 标准正态分布表数值的拟合

标准正态分布表在统计学应用中具有重要作用，通过上述10次多项式拟合，可以看出其拟合效果良好，因此可以直接使用上述类的类方法计算任意x点处的

y值，用于标准正态分布的工程计算。

5. 连续函数sin(x)的10次多项式拟合

图13 连续函数sin(x)的拟合

用过用10次多项式对连续函数sin(x)在[-2π,2π]的最小二乘拟合的上图看出，拟合效果良好，原函数曲线和拟合曲线基本重和。

通过对上述5种情况的分析，可以发现，高次多项式拟合对于各种复杂的函数都具有良好的适应性，而非多项式拟合效果则受到所选基函数的影响较大，究其原因是非多项式拟合的基函数相关函数特性不能够很好的适应拟合函数的函数特性，如：单调性，奇偶性，定义域等。

对于计算实习5.1进行从0~9次和0~30次的多项式拟合观察其误差变化，结果如下：

图14 多项式拟合最高次数与误差关系曲线1

图15 多项式拟合最高次数与误差关系曲线2

通过上图结果可以知道，在多项式最小二乘拟合中，并非最高项次数越高其拟合精度就越高，其和数据本身所服从的函数关系有很大关系，当最高项次数超过一定限制时会出现误差急速偏大的情况，并且随着最高项次数的继续增大，误差随之产生波动。

3. 非线性方程组的迭代解法（牛顿法）

牛顿法是求解非线性方程组十分重要的一种方法，其基本思想是将非线性函数逐次线性化，从而形成一个迭代过程。

将方程组每个方程在点x(k)=(x1(k),x2(k),…xn(k))T处按多元函数泰勒公式展开得到：

将上述内容以矩阵向量形式表示后可以得到此方程组的解x(k+1)。

令 ，可以得到迭代公式如下：

3.1程序的使用方法

1.def one_stage_Newton(n, s, start, equations=[], variable=[], initial_value=[], e=0.01, K=1000, pot=11, save=True):
2.    """
3.    直接一次性实现牛顿非线性方程组迭代求解。
4.    :param n: 变量数目 int型
5.    :param s: 变量标记 str型
6.    :param start: 0 or 1,起始编号 s0 or s1
7.    :param equations: 方程组，以列表新式传递，内部为标准等式组。如：
8.                           “x1**2+exp(x2)-10”……省略“=0”，str型
9.    :param variable:变量标志，如:"x1",'x2'……，str型
10.    :param initial_value:初始迭代值，与变量标志位置对应，默认全为0
11.    :param e:停止的相对误差限,默认0.01
12.    :param K:最大迭代次数
13.    :param pot:计算结果小数点位数，默认11，最大16位
14.    :param save:结果保存标志，默认True，保存
15.    :return: 无
16.    """
 
1.class Newton():
2.    def __init__(self, equations=[], variable=[], initial_value=[], e=None, K=None, pot=11):
3.        """
4.        这是牛顿法求解非线性方程组。
5.        :param equations: 方程组，以列表新式传递，内部为标准等式组。如：
6.                           x1**2+exp(x2)-10……省略“=0”
7.        :param variable:变量标志，如:x1,x2……
8.        :param initial_value:初始迭代值，与变量标志位置对应，默认全为0
9.        :param e:停止的相对误差限
10.        :param K:最大迭代次数
11.        :param pot:计算结果小数点位数，默认11，最大16位
12.        """
 
1.def define_variables(n, s, start):
2.    """
3.    定义变量。直接定义n个以s为基本标记的变量，如s1 s2 s3……
4.    :param n: 变量数目 int型
5.    :param s: 变量标记 str型
6.    :param start:0 or 1,起始编号 s0 or s1
7.    :return:
8.    """