NumAnalysis

数值分析/NumAnalysis

我们希望将《计算方法》课程中的理论内容通过Python语言进行实现，并将其运用到实际问题之中。 作为一个对计算机稍微有一定了解的学习者，深知将理论转化为实际往往是存在困难的。通常发展了一门新的理论，并将其转化为技术，这期间需要很长的时间来逐步完善并提高效率。

更重要的是，我们认为这门课程中的一些方法是形象生动而有趣的。可能在更高级的模块中集成了比这些先进得多的算法，但是这并不重要，我们希望用这本书中的方法来解决一些适合于我们的需求的问题。

This repository is used to create a project on lesson "Numerical Analysis" and practically use the method we have learned in problem solving. As a learner who knows only a little about computer and programming, I know that there must exist difficulties in the transformation from theory to practice. It takes long time when the theory generally improves efficiency as technical tools.

From our perspective, some methods in this lesson are really interesting and fantastic. Maybe some much advanced technique has been used in popular modules, however, it is not important while we hope that the method in our textbook can solve the problem we faced and adapted to our demand.

目录/Table of Contents

• 非线性方程数值解/Numerical Solution to Nonlinear Equation
  • 迭代法/Recursive Method
  • 牛顿法/Newton Method
  • 牛顿-拉夫逊法/Newton-Ralfsnn's method
  • 双点快速截弦法/
• 线性方程组LU分解/System of Linear Equation——LU Decomposition
  • 杜利特尔分解/Doolittle Decomposition
  • 性态分析/Property Analysis
• 线性方程组迭代法/Recursive way to solve the System of Linear Equation
  • Jacobi迭代/Jacobi Recursion
  • Gauss-Seidel迭代/Gauss-Seidel Recursion
• 插值与拟合/Interpolation and Fitting
  • 多项式插值/Polynomial Interpolation
  • Lagrange插值/Lagrangian Interpolation
  • Newton插值/Newton Interpolation
  • 最小二乘拟合/Least Square Fitting
  • 插值与拟合绘图/Plot Illustration
• 数值积分/Numerical Integral
  • 插值型数值积分
  • 牛顿科特斯公式/Newton-Cotes Method
    • 梯形公式/Trapezium
    • Sinpson公式/Sinpson
  • 复化求积公式/Composite Integral
    • 复化梯形公式/Composite Trapezium
    • 复化Sinpson公式/Composite Sinpson
• 常微分方程数值解/Numerical Solution in Ordinary Differential Equation
  • Euler公式/Euler Method
  • 改进的Euler法/Improved Euler Method

非线性方程数值解/Numerical Solution to Nonlinear Equation

迭代法

Recursive Method

基本原理 在这个部分中我们将利用循环完成迭代的操作，需要手动输入迭代公式，具体代码如下：

import numpy as np
import sympy as sp
g = lambda x: x**2-3*x**(1/2)+4 # Here input your recursive function
x0 = 1 # Your initial point
n = 10 # Your iteration times
residual = 0.001 # Your tolerance 

def Recursive_Method(g,x0,n,residual):
    x = np.zeros(n+1)
    x[0] = x0
    for index in range(1,n+1):
        x[index] = g(x[index-1])
        difference = x[index] - x[index-1]
        if np.abs(x[index]-x[index-1])<residual:
            print("Iteration %r: x=%r, difference=%r <= Residual=%r" %(index,x[index],difference,residual))
            break
        else:
            print("Iteration %r: x=%r, difference=%r" %(index,x[index],difference))
    print("Terminal: The final result is %r" % x[index])
    return x[index]

Recursive_Method(g,x0,n,residual)

例如我们要使用求解在上的根，迭代10次或两次误差小于，则使用以下的代码

g = lambda x: 0.5*(10-x**3)**0.5
Recursive_Method(g,1.5,10,1e-5)

得到的结果为：

Iteration 1: x=1.2869537676233751, difference=-0.2130462323766249
Iteration 2: x=1.4025408035395783, difference=0.11558703591620323
Iteration 3: x=1.3454583740232942, difference=-0.057082429516284172
Iteration 4: x=1.3751702528160383, difference=0.029711878792744173
Iteration 5: x=1.3600941927617329, difference=-0.015076060054305396
Iteration 6: x=1.3678469675921328, difference=0.0077527748303998223
Iteration 7: x=1.3638870038840212, difference=-0.0039599637081115802
Iteration 8: x=1.3659167333900399, difference=0.0020297295060187626
Iteration 9: x=1.3648782171936771, difference=-0.0010385161963628597
Iteration 10: x=1.3654100611699569, difference=0.00053184397627981106
Terminal: The final result is 1.3654100611699569

注意

1. 本部分尚未加入是否发散的判断，因此要根据迭代结果进行发散的判断
2. 请自行计算迭代式，在输入过程中注意指数与根号的输入规范

牛顿法

Newton Method

具有特殊迭代格式的迭代法，对于的求解，采用格式：

能保证更快的收敛速度。这里，我们为了得到更为精确的结果，避免数值微分带来的误差和符号计算可能碰到的不确定性，要求手动输入函数，并同理使用迭代法中的方法：

f = lambda x: x**2+3*x+1 # Your function f(x)
f1 = lambda x: 2*x+3 # First order deviation of f(x)
x0 = 1 # The initial point
n = 100 # Iteration time
residual = 1e-5 # Tolerance

