第1章绪论

Contents

第1章绪论#

import numpy as np
import pandas as pd

例子1：计算\(e^{-1} \int_{0}^1 x^n e^n\ dx\)#

def In_almost_exact(n):
    # use quadrature to compute integral
    m = 10**6
    x = np.arange(0, stop=1, step=1/m)
    return np.trapz(np.exp(-1) * x**n * np.exp(x), x=x, dx=1/m)
    
def In_approx(n, r=4):
    # 课本里使用的迭代计算法
    # r是有效数字的位数
    value = np.zeros(n)
    value[0] = np.round((1-np.exp(-1))*10**r)*10**(-r)
    print("I_0 使用近似值 {:.10f}".format(value[0]))
    for j in np.arange(1, stop=n):
        value[j] = 1 - j * value[j-1]
    return value

n = 15
In_exact_value = np.zeros(n) # 可以被认为真实的积分值
for j in range(n):
    In_exact_value[j] = In_almost_exact(j)
    
# 使用4位有效数字
r = 4
approx_vec = In_approx(n, r=r) # 课本使用的近似法
err_vec = np.abs(approx_vec - In_exact_value) # 误差
pd.DataFrame({"n": np.array(np.arange(n)), \
              "真实" : In_exact_value, "近似": approx_vec, "误差": err_vec})

I_0 使用近似值 0.6321000000

	n	真实	近似	误差
0	0	0.632120	6.321000e-01	1.955883e-05
1	1	0.367878	3.679000e-01	2.155883e-05
2	2	0.264240	2.642000e-01	4.011766e-05
3	3	0.207276	2.074000e-01	1.243530e-04
4	4	0.170892	1.704000e-01	4.924119e-04
5	5	0.145532	1.480000e-01	2.468059e-03
6	6	0.126801	1.120000e-01	1.480136e-02
7	7	0.112383	2.160000e-01	1.036175e-01
8	8	0.100931	-7.280000e-01	8.289310e-01
9	9	0.091611	7.552000e+00	7.460389e+00
10	10	0.083876	-7.452000e+01	7.460388e+01
11	11	0.077351	8.207200e+02	8.206426e+02
12	12	0.071772	-9.847640e+03	9.847712e+03
13	13	0.066947	1.280203e+05	1.280203e+05
14	14	0.062731	-1.792283e+06	1.792284e+06

r = 6 #有效位数变成6
approx_vec = In_approx(n, r=r)
err_vec = np.abs(approx_vec - In_exact_value)
pd.DataFrame({"n": np.array(np.arange(n)), \
              "真实" : In_exact_value, "近似": approx_vec, "误差": err_vec})

I_0 使用近似值 0.6321210000

	n	真实	近似	误差
0	0	0.632120	0.632121	1.441171e-06
1	1	0.367878	0.367879	5.588274e-07
2	2	0.264240	0.264242	1.882341e-06
3	3	0.207276	0.207274	1.647031e-06
4	4	0.170892	0.170904	1.158811e-05
5	5	0.145532	0.145480	5.194058e-05
6	6	0.126801	0.127120	3.186434e-04
7	7	0.112383	0.110160	2.222504e-03
8	8	0.100931	0.118720	1.778903e-02
9	9	0.091611	-0.068480	1.600913e-01
10	10	0.083876	1.684800	1.600924e+00
11	11	0.077351	-17.532800	1.761015e+01
12	12	0.071772	211.393600	2.113218e+02
13	13	0.066947	-2747.116800	2.747184e+03
14	14	0.062731	38460.635194	3.846057e+04

例子2：大数吃小数#

例子来源于潘建瑜教授的讲义

(10**9 + 10**-9 - 10**9)/(10**-9)

0.0

(10**9 + 10**-9) == 10**9

True

例子3: 解线性方程组#

alpha = 0.99
A = np.matrix([[1, alpha],[alpha,1]])
b = np.array([1,0])
np.linalg.solve(A,b)

array([ 50.25125628, -49.74874372])

对于\(\alpha\)，一个小的扰动就可能造成解的误差很大。

该问题是源自线性方程组自身的问题，即矩阵的条件数过大（又叫做病态问题）。

alpha = 0.991
A = np.matrix([[1, alpha],[alpha,1]])
b = np.array([1,0])
np.linalg.solve(A,b)

array([ 55.80668564, -55.30442547])