《工程数值计算Python教程》笔记

文章目录

• • @[toc]
  • 第一章：绪论
  • • $1.1$|数值计算在工程科学中的重要性
    • $1.2$|数值计算方法
    • $1.3$|程序设计
    • • 盒图
      • 计算方法的选取
      • • 减少运算次数
        • 避免相近的数相减
    • $1.4$|误差的来源、表示及传递
    • • 误差的来源和分类
      • • 模型误差
        • 观测误差
        • 截断误差
        • 舍入误差
      • 误差的表示
      • • 绝对误差
        • 相对误差
        • 平均误差
        • 标准误差
      • 误差的传递
      • • 误差在和、差计算中的传递
        • • 绝对误差
          • 相对误差
        • 误差在积、商计算中的传递
        • • 乘积的绝对误差
          • 乘积的相对误差
          • 商的绝对误差
          • 商的相对误差
  • 第二章：`Python`基础
  • • $2.1$|概述
  • 第三章：方程（组）的求解
  • • $3.1$|非线性代数方程的求根
    • • 二分法
      • • 原理
        • 步骤
        • `Python`实现

第一章：绪论

$1.1$|数值计算在工程科学中的重要性

$1.2$|数值计算方法

$1.3$|程序设计

盒图

• 盒图也称为$N-S$结构化流程图
• 如果一个算法可以用$N-S$流程图描述，则说明它是结构化的
• 盒图示例

计算方法的选取

减少运算次数

