数值计算——浮点数，误差，范数，并行计算及python实现

一、计算机中数的浮点表示

1、数的浮点表示

允许小数点位置浮动的表示方法称为数的浮点表示，一个数有多种浮点表示。

• 规格化的浮点数表示为
  $$x=\pm 0.a_1{a}_2...a_s\times {\beta }^c,$$
  其中，$$1\le a_1<\beta ,0\le a_i<\beta ,i=2,3,\cdots ,s.$$称$0.a_1{a}_2...a_s$为浮点数的尾数部分，$s\ge 1$，大小不限制，它为$\beta$进制整数

---

2、计算机中的浮点数

计算机均采用二进制存储数据，IEEE 浮点数标准用三元组 { S , E , M } \{S,E,M\} {S,E,M}来表示一个数$V$的，即
$$V=(-1{)}^S\times M\times 2^E,$$
其中：

• $S$(Sign)为符号位
• $E$(Exponent)为指数位
• $M$(Significant)为有效数字位

---

3、浮点数的计算机表示

计算机位数有限，所以浮点数尾数部分的位数是固定的，如$t$位，也称为计算机字长

• 机器数：受具体计算机限制的浮点数，它们构成机器的数系，记为$F(\beta ,t,m,M)$,$m\le c\le M$，设$f\in F(\beta ,t,m,M)$，有
  $${\beta }^{m-1}\le |f|\le (1-{\beta }^{-t}{\beta }^M).$$
  小于下限，机器自动作为零处理，高于上限则中断运算。
• 若$x\notin F(\beta ,t,m,M)$，机器将自动用机器数$f$表示，记为$f=fl(x)$，计算机会按照一定的舍入或截断规则得到$fl(x)$

---

二、误差的基本概念

1、误差来源

模型误差：数学模型往往是抓住主要因素后进行数学概括，因此数学模型是近似的，存在模型误差。

• 观测误差：观测和实验中产生的误差。
• 截断误差：模型的精确值与数值方法精确解之差。
• 传播误差：初始误差在计算在计算中不断积累最终对计算结果产生影响。
• 舍入误差：计算机浮点数的字长有限，计算过程中按四舍五入或只舍不入截成有限位数的机器数，所引起的误差。

---

2、近似数的误差和有效数字

设数 $x$是某个量的精确值, 数 $x^{\ast }$是该量的已知近似值, 记
$$E(x)=x-x^{\ast }$$
称 $E(x)$为近似数$x^{\ast }$的绝对误差, 简称误差.

• 近似数$x^{\ast }$的绝对误差和精确值之比，即
  $$E_{\Gamma }(x)=\frac{x-x^{\ast }}{x}$$
  称为$x^{\ast }$的相对误差。

• 设精确值$x$的近似值为
  $$x^{\ast }=\pm 0.a_1{a}_2\cdots a_k\cdots a_n\times 10^m$$
  其中$a_i(i=1,2,\cdots ,n)$是$0$到$9$中的某一整数，并且$a_1\ne 0$。
  如果$x^{\ast }$的绝对误差限满足
  $$|x-x^{\ast }|\le \frac{1}{2}\times 10^{m-k}$$
  则称$a_k$为有效数字。
  若$a_k$为有效数字，那么$a_1,\cdots ,a_{k-1}$都是有效数字。
  若$$|x-x^{\ast }|\le \frac{1}{2}\times 10^{m-n}$$，那么称$x^{\ast }$称为具有$n$位有效数字的有效数。

• 按“四舍五入”规则来选取近似数
  $$x=\pi =10\times 0.314159266\cdots ,$$
  三位有效数字：$x_3^{\ast }=10\times 0.314$
  五位有效数字：$x_5^{\ast }=10\times 0.31416$

# 计算pi的有效数字
def significant_figures(x_appr,x_prec):
    p = 0 # the order of the approximated number
    n = 0
    # 以下把x_appr化成规格化浮点数，计算阶数 p9
    x = x_appr
    while x >=1:
        x /= 10
        p += 1
    while x < 0.1:
        x *= 10
        p -= 1
 
    d = abs(x_appr-x_prec)
    while True:
        n += 1
        if d > 0.5 * 10**(p-n):
            break
  
    return n-1
print("3.14相对于pi有效数字为",significant_figures(3.14,3.1415926)) 
print("3.1415相对于pi有效数字为",significant_figures(3.14159,3.1415926))     