def Newton_Recursive(f,f1,x0,n,residual):
    x = np.zeros(n+1)
    x[0] = x0
    for index in range(1,n+1):
        x[index] = x[index-1]-f(x[index-1])/f1(x[index-1])
        difference = x[index] - x[index-1]
        if np.abs(x[index]-x[index-1])<residual:
            print("Iteration %r: x=%r, difference=%r <= Residual=%r" %(index,x[index],difference,residual))
            break
        else:
            print("Iteration %r: x=%r, difference=%r" %(index,x[index],difference))
    print("Terminal: The final result is %r" % x[index])
    return x[index]

Newton_Recursive(f,f1,x0,n,residual)

作为对比，我们依旧解决迭代法中10次尚未达到精度的例子：

求解在上的根，迭代10次或两次误差小于。 此时，在牛顿法中，。对应代码为：

f = lambda x: x**3+4*x**2-10
f1 = lambda x: 3*x**2+8*x
Newton_Recursive(f,f1,1.5,10,1e-5)

仅通过三次迭代就得到收敛结果

Iteration 1: x=1.3733333333333333, difference=-0.12666666666666671
Iteration 2: x=1.3652620148746266, difference=-0.0080713184587066777
Iteration 3: x=1.3652300139161466, difference=-3.2000958479994068e-05
Iteration 4: x=1.3652300134140969, difference=-5.0204973511824846e-10 <= Residual=1e-05
Terminal: The final result is 1.3652300134140969

牛顿-拉夫逊法

Newton-Ralfsnn Method

非线性方程组的求解是一个更加困难的问题，非线性方程组的解往往具有多个，也可能不存在，并且对初值十分敏感，我们在这里使用牛顿-拉夫逊法进行非线性方程组的迭代求解。这种方法是牛顿法的推广，但存在很多的局限性。

，其中为Jacobi矩阵（全体一阶偏导矩阵）

从简化的角度，我们仅分析三个方程组三个未知数的情形，解决代码如下：

def Jacobi_Matrix(f1,f2,f3,ak,bk,ck):
    a = sp.Symbol('a')
    b = sp.Symbol('b')
    c = sp.Symbol('c')
    F = [f1,f2,f3]
    X = [a,b,c]
    n = len(F)
#    F_1 = np.array([[a,b,c],
#                    [a,b,c],
#                    [a,b,c]])
    F_1v = np.zeros((n,n))
    for i in range(n):
        for j in range(n):
#            F_1[i,j] = sp.diff(F[i],X[j])
            storage = sp.diff(F[i],X[j])
            F_1v[i,j] = storage.evalf(subs = {a:ak,b:bk,c:ck})
    #F1v = F_1.evalf(subs = {X:X0})
    return np.linalg.inv(F_1v)

def Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n):
    a = sp.Symbol('a')
    b = sp.Symbol('b')
    c = sp.Symbol('c')
    x = np.zeros((3,n+1))
    x[0,0] = a0
    x[1,0] = b0
    x[2,0] = c0
    for k in range(n):
        if np.min(np.abs(x[:,k+1]-x[:,k]))>100:
            return "Error with divergence"
#            print("Error")
        else:    
            if np.max(np.abs(x[:,k+1]-x[:,k]))>1e-5:
                F1 = f1.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
                F2 = f2.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
                F3 = f3.evalf(subs = {a:x[0,k],b:x[1,k],c:x[2,k]})
                F1 = float(F1)
                F2 = float(F2)
                F3 = float(F3)
                F = np.array([F1,F2,F3])        
                x[:,k+1] = x[:,k]-Jacobi_Matrix(f1,f2,f3,x[0,k],x[1,k],x[2,k])@F
            else:
                return x[:,:k+1]
#                print("Error")
    return x

在这里我们默认两次迭代之间最小的变化的绝对值超过100为收敛条件（这可能不是一个好的判断方法），当两次迭代最大变化绝对值小于1e-5则认为迭代法收敛。一个例子如下：

# f = [f1,f2,f3,f4,f5]'=0
# a* = 1 b* = 3 c* = 2 
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
f1 = a**2 + b - 2*c
f2 = 3*a + 2**b -c**3 - 3
f3 = a + b + c - 6
a0 = 10
b0 = 25
c0 = 80
n = 100
solution = Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n)

我们得到在初值为时的迭代结果为（尽管在我构造这个非线性方程组时设定的解是，下面我们再使用另一组初值进行迭代：

# f = [f1,f2,f3,f4,f5]'=0
# a* = 1 b* = 3 c* = 2 
a = sp.Symbol('a')
b = sp.Symbol('b')
c = sp.Symbol('c')
f1 = a**2 + b - 2*c
f2 = 3*a + 2**b -c**3 - 3
f3 = a + b + c - 6
a0 = 1
b0 = 2.5
c0 = 3.3
n = 100
solution = Newton_Recursive_Multi(f1,f2,f3,a0,b0,c0,n)

这种情况下我们得到了我们想要的收敛结果。

双点快速截弦法/

线性方程组LU分解/System of Linear Equation——LU Decomposition

本部分求解的线性方程组为行列式不为0的实方阵，用数学符号来表示，为：

其中

杜利特尔分解

Doolittle Decomposition

我们考虑将系数矩阵分解为一个同维度的下三角矩阵和一个同维度的上三角矩阵的乘积，即

之所以这样做，是因为对于的求解是相对容易的，可以