Skip to content

Latest commit

 

History

History
774 lines (640 loc) · 35 KB

File metadata and controls

774 lines (640 loc) · 35 KB

计算方法

初闻不识曲中意,再闻已是曲中人 https://zhuanlan.zhihu.com/p/456099201

考点

理解层面

  • 前置知识
  • 线性方程的数值解法
    1. 高斯消去法
    2. 矩阵分解
    3. 迭代法性质
    4. 雅可比迭代法 & 高斯-塞德尔迭代法
  • 插值与拟合方法
    1. 曲线拟合方法
    2. 最佳逼近理论
  • 优化方法
    1. 最优化问题
    2. 凸集/凸函数/凸优化

掌握层面

  • 非线性方程的数值解法
    1. 二分法
    2. 牛顿迭代法(切线法)
    3. 弦截法(离散牛顿法)
    4. 广义迭代法
  • 插值与拟合方法
    1. 多项式插值
    2. 分段插值
  • 优化方法
    1. 黄金分割法
    2. 二分法

实操层面

  • 迭代实现的计算器方法

PROPERTIES

前置知识

  我们的故事始于方程求根这件“小事”,尤其是对于一个复杂的非线性方程,我们往往只能借助于计算机来实现数值解的求解。而这一过程在计算机中又是以什么样的逻辑进行实现的呢?这也就是这本书想要和我们探讨的话题……

  首先,回顾一下我们已经拥有的toolbox,譬如:求根公式、因式分解、代入消元、复数表示等。而如上所述往往局限于低次方程的求解,而无法胜任五次及以上的代数方程的求解。在此基础上,更别提规模庞大的矩阵方程(多元方程组)的求解了。尽管当上升到方程组层面,我们已有的工具包括但不限于Cramer法则、各种消元法等。

  既然如此,计算机科学的发展便为我们toolbox的扩容承担起了开辟全新道路的光荣使命。换言之,一系列数值算法的引入极大地缓解了现代科学中所涉及到的一切有关方程(组)求解的计算量问题。

  那么……就让我们开始吧!


  首先,一些前置知识的铺垫对于接触一个全新的知识体系显得十分必要。而范数这一概念的学习,便是计算方法这门学科前置知识中的的前置知识。读者很可能想要问为什么,绝对值这个东西我们难道不是早就学透了吗?其实不然。这里面可有一些门道……

  我们都知道,在大物实验中衡量一次测量好坏的最好方法就是算出它的误差,我们也学过了一系列的误差表示方式以及衡量手段。而在计算方法中,我们统一规定了一套标准来衡量计算机在进行一次“观测”——求方程的某个数值解时的performance,这套标准就是范数,尤其是浓缩着海量信息的向量和矩阵的范数。

  类似于定义一个元素的范数,向量范数的概念几乎没什么大的不同,无非是满足一些基本的规定并具有一些相同的性态。但重点来了,它们的唯一区别在于向量范数的家族十分庞大,以至于要用编码的形式来将他们区分开来,譬如:1-范数、2-范数、$\infty$-范数等等。它们也服从着一套统一的定义形式,即$||x||p=(\sum^n{i=1}|x_i|^p)^{1/p}$,说白了就是把向量拆开成一个个单独的元素,进行一波操作后得到一个所谓的整体。

  哦对了,又到了激动人心的时刻,那就是谈谈具体的考法了。就像上面所说的那样,定义一种全新运算的同时往往意味着对现有运算提出了一次全盘考核,看看他们哪些已经具备了这些新的“时代精神”。(这个过程本身其实也是对于向量范数定义的一次推广)这三种新的时代精神也就是非负性齐次性三角不等式

  接下来我们要介绍一种向量范数的既独特又神奇性质——等价性。虽然向量范数的家族十分庞大,但他们之间却存在着密不可分的关系。大体上讲,对于不等式$c_1||x||\beta\leq ||x||\alpha\leq c_2||x||_\beta$而言,可以用夹逼定理理解,即被两个同质序列从两边限制住范围的序列的极限与两侧相同。它想描述的其实也是一种p-范数信息层面的冗余性问题。

  理解之后,我们还需要重点关注由等价性推出的三个不等式,掌握它们的证明方式和记住它们同样重要。 $$||x||_2\leq ||x||1\leq \sqrt{n}||x||2$$ $$||x||\infty\leq ||x||1\leq n||x||\infty$$ $$||x||\infty\leq||x||2\leq \sqrt{n}||x||\infty$$ 类似的,向量范数还有着收敛充分性这一性质,其证明方法与上文夹逼定理的使用方法相同。最后一个性质,也是起着承前启后作用的性质就是压缩性了,这里先不过多展开,具体会在迭代法性质一节中进行详细讨论。

  终于,我们已经盲人摸象般把向量范数这个新东西给大致捋清了。趁热打铁,我们将它推广至多维向量(矩阵)中看看会闹出哪些新名堂。

  在基本的规定与性态方面,二者大体相同,唯一一点就是矩阵范数多了一个矩阵乘法不等式。其实不难理解,两个形态相同的向量是不能作乘法的,而矩阵就可以。而这一不等式其实就是把三角不等式中的加号换成了乘号罢了。除此之外,矩阵范数的定义十分自然,恰好就等于这个矩阵所囊括的所有元素范数之和。

  在向量范数的编码体系的基础上,矩阵范数将其进一步与矩阵本身天然的“图形化界面”相融合。即1-范数、2-范数、$\infty$-范数分别对应起列(和)范数欧几里得范数(谱范数)和行(和)范数。最特殊的显然就是这个“欧几里得”了,它的定义如下$||A||{Euclid}=\sqrt{\lambda{max}(A^TA)}$,看到这里的 $\lambda$ 是否有一种亲切感?有就对了!它就是我们在线性代数中学过的特征值。

  类似的考法在矩阵范数这里重现——用四条规则去“卡”一个不知道从哪里蹦出来的生造的矩阵范数。正是因为如此相似,矩阵范数与向量范数间的相互转换也是具有合法性的,它甚至比矩阵和向量间的相互转换还要来得更直接。

  一个新的数学工具要想踏上实用化的道路,十分矛盾的是,要做的最后一个工作就是把它给抽象化。矩阵范数也不例外。算子范数就是这样一个例子:顾名思义,它所表征的就是一个作为线性算子的矩阵A对于任意一个向量经过线性映射后在前后两种线性空间中的某种范数变化的“拉伸效应”的一种度量。可以说,算子范数搭建起了向量范数与矩阵范数间的桥梁。

  同样,这里我们也需要重点关注几个不等式: $$\frac{1}{\sqrt{n}}||A||2\leq ||A||\infty\leq \sqrt{n}||A||2$$ $$\frac{1}{n}||A||\infty\leq ||A||1\leq n||A||\infty$$ 这里须格外注意避免与向量范数混淆。

  最后,作为范数概念的一个重要补充,并且也是后文函数的内积与正交多项式这一章节的前置知识点,连续函数的范数的概念也很容易类比记忆。此外,向量/矩阵序列的极限为范数的等价性上了最后一层保障——即不同“编码”下收敛的性态是一致的。而谱半径 $\rho$ 也为我们证明矩阵序列收敛提供了良好抓手。它的定义是这样的——给定矩阵全体特征值中的范数最值。(类似于一种“缩放因子”)不难猜到,当这个“缩放因子”小于1时,就注定着这个矩阵序列是收敛的。

  • 谱半径不大于给定矩阵的行列式
  • 实对称矩阵的欧式范数就是它的谱半径

  话题终于来到一切范数的必然归宿——误差。这里不妨提前借用一个之后会学到的知识,不动点迭代的收敛与发散。而发散的迭代序列必然会招致误差的无限扩大。这时,对误差的控制就显得尤为重要。 $$(绝对)误差:E(\tilde x)=x^-\tilde x$$ $$绝对误差极限:\delta\geq |x^-\tilde x|$$ $$ps.The\ error\ limit\ is\ not\ unique.\ It\ is\ usually\ the\ best\ value\ that\ satisfies\ the\ inequality\ relationship.$$ $$相对误差:E_r(\tilde x)=\frac{E(\tilde x)}{\tilde x}$$ $$相对误差极限:\delta_r\geq |\frac{E(\tilde x)}{\tilde x}|$$ $$ps.Since\ the\ exact\ value\ is\ unknown,\ an\ approximate\ value\ is\ used\ instead.$$ $$绝对误差限:\delta(\tilde x)=|x^-\tilde x|$$ $$相对误差限:\delta_r(\tilde x)=\frac{\delta(\tilde x)}{|\tilde x|}$$   在高度依附于计算机的计算方法中,误差主要表现为截断误差(也叫方法误差)与舍入误差两方面。而舍入误差不在我们的考虑范围之内,那是体系结构这门课应该考虑的事情。与大学物理*这门课相同的是,此处的绝对误差、相对误差的定义基本不变,唯一比较玄乎的就是有效数字这个概念的重新定义。我尽可能地把它讲得生动形象些(

  有效数字在计算方法中的定义是这样的:如果一个测量值与精确值间的误差的数量级不超过数位某位的半个单位$(e.g.\ 0.005即百分位的半个单位)$,那么我们就把这个数位在测量值中的某位及其之前的所有非零数字称为该测量值的有效数字。乍听起来有些抽象,这里我们不妨把它与绝对误差限和相对误差限建立起联系:若绝对误差限$\delta(\tilde x)\leq \frac{1}{2}\times10^{k-n}$,则$\tilde x$具有$n$位有效数字。类似地,若相对误差限$\delta_r(\tilde x)\leq \frac{1}{2(a_1+1)}\times10^{1-n}$,则$\tilde x$具有$n$位有效数字;若$\tilde x$具有$n$位有效数字,则$\delta_r(\tilde x)\leq \frac{1}{2a_1}\times10^{1-n}$。

  当然,读者可能会认识到,我们一直在拿单一的数据摆弄来摆弄去,然而迎面而来的却是复杂的函数与方程,那么误差在其中又会有着怎么样的新变化呢?这就不得不谈到误差的$propagation$了。这里我建议大家只用记住一个核心的式子即可,万变不离其宗,$\delta f(\tilde x)\leq |f'(\tilde x)|\cdot\delta(\tilde x)$,在后面我们也会见到体现类似思想的中值定理相关的应用。在此基础上,我们不难推得更为一般的多元函数的误差传播公式,聪明的读者可能已经猜到了不过是把$f(x)$的导数替换为了分别对多个变量的偏导……这里为大家提供几个常用的多元函数模板(四则运算),以供在更复杂的函数中作为基本模块实现多次复用。 $$然而并没有提供……$$

  话不多说,我们言归正传,为了在数值算法的运作过程中始终保持它的数值稳定性,哦对,我们这里要先讲一下数值稳定性这个概念……像那种不动点迭代中发散的迭代公式,我们称之为数值不稳定。相反的,数值稳定这个属性又被细分为了条件稳定(相对稳定)和无条件稳定(绝对稳定),顾名思义,条件稳定往往跟一些具体的条件息息相关——初始值条件、输入数据的性质、算法参数诸如此类。除了要尽可能的使用数值稳定的算法之外,我们也与一些人为可操作的小技巧来帮助我们防止误差扩大过快:

  1. 避免相近数间作差
  2. 避免“小不点”作分母
  3. 避免舍入误差(大数吃掉小数)
  4. 简化计算步骤、改写迭代方程

高斯消去法

  • 精确解法
  • 进行条件:主元素保持非零 (进阶方案:主元素高斯消元法)
  • 误差控制:主元素绝对值过小的问题 (进阶阶方案:多次交换方程顺序)
一般步骤
  1. 向“上三角”靠拢——消元过程 (抽象之处:矩阵表示)
  2. 相反顺序的回代过程 (抽象之处:写出递推公式)
计算量
  1. 乘法
  2. 除法
  3. 总计算量
衍生物

高斯-约当消元法:系数矩阵对角化(不需要回代过程的高斯消去法)

矩阵分解

  • 精确解法
解法分类
  1. 直接三角分解法
  2. Cholesky分解法(对称正定阵可以被分解为一个下三角阵及其转置的乘积) $$A=LL^T$$

不作考察要求? 详见线性方程组-直接法 2:Cholesky分解

  1. 追赶法
    1. Doolittle分解(单位下三角阵与上三角阵的乘积形式)
    2. Crout分解(下三角阵与单位上三角阵的乘积形式)

假迭代,真误差分析

区分一组概念

病态的范围更大 ——即使一个问题是适定的,它也可能仍是病态的

  1. 适定问题(一一映射)
    1. 解存在
    2. 解唯一(一个输入对应且仅对应一个输出)
    3. 输出随输入连续改变(一个输出对应且仅对应一个输入)
  2. 不适定问题(上述违反其一)
  3. 病态问题 ——输出结果对输入数据非常敏感的问题,如果输入数据有微小误差,则引起问题解的相对误差很大
    1. 病态方程组的鉴定及解决
核心问题

条件数

  • 性质

迭代法性质

  1. 收敛性(不动点迭代)
  2. 历史性(单步迭代)
  3. (全局)收敛性基本定理 条件充足时的“强定理”

    前提:闭区间连续,开区间可导 映内性:迭代函数值不会跳出定义域 压缩性:迭代函数的导函数的绝对值不大于一个小于1的常数(每一次映射都是一次聚焦的过程,且这个聚焦的倍率还有更严格的限制)

    1. 唯一性:区间内不动点唯一
    2. 任意性:区间内初值任意,不影响迭代序列收敛
    3. 误差不等式:$|x_k-x^*|\leq\frac{L}{1-L}|x_k-x_{k-1}|$ $ps.$可推广至迭代初值
  4. 局部收敛性定理 条件不足时的“弱定理”

    迭代函数的导函数的绝对值在不动点邻域内连续,且小于1即可

    1. 该迭代函数局部收敛
  5. 收敛阶 ——误差递减的速率

类比最优化方法的收敛速度

  • 线性收敛
  • 超线性收敛
    • 二次收敛(平方收敛)
    • 1.618阶收敛
  • 整数阶超线性收敛定理(局部收敛的推广)
  1. 收敛迭代/最优化方法的终止准则
  2. 迭代法收敛性基本定理(方程组/矩阵)
  • 充要条件:迭代矩阵的谱半径$\rho(B)<1$

系数矩阵 -> 迭代矩阵 -> 收敛性(与初始向量和方程组右端项无关)

雅可比迭代法

迭代矩阵的构造路径
  1. 雅可比迭代法的的矩阵形式(掌握推导) $$Ax=b$$$$A=D-L-U$$$$(D-L-U)x=b$$$$Dx=(L+U)x+b$$$$x^{(k+1)}=D^{-1}(L+U)x+D^{-1}b=D^{-1}(D-A)x+D^{-1}b$$$$=(I-D^{-1}A)x+D^{-1}b$$$$=B_Jx^{(k)}+f_J$$ $$B_J=(I-D^{-1}A)=D^{(-1)}(L+U)$$$$f_J=D^{-1}b$$

  2. 雅可比迭代法的分量形式(掌握计算)

import numpy as np
# 主函数
def Jacobi(A, B, x0, x, eps, k):
    """雅可比迭代法
        Args:
            A:为方程组的系数矩阵
            B:为方程组右端的列向量
            x:为迭代初值构成的列向量
            x:迭代向量
            eps :精度误差
            k:最大迭代次数
        Returns:
            times:迭代次数
            x:Array数组,方程的解
        Raises:
    """
    n, m = A.shape
    if n != m:
        print("invalid data!")
        return None
    times = 0
    while times < k:
        for i in range(n):
            temp = 0
            for j in range(n):
                if i != j:
                    temp += x0[j] * A[i][j]
            x[i] = ((B[i] - temp) / A[i][i])
            print(x[i])
        error = max(abs(x - x0))
        times += 1
        if error < eps:
            print("精确度等于", eps, "时,需要迭代", times, "次收敛")
            return (x, times)
        else:
            x0 = x.copy()
    print("在最大迭代次数内不收敛", "最大迭代次数后的结果为", x)
    return None
 
if __name__ == '__main__':
    a = np.array([[5,2,1],[-1,4,2],[2,-3,10]])
    b = np.array([-12,20,3])
    x = np.array([0.0, 0.0, 0.0])   # 迭代初始值
    x0 = x.copy()
    g = 1e-4
    Jacobi(a, b, x0, x, g, 100)
    print(x)


# 利用Python实现雅可比迭代法的核心算法,求解线性方程组,并利用一个实例对结果进行验证
# import numpy as np
# # x_{k+1} = x_k + A^{-1}b - A^{-1}Ax_k = A^{-1}b
# A = np.array([[2,1,0,0],[1,2,-3,0],[0,3,-7,4],[0,0,2,5]])
# b = np.array([3,-3,-10,2])
# def F(X):
#     return A.dot(X) - b

# J = A
# x0 = np.zeros(4)
# x = x0 - np.linalg.inv(J).dot(F(x0))
# # linalg是一种内置的封装好的方程解决方案
# print(x)

高斯-塞德尔迭代法

迭代矩阵的构造路径
  1. 高斯-塞德尔迭代法的矩阵形式(掌握推导) $$Ax=b$$$$A=D-L-U$$$$(D-L-U)x=b$$$$Dx=Lx+Ux+b$$$$Dx^{(k+1)}=Lx^{(k+1)}+Ux^{(k)}+b$$$$(D-L)x^{(k+1)}=Ux^{(k)}+b$$$$x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b$$ $$B_{G-S}=(D-L)^{-1}U$$$$f_{G-S}=(D-L)^{-1}b$$

  2. 高斯-塞德尔迭代法的分量形式(掌握计算)

import numpy as np
 
def Jacobi(a,b,x0,x,eps,k):#k迭代次数
    m,n=a.shape
    if(m!=n):
        print("NO answer")
    times=0
    while times<k:
        for i in range(n):
            sum=0
            for j in range(n):
                if(i!=j):
                  sum=sum+a[i][j]*x[j]
            x[i]=(b[i]-sum)/a[i][i]
            x[i]=x[i].copy() #每算出一个x[i],就更新一个X[i],
        print(x)
        error = max(abs(x - x0))
        times += 1
        if error < eps:
            print("精确度等于", eps, "时,需要迭代", times, "次收敛")
            return (x, times)
        else:
            x0 = x.copy()
    print("在最大迭代次数内不收敛", "最大迭代次数后的结果为", x)
    return None
if __name__ == '__main__':
    a=np.array([[5.0,2.0,1.0],[-1.0,4.0,2.0],[2.0,-3.0,10.0]])
    b = np.array([-12.0,20.0,3.0])
    x = np.array([0.0, 0.0, 0.0])   # 迭代初始值
    x0 = x.copy()
    g = 1e-4
    Jacobi(a,b,x0,x,g,100)
  • 是否考察证明(待定)
  • 两种方法的关联:当雅可比迭代矩阵为非负矩阵(全体元素均非负数的矩阵)时,雅可比方法与高斯-塞德尔方法同收敛,同发散;若同收敛,则高斯-塞德尔方法收敛更快。

曲线拟合方法

引入
插值基函数的由来
  1. 原则上,按照旧有的方法,对于有$n+1$个样本点的原始数据集,我们需要求解一个 $n+1$$n$ 次方程组,才能得到这条多项式拟合曲线。
  2. 灵感 从众多簇base function中进行整流的激活函数 (详见如何直观地理解拉格朗日插值法
  3. 构造 严谨表达:利用乘法中0因子的湮灭效应
评价
  1. 插值多项式存在且唯一
  2. 余项定理对所有插值多项式均恒成立

最佳逼近理论

概念引入
再说一遍,我不是插值!

  插值和函数逼近有很大的相似性,但是插值得到的函数是要过所有已知数据点的,而函数逼近得到的函数与已知的函数/数据点是存在一定误差的。两者没有绝对的优劣,根据实际情况自行选用即可。(需要过已知数据点的,就用插值,不需要的,就用函数逼近。)

偏差
  • $$||f-p||\infty=\max{a\leq x\leq b}|p(x)-f(x)|=\mu$$

类比函数的范数定义

正式定义
  • 使得偏差最小的多项式即为最佳一致逼近多项式
求解最佳一致逼近多项式
  • 线性 STEP 1:区间定域 $$a_1=\frac{f(b)-f(a)}{b-a}$$ STEP 2:反向求导 $$f'(x_2)=a_1$$ STEP 3:代入通项 $$a_0=\frac{1}{2}[f(a)-f(x_2)]-\frac{a+x_2}{2}\ a_1$$

  • 高次 STEP 1:区间变换 通过换元法将新变元区间映射到$[-1,1]$上 STEP 2:整体代入并化简求得标准表达式 STEP 3:依次数决定切比雪夫多项式并作差比较 STEP 4:回代过程

求解最佳平方逼近多项式

希尔伯特矩阵(高度病态)和勒让德多项式是一对好兄弟

KEY POINTS

最优化问题

凸集/凸函数/凸优化

概念考察篇
  • 凸集:点集上任意两点间连线上的点仍属于该点集的集合
  • 凸函数:函数图像上任意两点间连线上的点在该函数图像的上方的函数本身
  • 严格凸函数
  • 可行域
  • 整体最优解
  • 严格整体最优解
  • 局部最优解
  • 严格局部最优解
  • 局部极小点
  • 严格局部极小点
  • 全局极小点
  • 严格全局极小点

以下类比高数

  • 一阶必要条件
  • 二阶必要条件
  • 二阶充分条件
  • 凸最优性定理

INSTANCES

二分法

一般步骤

STEP 1:捕获一个满足条件的初始“隔离区间”

  • 区间足够小(以至于函数图象单调)
  • 区间足够一般(以至于不包含间断点)
  • 区间足够特殊(以至于函数图象恰好穿轴)

STEP 2:进行一波递归操作(过程中时不时比较精度)

STEP 3:精度达标,停止操作,读取结果

定参

$$二分次数:k(from\ 0)$$ $$区间长度:\frac{b-a}{2^k}$$ $$区间中点:x^k=\frac{a_k+b_k}{2}$$ $$理论解:x^*$$

KEY POINTS
  • 误差应严格小于精度(偏离阈)
  • 用所得结果$k>?$进行取整约束)
  • 迭代步数具有初值依赖性
  • 收缩速度:与等比级数同阶(线性收敛)

牛顿迭代法(切线法)

一般步骤

STEP 1:捕获一个满足条件的初始区间

  • 区间足够小(以至于函数的二阶导不变号)
  • 区间足够一般(以至于函数“增长”不停滞)(一阶导作分母)
  • 区间足够特殊(以至于函数图象恰好穿轴)
  • 初值足够特殊(以至于初始函数值与其二阶导同号)

STEP 2:始于迭代初值的一系列一阶泰勒展开(切线方程)

STEP 3:

import sympy as sp

# 牛顿迭代函数
def NTM(f, p0, tol, maxK): # 入四出三
    x = sp.symbols('x')
    P = [p0] # 当前迭代初始值
    k = 2 # 当前迭代次数
    df = sp.diff(f, x) # 求f关于x的导数

    while k <= maxK:
        P.append(P[k - 2] - f.subs(x, P[k - 2]) / df.subs(x, P[k - 2]))

        err = abs(P[k - 1] - P[k - 2]) # 计算绝对误差
        if err < tol:
            print(f'迭代{k - 1}次即可满足允许误差值退出')
            break

        k += 1

    if k - 1 == maxK:
        print("超过最大迭代次数!")

    p = P[k - 2]
    k -= 1
    Y = P

    return p, k, Y

# 主程序
x = sp.symbols('x')
f = input("请输入需求零解的方程f(x)=(自变量为x,如x^3-x^2-5):   ")
f = sp.sympify(f)
p0 = float(input("请输入牛顿迭代法的初始值p_0:  "))
tol = float(input("请输入精度E:  "))
maxK = int(input("请输入最大迭代次数:  "))

p, k, Y = NTM(f, p0, tol, maxK)

DP = "使用牛顿迭代法迭代{}次,计算{}=0以{}为迭代初始值的解为:{}".format(k, f, p0, p)
print(DP)
print("迭代值如下:")
print(Y)
定参

$$迭代次数:k(from\ 0)$$ $$迭代公式:x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$$ $$理论解:x^*$$

KEY POINTS
  • 是否考察证明(待定)
  • 收敛速度:具有“平方”收敛速度(至少平方局部收敛)证明(待定)
  • 二重近似:泰勒公式+有限次迭代
  • 能求重根
  • 对函数性态有一定要求(结构简单且可导)
  • 牛顿迭代法属于不动点迭代法吗?(待定)

弦截法(离散牛顿法)

一般步骤

STEP 1:捕获一个满足条件的初始“隔离区间”

  • 初始要求同 二分法

STEP 2:始于迭代初值的一系列割线方程

  • 确保两次相邻迭代位点同号

STEP 3:

import sympy as sp

# 弦截函数
def SM(f, p0, p1, tol, maxK): # 入五出三
    x = sp.symbols('x') # 创建一个符号变量x,用于表示方程中的自变量
    P = [p0, p1] # 创建一个列表,用于存储弦截法的迭代值
    k = 2 # 初始化一个迭代计数器k为2,因为已经有两个初始值

    while k <= (maxK + 1):
        df = (f.subs(x, P[k - 1]) - f.subs(x, P[k - 2])) / (P[k - 1] - P[k - 2])
        P.append(P[k - 1] - f.subs(x, P[k - 1]) / df)

        err = abs(P[k] - P[k - 1])
        if err < tol:
            print(f'迭代{k - 1}次即可满足允许误差值退出')
            break

        k += 1

    if k - 1 == maxK:
        print("超过最大迭代次数!")

    p = P[k - 1]
    k -= 1
    Y = P

    return p, k, Y

# 主程序
x = sp.symbols('x')
f = input("请输入求零解的方程f(x)=(自变量为x,如x^3-x^2-5):   ")
f = sp.sympify(f) # 将用户输入的方程字符串转换为SymPy表达式
p0 = float(input("请输入弦截法的初始值p_0:  "))
p1 = float(input("请输入弦截法的初始值p_1:  "))
tol = float(input("请输入精度E:  "))
maxK = int(input("请输入最大迭代次数:  "))

p, k, Y = SM(f, p0, p1, tol, maxK)

DP = "使用弦截法迭代{}次,计算{}=0 以x0={}, x1={}为迭代初始值的解为:{}".format(k, f, p0, p1, p)
print(DP)
print("迭代值如下:")
print(Y)
定参

$$迭代次数:k(from\ 1)$$ $$迭代公式:x_{k+1}=x_k-\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}f(x_k)$$ $$ps.单点弦截法的情形下x_{k-1}项退化为x_0$$ $$理论解:x^*$$