$x$的近似数$x^{\ast }$的每一位不一定都是有效数字，绝对误差要满足一定条件，不能超过末尾数字的半个单位

• 有效数字越多，相对误差限越小

---

3、初始误差在运算中传播

不考虑计算是按照什么方式，不考虑计算过程中舍入误差的积累,计算产生的误差用泰勒展开来估计

• 加减
  $$f(x_1,x_2)=x_1+x_2$$
  $$E(x_1+x_2)\approx \frac{\partial (x_1+x_2)}{\partial x_1}E(x_1)+\frac{\partial (x_1+x_2)}{\partial x_2}E(x_2)=E(x_1)+E(x_2)$$
  $$E_{\Gamma }(x_1+x_2)\approx \frac{x_1{E}_{\Gamma }(x_1)}{x_1+x_2}+\frac{x_2{E}_{\Gamma }(x_2)}{x_1+x_2}$$
  可以得到两近似数作加法运算时计算结果的绝对误差限和相对误差限均不超过各项绝对误差限和相对误差限之和
  避免两个接近的数相减时引起精度下降：

$x_1$很接近$x_2$ ：$lg{x}_1-la{x}_2=lg\frac{x_1}{x_2}$
$x$接近于0时变换：$\frac{1-cosx}{sinx}=\frac{sinx}{1+cos}$
$x$充分大时变换：$arctan(x+1)-arctanx=arctan\frac{1}{1+x(x+1)}$

第一种计算方式为什么可以解决精度下降的问题？
带入上面的公式计算误差

• 乘除
  $$f(x_1,x_2)=x_1\cdot x_2$$
  $$E(x_1\cdot x_2)\approx \frac{\partial (x_1\cdot x_2)}{\partial x_1}E(x_1)+\frac{\partial (x_1\cdot x_2)}{\partial x_2}E(x_2)=x_{2}E(x_1)+x_{1}E(x_2)$$
  $$E_{\Gamma }(x_1\cdot x_2)\approx \frac{1}{x_1\cdot x_2}[x_{2}E(x_1)+x_{1}E(x_2)]=E_{\Gamma }(x_1)+E_{\Gamma }(x_2)$$
  当 $x_1$或$x_2$ 的绝对值很大时$|E(x_1\cdot x_2)|$的值可能很大, 即导致 $x_1$$乘$x_2$ 的绝对误差放大;
  在实际计算中，避免两个相近数的相减、大的乘数作乘法以及接近于零的除数作除法。

---

4、浮点运算舍入误差

规格化的浮点机器数经过算数运算后的精确结果不一定是规格化的机器数

• 计算机上进行浮点数的四则运算时，总是先进行运算并将数据尾数保留2t位，然后进行舍入或截断成t位浮点机器数
• 每次运算会产生相对少的误差，若$nu\le 0.01$，则
  $$fl(\sum _{j=1}^n{x}_i\cdot x_j)=\sum _{j=1}^n(x_i\cdot x_j)[1+1.01(n+2-j){\theta }_{j}u],|{\theta }_j|\le 1,j=1,2,\cdots ,n.$$

---

5、向前误差分析和向后误差分析

给定一个输入 $x$，我们想通过一个计算过程 $f(x)$获得 $y$， 但是由于计算存在误差，我们实际计算出了 $y^{\ast }$，即
$y^{\ast }\ne f(x)$
此时$y^{\ast }-y$ 称为向前误差或前向误差（forward error)，对应的误差估计方法称为向前误差分析。
假设存在另一个输入 $x^{\ast }$，使得
$$y^{\ast }=f(x^{\ast })$$
成立，则称 $x^{\ast }-x$为实际计算结果 $y^{\ast }$的向后误差或后向误差（backward error)。

两种误差分析要怎么选择？
向前误差分析方法只能应用于十分简单的情形。
向后误差分析方法则不同了,它最大的优点则是把浮点运算的舍入误差估计转化为通常的实数运算,可以在误差分析中运用实数运算法则。

---

三、范数计算

1、范数定义

记$\mathbb{R}$为n维欧式空间，定义${\mathbb{R}}^n$中的$p$范数
$${\left|\left|x\right|\right|}_p=(\sum _{i=1}^n{\left|x_i\right|}^p{)}^{\frac{1}{p}},\forall x=(x_1,x_2,\cdots ,x_n{)}^{\mathrm{T}}\in {\mathbb{R}}^n,p\in [1,\infty ),$$
则${\mathbb{R}}^n$为一赋范线性空间