• 计算多项式$P_n(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$的值
  • 如果逐项计算，所需乘法次数为$1+2+\cdots +n=\frac{1}{2}n(n+1)$，所需加法次数为$n$

$$\begin{gathered} S_0=a_0 \\ {S}_k=a_k{x}^k(k=1,2,\cdots ,n) \\ {P}_n(x)=\sum _{i=0}^n{S}_i \end{gathered}$$

• • 如果采用递推法（秦九韶算法），则只需$n$次乘法和$n$次加法

$$\begin{gathered} S_n=a_n \\ {S}_k=x{S}_{k-1}+a_k(k=n-1,n-2,\cdots ,0) \\ {P}_n(x)=S_0 \end{gathered}$$

避免相近的数相减

• 在计算过程中，应避免相近的数相减，不然会使有效数字的位数大大减少

$1.4$|误差的来源、表示及传递

误差的来源和分类

模型误差

• 模型与实际问题的误差

观测误差

• 通常数学模型中都包含一些需要实验测定的参数，这些参数存在测定值与真实值之间的误差

截断误差

• 数值计算方法中，常用收敛无穷级数的前几项代替无穷级数
• 截断误差与算法有关，常用截断误差限或截断误差的阶来判断某种算法的优劣

舍入误差

• 由于计算机位数有限引起的误差称为舍入误差
• 舍入误差与计算机有关，也与数学表达式有关

误差的表示

绝对误差

• 用$x$表示准确值$x^{\ast }$的近似值，绝对误差为$e=x-x^{\ast }$
• 通常$x^{\ast }$是未知的，故$e$的真实值也是不知道的，常根据实际情况估算它的上限，$|e|=|x-x^{\ast }|\le \epsilon$，$x-\epsilon \le x^{\ast }\le x+\epsilon$，$\epsilon$称为$x$的绝对误差限

相对误差

• $e_r=\frac{e}{x}=\frac{x-x^{\ast }}{x}$
• $|e_r|=|\frac{x-x^{\ast }}{x}|\le {\epsilon }_r$

平均误差

• 算数平均值$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$
• 平均误差为$\delta =\frac{1}{n}\sum _{i=1}^n{d}_i=\frac{1}{n}\sum _{i=1}^n|\bar{x}-x_i|$

标准误差

• $\sigma =\sqrt{\frac{1}{n-1}\sum _{i=1}^n{d}_i^2}$
• 标准误差也称为均方根误差

误差的传递

误差在和、差计算中的传递

• 设$x^{\ast }$、$y^{\ast }$的近似值为$x$、$y$，$z^{\ast }=x^{\ast }+y^{\ast }$

绝对误差

• $e_z=e_x+e_y$
• $|e_z|\le |e_x|+|e_y|$

相对误差

• $e_{rz}=\frac{(x+y)-(x^{\ast }+y^{\ast })}{x+y}=\frac{x-x^{\ast }}{x}\times \frac{x}{x+y}+\frac{y-y^{\ast }}{y}\times \frac{y}{x+y}$
  • 当$x$与$y$同号时，$|\frac{x}{x+y}|\le 1$，$|\frac{y}{x+y}|\le 1$，$e_{rz}\le |e_{rx}|+|e_{ry}|$
  • 当$x+y\approx 0$时，可能有$e_{rz}\gg |e_{rx}|+|e_{ry}|$

误差在积、商计算中的传递

• 设$x^{\ast }$，$y^{\ast }$均为正数，近似值分别为$x$、$y$，绝对误差为$dx=x-x^{\ast }$，$dy=y-y^{\ast }$，相对误差为$\frac{x-x^{\ast }}{x}=\frac{dx}{x}=d\ln x$，$\frac{y-y^{\ast }}{y}=\frac{dy}{y}=d\ln y$

乘积的绝对误差

• $d(xy)=xdy+ydx$

乘积的相对误差

• $d\ln xy=d(\ln x+\ln y)=\frac{dx}{x}+\frac{dy}{y}$

商的绝对误差

• $d(x/y)=\frac{ydx-xdy}{y^2}$

商的相对误差

• $d\ln x/y=d(\ln x-\ln y)=\frac{dx}{x}-\frac{dy}{y}$

第二章：Python基础

$2.1$|概述

第三章：方程（组）的求解

$3.1$|非线性代数方程的求根

二分法

原理

• 设$f$是一个连续函数，假设存在一个区间$[a,b]$，$f$在区间的两端点处的值符号相反，即$f(a)f(b)<0$，于是，$f$在区间$[a,b]$有一个根，这可由函数的中值定理得到

步骤

• 有一个区间$[a,b]$和值$u=f(a)$、$v=f(b)$，满足$uv<0$

• 构造区间的中点$c=(a+b)/2$并计算$w=f(c)$，可能碰巧$w=0$，则算法结束

• 通常情况下$w\ne 0$，如果$wu<0$，则可断定在区间$[a,c]$中存在$f$的一个根，因而将$c$的值存储到$b$中，将$w$的值存储到$v$中，如果$wv<0$，则$f$在区间$[c,b]$中有一个根，因而将$c$的值存储到$a$中，将$w$的值存储到$u$中

• 重复进行，每次区间长度减少一半，直到区间达到需要的精度

• 结束时，根的最佳估计是$(a+b)/2$

• 在程序运行过程中，中点处函数值$w=0$的概率是极低的，为此增加一个判断语句是不值得的，对任意区间$[a,b]$，只要$uv\le 0$即可保证区间中包含零点

Python实现

• 求解$f(x)=x^3-3x+1$在区间$[0,1]$上的零点

def bisect(f, a, b, eps=1e-9, args=()):
    u, v = f(a, *args), f(b, *args)

    if u * v <= 0:
        while abs(a - b) > eps:
            c = (a + b) / 2
            w = f(c, *args)

            if w * u <= 0:
                b, v = c, w
            else:
                a, u = c, w

        return (a + b) / 2
    else:
        print(f'The function values at {
     a} and {
     b} have the same sign')


def f(x): return x * (x * x - 3) + 1


if __name__ == '__main__':
    a, b = 0, 1
    eps = 1e-9

    res = bisect(f, a, b, eps)

    print('The result is:', res)

• 求解$f(x)=\cosh (\sqrt{x^2+1}-e^x)+\log |\sin x|$在区间$[0,1]$上的零点
  • 该式在$0$点没有定义，计算会发生溢出，此时可以选择一个很小的值作为左边界，如$10^{-10}$

import numpy as np
from bisection import bisect


def f(x): return np.cosh(np.sqrt(x * x + 1) - np.exp(x)) + np.log10(np.abs(np.sin(x)))


if __name__ == '__main__':
    a, b = 1e-10, 1

    res = bisect(f, a, b)

    print(f'The result is: {
     res}')