KEY POINTS
  • 收敛速度:单点弦截法具有“线性”收敛速度;双点弦截法具有“1.618阶”收敛速度
  • 初值依赖性同 二分法 (详见收敛速度的三种形式
  • 不属于不动点迭代

广义迭代法

一般步骤

STEP 1:捕获一个满足条件的初始“隔离区间”

  • 初始要求同 二分法
  • 原方程能够转化为独立出方程的根x的形式(可能存在多种形式)
  • 选择上述形式中满足收敛条件的迭代序列

STEP 2:始于迭代初值的一波递归操作(过程中时不时比较精度)

STEP 3:精度达标,停止操作,读取结果

定参

$$迭代次数:k(from\ 0)$$ $$迭代公式:x_{k+1}=\phi(x_k)$$ $$理论解:x^*$$

KEY POINTS
  • 选取迭代函数$\phi(x)$ ,使迭代公式$x_{k+1}=\phi(x_k)$收敛

拉格朗日插值(多项式)

一般步骤

STEP 1:多项式的求解 $$L_n(x)=\sum_{i=0}^n(\prod_{j=0,\ j\neq i}^n\frac{x-x_j}{x_i-x_j})y_i$$

STEP 2:代数化简(特殊化)

STEP 3:摁计算器

STEP 4:误差估计(通项与某点)

  • 算出$M_{n+1}$(原函数的 $n+1$ 阶导在定义域下的最大值)(前提:高阶导存在)
  • 代入余项定理并建立误差不等式(乘以一个插值多项式的分子部分,再除以一个 $(n+1)!$
  • 若要求某点的误差,则代数计算
  • 记住带上绝对值与精度要求进行比较
函数部分
function yy = Lagrange(x,y,xx)
n = length(x);
% 有 n 个数据点,就有 n 个基函数
% 下面计算基函数 l。基函数计算公式如下

lx = sym(ones(1,n));
syms X
for i = 1:n
    for j = 1:n
        if i==j
            continue
        else
            lx(i)=lx(i)*(X-x(j))/(x(i)-x(j));
        end
    end
end
% 利用得到的基函数,计算目标点函数值,计算公式如下。


Lx = sum(y.*lx);
Lx = matlabFunction(Lx);
yy = Lx(xx);
end
定参

$$互异点的个数:n+1$$ $$插值多项式的次数/插值函数的阶数:n$$ $$插值基函数:l_i(x)\ (i=0,1,...,n)$$ $$插值多项式:L_n(x)$$ $$M_{n+1}=\max_{a\leq x\leq b}|f^{(n+1)}(x)|$$ $$插值余项(误差):R_n(x)$$

KEY POINTS
  • 当阶数为1时,拉格朗日插值退化为线性插值
  • 当阶数为2时,拉格朗日插值退化为抛物插值
  • 当原函数为不超过n次的多项式时,插值多项式就是原函数
  • 插值点的选取:尽量靠近(误差)所求点(用于规避冗余条件)
  • 事后误差估计(仅供参考的实用技巧) (详见整活用
  • 高次插值 = 精度高 +(数值不稳定/不收敛/龙格现象/过拟合)
  • 全部节点依赖性(牵一发而动全身)-> 适合求固定节点的函数值

拉格朗日插值(分段线性)

一般步骤

STEP 1:捕获一个所在区间满足如下条件的函数

  • 区间足够小(使得该函数是每个小区间上的线性函数)
  • 区间足够一般(使得该函数在该区间内不存在间断点)
  • 区间足够特殊(使得每个小区间的左端点的函数值恰好等于原函数值)
定参

$$h=\max_{0\leq i\leq n-1}|x_{i+1}-x_i|$$ $$|R(x)|\leq \frac{h^2}{8}M$$

KEY POINTS
  • 相邻近点相关性,远则无关(信息缺失)
  • 其他类比拉格朗日插值多项式
  • 整个区间充分接近被插函数,但光滑性差

牛顿插值(多项式)

一般步骤

STEP 1:列差商(均差)表 STEP 2:构造牛顿插值公式 STEP 3:代数化简(特殊化)

定参
KEY POINTS
  • 核心思想:通过提前计算出一系列修正项,来规避拉氏插值多项式的重新计算问题所带来的负担
  • 根据精度要求逐步增加节点且计算量小

埃尔米特插值(多项式)

一般步骤(基函数法)

STEP 1:判定前提条件

  • 一阶导数存在 STEP 2:将4个基函数整合成$H(x)$ STEP 3:代数化简(特殊化)
二般步骤(两点三次)(待定系数法)

STEP 1:构造三次多项式形式 STEP 2:对上式求导 STEP 3:指定系数,求解方程组

定参

$$基函数:h_0(x),h_1(x),H_0(x),H_1(x)$$ $$h_0(x)=(1+2\frac{x-x_0}{x_1-x_0})(\frac{x-x_1}{x_0-x_1})^2$$ $$h_1(x)=(1+2\frac{x-x_1}{x_0-x_1})(\frac{x-x_0}{x_1-x_0})^2$$ $$H_0(x)=(x-x_0)(\frac{x-x_1}{x_0-x_1})^2$$ $$H_1(x)=(x-x_1)(\frac{x-x_0}{x_1-x_0})^2$$

KEY POINTS
  • 公式记忆可采取对称性互换原理
  • 补充:余项定理$R(x)=\frac{f^{(4)}(\epsilon)}{4!}(x-x_0)^2(x-x_1)^2$
  • 适合于导数值已知的情形

样条函数插值(多项式)

一般步骤(基函数法)

STEP 1:区间满足以下条件

  • 二阶导数连续
  • 三次多项式(三次样条函数)
  • 插值条件(使得每个小区间的左端点的插值函数值恰好等于原函数值)

STEP 2:将4个基函数整合成$H(x)$ STEP 3:代数化简(特殊化)

二般步骤(待定系数法)

STEP 1:构造$(n-1)$个三次多项式 STEP 2:多次求导,多次代入插值条件 STEP 3:求解方程组

定参

$$S''(x)=M_j\frac{x_{j+1}-x}{h_j}+M_{j+1}\frac{x-x_j}{h_j},\quad x\in[x_j,x_{j+1}],j=0,1,\cdots,n-1$$ $$\begin{aligned} S(x)&=M_j\frac{(x_{j+1}-x)^3}{6h_j}+M_{j+1}\frac{(x-x_j)^3}{6h_j}+\left(y_j-\frac{M_jh_J^2}{6}\right)\frac{x_{j+1}-x}{h_j}+\left(y_{j+1}-\frac{M_{j+1}h_j^2}{6}\right)\frac{x-x_j}{h_j},\quad j=0,1,\cdots,n-1 \end{aligned}$$ $$第一边界条件:节点一阶导相等$$ $$第二边界条件:节点二阶导相等$$ $$第三边界条件:周期性条件$$ $$\lambda$$ $$\mu$$ $$d$$

KEY POINTS
  1. 样条插值要求二阶导数连续,而hermite插值只需要一阶导数存在。
  2. 样条插值用到三个点,而hermite是两个。
  3. 样条函数计算一次得出两个三次多项式,而hermite则是一个。

黄金分割法

一般步骤

STEP 1:确定黄金分割法的终止准则

  • 搜索区间的长度小于预定的容限为止
  • 函数值差值:构造梯度迭代凸函数

STEP 2: 根据区间边界设定两个对称点(吸引子) STEP 3:比较大小,逐次迭代

clc;
clear;
a = -1;b = 1; %初始区间
f = @(x) 2 * x ^ 2 -  x - 1; %创建题目要求匿名函数,方便使用
eps = 0.08; %区间精度
while((b - a) >= eps)
    x1 = a + 0.382 * (b - a);
    x2 = a + 0.618 * (b - a); %黄金分割法主要步骤
    if f(x1) < f(x2) %两种情形的判断
        b = x2;
    else
        a = x1;
    end
end
x = (a + b) / 2; %得到满足条件的最优解
disp(['最优解: x = ',num2str(x)]);
disp(['此时: f(x) = ',num2str(f(x))]);%使用disp函数和num2str()进行输出
定参
KEY POINTS
  • 一维搜索算法
  • 适用于单峰函数搜索极值
  • 收敛速度为0.618

二分法

一般步骤

STEP 1:确定二分法的初始搜索区间与终止区间长度

  • 迭代过程中始终满足区间左右端点导数值异号

STEP 2:区间取中点,并求对应函数的导数值,分类讨论

  • 向左归并
  • 向右归并

STEP 3:建立新的搜索区间,直到满足终止条件

function solve = dichotomy(f,a,b,delta0)
n = 100;
x = (a+b)/2;
for j = 1:n
    if f(a)*f(x) < 0
        b = x;
    elseif f(a)*f(x) > 0
        a = x;
    end
    x = (b+a)/2;
    if (b-a) < delta0
        break
    end
end
solve = x;
end
定参

$$初始区间:[a_1,b_1]$$ $$终止区间长度上限:\delta$$ $$迭代次数n满足:(\frac{1}{2})^n\leq \frac{\delta}{b_1-a_1}$$

KEY POINTS