• 向量的1-范数：各个元素的绝对值之和
  $$||X|{|}_1=\sum _{i=1}^n|x_i|$$
• 向量的2-范数：各个元素的平方和再开方
  $$||X|{|}_2=(\sum _{i=1}^n{x}_i^2{)}^{\frac{1}{2}}=\sqrt{\sum _{i=1}^n{x}_i^2}$$
• 向量的无穷范数：取绝对值最大的元素
  $$||X|{|}_{\infty }=\underset{1\le i\le n}{\max }=|x_i|$$
  用Python实现

import numpy as np
import math
arr = np.array([1,2,3,4])

#计算L1范数
L1=0
for i in arr:
    L1=L1+abs(i)
print("数组的L1范数是",L1)
 
#计算L2范数
L2=0
for i in arr:
    L2=L2+i^2
L2 = math.sqrt(L2)
print("数组的L2范数是",L2)

#计算无穷范数
L=abs(arr[0])
for i in arr:
    if abs(i)>L:
        L = abs(i)
print("数组的无穷范数是",L)
- Chebyshev 范数
$$\Vert f\Vert=\mathop{\max}\limits_{a\le x\le b}|f(x)|,\quad \forall f\in C[a,b]$$。
from scipy import optimize
def f(x):
    return x**3
def Chebyshev(min,max):
    a = abs(minimum)
    b = abs(maximum)
    if a>b:
        print("范数为：",a)
    if a<b:
        print("范数为：",b)
minimum = optimize.fminbound(f, -1, 2)
min = f(minimum)
maximum = optimize.fminbound(lambda x: -f(x), -1, 2)
max = f(maximum)
print (min)
print (max)
Chebyshev(min,max)     

---

2、内积

最简单的内积空间是$n$维欧氏空间${\mathbb{R}}^n$，其中的内积定义为
$$(x,y)=x^{\mathrm{T}}y=\sum _{i=1}^n{x}_i{y}_i$$
$x=(x_1,x_2,\cdots ,x_n{)}^{\mathrm{T}},y=(y_1,y_2,\cdots ,x_n{)}^{\mathrm{T}}$为${\mathbb{R}}^n$中的任意向量

#计算两个向量内积
import numpy as np
vect1 = np.array([1,2,3])
vect2 = np.array([2,3,4])
dotx = np.dot(vect1,vect2)
print("dotx: {}".format(dotx))

• 性质：
  Cauchy-Schwarz（柯西–施瓦茨）不等式
  $$|(x,y)|\le \sqrt{(x,x)\cdot (y,y)},\forall x\in X$$
  平行四边形等式
  $$\Vert x+y{\Vert }^2+\Vert x-y{\Vert }^2=2(\Vert x{\Vert }^2+\Vert y{\Vert }^2)$$
• 内积空间中的向量夹角
  $$\forall x,y\in X,\theta =\arccos \frac{(x,y)}{\Vert x\Vert \cdot \Vert y\Vert }$$
  $\theta$定义为$x,y$的夹角。

#计算向量夹角
import numpy
v3 = np.array([1,0,0])
v4 = np.array([0,1,0])

dotx = np.dot(vect1,vect2)
dot2 = np.dot(v3,v4)
theta = np.arccos(dot2/(np.linalg.norm(v3)*np.linalg.norm(v4)))*180/np.pi

print("theta: {}".format(theta))

• A范数：
  $$||X|{|}_A=\sqrt{(x,x{)}_A},$$
  $A$为任意对称正定矩阵，而$A$内积定义为
  $$(x,y{)}_A=(Ax,y),\forall x,y\in {\mathbb{R}}^n,$$

---

四、几何断言


为什么用这种方式可以判断？

• 判断点与线的关系
  注意两点重合的情况

import numpy as np                                                             
from sympy import *
#点a，b在直线上，判断点c与直线的位置关系
def orientation(a_x,a_y,b_x,b_y,c_x,c_y):
    if a_x != b_x or a_y !=b_y:
  
        M = Matrix([[a_x,a_y,1],[b_x,b_y,1],[c_x,c_y,1]])
        result = M.det()
        #print() 

        if result > 0:
            print("点c在直线左边的区域")
        if result <0:
            print("点c在直线右边的区域")
        if result == 0:
            print("点c在直线上")
    else:   
        print("无法判断")
orientation(2,4,1,7,0,2)

• 判断点是否在圆内
  三点共线确定不了圆的位置

import numpy as np
from sympy import *
 #判断点d与圆的关系
   
def incircle(a_x,a_y,b_x,b_y,c_x,c_y,d_x,d_y):
    if a_x != b_x or a_y != b_y:                                                                                                                                         
        if a_x != c_x or a_y != c_y:
            if b_x != c_x or b_y != c_y:
                M = Matrix([[a_x,a_y,a_x^2+a_y^2,1],
                    [b_x,b_y,b_x^2+b_y^2,1],
                    [c_x,c_y,c_x^2+c_y^2,1],
                    [d_x,d_y,d_x^2+d_y^2,1]])
                result = M.det()
                print (result)
                if result==0:
                    print ("点d在圆上")
                if result > 0:
                     print("点d在圆外")
                if result < 0:
                    print("点d在圆内")
    else:
        print("无法判断")
#点a坐标为(a_x,a_y)，点b坐标为(b_x,b_y)，点c坐标为(c_x,c_y)，点d坐标为(c_x,c_y)
incircle(3,4,2,5,7,2,4,0)
---

五、并行计算

1、并行计算的基本概念

1.1 并行计算的定义

同时利用多个计算资源协同完成某个计算任务的过程称为并行计算。

1.2 进程与线程

进程：操作系统进行资源分配和调度的基本单位，是操作系统结构的基础。

• 线程：处理器调度的基本单位。
• 关系
  一个进程内有至少一个线程，像是线程的“容器”，它被操作系统分配了单独的、与其它进程独立的内存空间，其资源被它的所有线程共享。

---

1.3 多线程和多进程

多线程：用到模块threading，利用CPU 和 IO 可以同时执行的原理，让CPU 不会等待IO 完成。

• 多进程：模块multiprocessing，利用多核CPU的能力，真正并行执行任务。

---

1.4 CPU 密集型与 I/O 密集型

CPU 密集型（CPU-bound）：CPU 密集型也叫做计算密集型，是指I/O在很短时间内就可以完成，CPU 需要大量的计算和处理，特点是CPU 占用率相当高；例如：压缩解压缩、加密解密、正则表达式搜索

• I/O密集型（I/O bound）：IO 密集型指的是系统运作大部分的状况是CPU 在等I/O（硬盘，内存）的读写操作，CPU 占用率较低，例如：文件处理程序，网络爬虫程序，读写数据库程序；

选择多进程还是多线程？
CPU 密集型：用多进程，用多线程反而会运行得慢
I/O 密集型：用多线程

2、多进程计算圆周率 Pi (π) 的值

2.1 用公式求圆周率的近似值

$$\pi =\sum _{n=0}^{\infty }[\frac{1}{16^n}(\frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6})]$$

---

2.2 Concurrent.futures 模块

它提供 ThreadPoolExecutor (线程池) 和 ProcessPoolExecutor (进程池) 两个类，实现了对 threading 和 multiprocessing 的更高级的抽象，对编写 线程池/进程池 提供支持。

concurrent.futures 官方文档：https://docs.python.org/3/library/concurrent.futures.html

2.3 分别用单线程和多进程实现对圆周率的近似

import concurrent.futures as cf                                           
import time
def cal(start,end):
    PI = 0
    for n in range(start, end):
        PI += 1/pow(16, n) * (4/(8*n+1) - 2 / (8*n+4) - 1/(8*n+5) - 1/(8*n    +6))
    return PI
def single_thread(num_iterations):
    start = 0
    end = num_iterations
    cal(start,end)
    
 def multi_process(max_workers,num_iterations):
    PI = 0
    with cf.ProcessPoolExecutor(max_workers) as p:
 
        futures = []
        for i in range(max_workers):
            start = i *num_iterations//max_workers
            end = (i+1)*num_iterations//max_workers
            futures.append(p.submit(cal,start,end))
        for future in cf.as_completed(futures):
            PI += future.result()
        return PI
if __name__ == '__main__':
    max_workers = 3#进程数
    num_iterations = 100000#迭代次数
 
    start_time = time.time()
    single_thread(num_iterations)
    end_time = time.time()
    print("single_thread,cost:",end_time-start_time,"seconds")
 
    start_time=time.time()
    multi_process(max_workers,num_iterations)
    end_time = time.time()
    print(PI)
    print("multi_process,cost:",end_time-start_time,"seconds")