摘要
数值分析是现代科学计算的基石，它研究如何使用数值近似方法解决数学分析问题。Python凭借其强大的科学生态系统（NumPy、SciPy、Matplotlib等），已成为数值计算的首选工具之一。本文将全面介绍数值分析的数学基础原理，并通过丰富的Python实战示例，详细讲解数值微积分、方程求解、线性代数、插值与拟合、优化方法等核心主题。文章包含大量表格对比和代码实例，旨在为读者提供一份系统、实用的数值分析学习指南。
1. 引言：为什么需要数值分析？
1.1 数值分析的意义
在科学和工程领域，我们经常遇到无法解析求解的数学问题：
-
无解析解的问题：大多数微分方程、高阶多项式方程和复杂积分没有解析解
-
大规模数据处理：现代科学实验产生的海量数据需要数值方法进行处理分析
-
实时计算需求：许多应用场景需要快速数值解而非精确的解析表达式
1.2 Python在数值分析中的优势
Python已成为数值分析的首选语言，主要原因包括：
| 优势 | 说明 | 相关库 |
|---|---|---|
| 丰富的科学生态 | 提供从基础计算到高级建模的完整工具链 | NumPy, SciPy, Pandas |
| 简洁的语法 | 代码更接近数学表达式，易于理解和维护 | 原生Python语法 |
| 高性能计算 | 底层使用C/Fortran实现，速度接近编译语言 | NumPy, Numba, Cython |
| 可视化能力 | 强大的数据可视化工具 | Matplotlib, Seaborn, Plotly |
| 社区支持 | 庞大的用户社区和丰富的学习资源 | 开源社区 |
1.3 本文结构概述
本文将按照以下结构组织内容：
-
数值分析基础与误差分析
-
非线性方程求解
-
线性代数计算方法
-
插值与多项式逼近
-
数值微积分
-
常微分方程数值解
-
数值优化方法
-
综合实战案例
每个部分都将包含数学原理讲解、算法实现和Python实战示例。
2. 数值分析基础与误差分析
2.1 数值误差的类型与来源
在数值计算中，误差是不可避免的。了解误差来源对于提高计算精度至关重要。
| 误差类型 | 描述 | 示例 |
|---|---|---|
| 模型误差 | 数学模型与实际问题之间的差异 | 忽略空气阻力的自由落体模型 |
| 观测误差 | 实验测量数据的不精确性 | 仪器测量精度限制 |
| 截断误差 | 用有限过程近似无限过程引起的误差 | 泰勒级数截断 |
| 舍入误差 | 计算机有限精度表示数字导致的误差 | π的浮点数表示 |
2.2 误差的度量指标
衡量数值计算准确度的常用指标：
| 指标 | 公式 | 描述 |
|---|---|---|
| 绝对误差 | Eabs=∥x^−x∥Eabs=∥x^−x∥ | 估计值与真实值之差的绝对值 |
| 相对误差 | Erel=∥x^−x∥∥x∥Erel=∥x∥∥x^−x∥ | 绝对误差与真实值的比值 |
| 百分比误差 | Epercent=Erel×100%Epercent=Erel×100% | 相对误差的百分比表示 |
2.3 Python中的浮点数精度
Python使用IEEE 754标准表示浮点数，了解其特性对数值计算至关重要：
python
import numpy as np
import sys
# 查看浮点数精度信息
print("Python浮点数精度信息:")
print(f"机器精度: {np.finfo(float).eps}")
print(f"最大浮点数: {np.finfo(float).max}")
print(f"最小正浮点数: {np.finfo(float).tiny}")
# 演示浮点数精度问题
a = 1.0
b = 1e-16
print(f"\n精度问题演示:")
print(f"1.0 + 1e-16 == 1.0? {1.0 + 1e-16 == 1.0}")
# catastrophic cancellation示例
def catastrophic_cancellation():
x = 1e10
y = 1e-10
z = -1e10
# 直接计算 (x + y + z) 应该是 y
result1 = x + y + z
result2 = y + x + z
print(f"直接计算 (x + y + z) = {result1}")
print(f"重新排列 (y + x + z) = {result2}")
print(f"期望结果 = {y}")
catastrophic_cancellation()
2.4 数值稳定性和条件数
数值稳定性是算法设计中的重要概念：
python
# 条件数示例：矩阵求逆的敏感性
def matrix_condition_demo():
# 良条件矩阵
A_well = np.array([[1, 0.5], [0.5, 1.33]])
# 病条件矩阵
A_ill = np.array([[1, 0.99], [0.99, 0.98]])
# 计算条件数
cond_well = np.linalg.cond(A_well)
cond_ill = np.linalg.cond(A_ill)
print(f"良条件矩阵的条件数: {cond_well:.4f}")
print(f"病条件矩阵的条件数: {cond_ill:.4f}")
# 演示微小扰动对解的影响
b = np.array([1, 1])
# 原始解
x_well = np.linalg.solve(A_well, b)
x_ill = np.linalg.solve(A_ill, b)
# 添加微小扰动
A_ill_perturbed = A_ill + np.random.normal(0, 1e-6, A_ill.shape)
x_ill_perturbed = np.linalg.solve(A_ill_perturbed, b)
error = np.linalg.norm(x_ill - x_ill_perturbed)
print(f"病条件矩阵解的误差: {error:.6f}")
matrix_condition_demo()
2.5 数值计算最佳实践
提高数值计算精度的实用技巧：
| 技巧 | 描述 | 示例 |
|---|---|---|
| 避免相近数相减 | 使用代数变换避免精度损失 | x+1−x=1x+1+xx+1 |
−x=x+1+x
| 1 | ||
| 避免大数吃小数 | 调整计算顺序，先处理小数值 | 求和时从小到大排序 |
| 使用更高精度 | 必要时使用高精度数据类型 | Python的decimal模块 |
| 稳定性算法选择 | 选择数值稳定的算法 | 使用QR分解而非正规方程 |
3. 非线性方程求解
非线性方程求解是数值分析的核心问题之一，形式为 f(x)=0f(x)=0。
3.1 二分法
二分法是最简单直观的求根方法，基于介值定理。
算法原理：
-
选择初始区间 [a,b][a,b]，确保 f(a)⋅f(b)<0f(a)⋅f(b)<0
-
计算中点 c=a+b2c=2a+b
-
根据 f(c)f(c) 的符号确定新的区间
-
重复直到满足精度要求
python
def bisection_method(f, a, b, tol=1e-6, max_iter=100):
"""
二分法求根实现
参数:
f: 目标函数
a, b: 初始区间
tol: 容差
max_iter: 最大迭代次数
返回:
根的近似值
"""
if f(a) * f(b) >= 0:
raise ValueError("函数在区间端点必须异号")
history = [] # 记录迭代历史
for i in range(max_iter):
c = (a + b) / 2
history.append((i, a, b, c, abs(b - a)))
if abs(b - a) < tol or abs(f(c)) < tol:
break
if f(a) * f(c) < 0:
b = c
else:
a = c
return c, history
# 示例：求解 f(x) = x^3 - 2x - 5 = 0
def example_function(x):
return x**3 - 2*x - 5
# 应用二分法
root, history = bisection_method(example_function, 2, 3, tol=1e-8)
print(f"二分法求得的根: {root:.8f}")
print(f"f({root:.8f}) = {example_function(root):.12f}")
# 输出迭代历史
print("\n迭代历史:")
print("Iteration | a | b | c | Error ")
print("-" * 55)
for iter, a, b, c, error in history:
print(f"{iter:9} | {a:.6f} | {b:.6f} | {c:.6f} | {error:.6e}")
3.2 牛顿-拉弗森方法
牛顿法利用函数的一阶导数信息，具有更快的收敛速度。
算法原理：
-
选择初始近似值 x0x0
-
迭代公式: xn+1=xn−f(xn)f′(xn)xn+1=xn−f′(xn)f(xn)
-
重复直到收敛
python
def newton_method(f, f_prime, x0, tol=1e-6, max_iter=100):
"""
牛顿法求根实现
参数:
f: 目标函数
f_prime: 目标函数的导数
x0: 初始猜测值
tol: 容差
max_iter: 最大迭代次数
返回:
根的近似值
"""
x = x0
history = []
for i in range(max_iter):
fx = f(x)
fpx = f_prime(x)
if abs(fpx) < tol: # 避免除零错误
break
x_new = x - fx / fpx
error = abs(x_new - x)
history.append((i, x, fx, fpx, error))
if error < tol or abs(fx) < tol:
break
x = x_new
return x, history
# 定义函数和导数
def f_example(x):
return x**3 - 2*x - 5
def f_prime_example(x):
return 3*x**2 - 2
# 应用牛顿法
root, history = newton_method(f_example, f_prime_example, 2.0, tol=1e-10)
print(f"牛顿法求得的根: {root:.10f}")
print(f"f({root:.10f}) = {f_example(root):.12f}")
# 输出迭代历史
print("\n迭代历史:")
print("Iteration | x_n | f(x_n) | f'(x_n) | Error ")
print("-" * 65)
for iter, x, fx, fpx, error in history:
print(f"{iter:9} | {x:.8f} | {fx:.8f} | {fpx:.8f} | {error:.6e}")
3.3 割线法
当导数难以计算时，可以使用割线法，它用差商近似导数。
算法原理：
-
选择两个初始近似值 x0,x1x0,x1
-
迭代公式: xn+1=xn−f(xn)xn−xn−1f(xn)−f(xn−1)xn+1=xn−f(xn)f(xn)−f(xn−1)xn−xn−1
-
重复直到收敛
python
def secant_method(f, x0, x1, tol=1e-6, max_iter=100):
"""
割线法求根实现
参数:
f: 目标函数
x0, x1: 初始猜测值
tol: 容差
max_iter: 最大迭代次数
返回:
根的近似值
"""
history = []
fx0 = f(x0)
fx1 = f(x1)
for i in range(max_iter):
if abs(fx1 - fx0) < tol: # 避免除零错误
break
# 割线法公式
x_new = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
error = abs(x_new - x1)
history.append((i, x1, fx1, error))
if error < tol or abs(fx1) < tol:
break
# 更新变量
x0, x1 = x1, x_new
fx0, fx1 = fx1, f(x_new)
return x1, history
# 应用割线法
root, history = secant_method(f_example, 2.0, 2.1, tol=1e-10)
print(f"割线法求得的根: {root:.10f}")
print(f"f({root:.10f}) = {f_example(root):.12f}")
# 输出迭代历史
print("\n迭代历史:")
print("Iteration | x_n | f(x_n) | Error ")
print("-" * 50)
for iter, x, fx, error in history:
print(f"{iter:9} | {x:.8f} | {fx:.8f} | {error:.6e}")
3.4 方法比较与选择
不同求根方法的特性比较：
| 方法 | 收敛速度 | 需要导数 | 初始条件 | 稳定性 | 适用场景 |
|---|---|---|---|---|---|
| 二分法 | 线性(1) | 不需要 | 区间[a,b], f(a)f(b)<0 | 高 | 简单问题，保证收敛 |
| 牛顿法 | 二次(2) | 需要 | 初始近似值x₀ | 中等 | 导数易求，初始值接近根 |
| 割线法 | 超线性(1.618) | 不需要 | 两个初始近似值 | 中等 | 导数难求或计算成本高 |
| 混合方法 | 可变 | 可选 | 依赖具体实现 | 高 | 实际应用，结合多种方法优点 |
4. 线性代数计算方法
线性代数问题是数值分析的核心内容，包括线性方程组求解、特征值计算等。
4.1 直接法：高斯消元与LU分解
4.1.1 高斯消元法
高斯消元是解决线性方程组的基础方法。
python
def gaussian_elimination(A, b):
"""
高斯消元法求解线性方程组 Ax = b
参数:
A: 系数矩阵 (n x n)
b: 右侧向量 (n)
返回:
解向量 x
"""
n = len(b)
# 构造增广矩阵
Ab = np.hstack([A.astype(float), b.astype(float).reshape(-1, 1)])
# 前向消元
for i in range(n):
# 部分主元选取
max_row = i + np.argmax(np.abs(Ab[i:, i]))
Ab[[i, max_row]] = Ab[[max_row, i]]
# 消元
for j in range(i + 1, n):
factor = Ab[j, i] / Ab[i, i]
Ab[j, i:] -= factor * Ab[i, i:]
# 回代
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = (Ab[i, n] - np.dot(Ab[i, i+1:n], x[i+1:])) / Ab[i, i]
return x
# 测试高斯消元法
A = np.array([[3, 2, -1],
[2, -2, 4],
[-1, 0.5, -1]], dtype=float)
b = np.array([1, -2, 0], dtype=float)
x = gaussian_elimination(A, b)
print("高斯消元法解:")
print(f"x = {x}")
print(f"验证 Ax = b: {np.allclose(np.dot(A, x), b)}")
4.1.2 LU分解法
LU分解将矩阵分解为下三角和上三角矩阵的乘积，适合多次求解相同系数矩阵的问题。
python
def lu_decomposition(A):
"""
LU分解 (Doolittle算法)
参数:
A: 方阵
返回:
L: 下三角矩阵
U: 上三角矩阵
"""
n = A.shape[0]
L = np.eye(n) # 单位下三角矩阵
U = np.zeros_like(A)
for i in range(n):
# 计算U的第i行
for k in range(i, n):
U[i, k] = A[i, k] - np.dot(L[i, :i], U[:i, k])
# 计算L的第i列
for k in range(i+1, n):
L[k, i] = (A[k, i] - np.dot(L[k, :i], U[:i, i])) / U[i, i]
return L, U
def solve_lu(L, U, b):
"""
使用LU分解求解线性方程组
参数:
L: 下三角矩阵
U: 上三角矩阵
b: 右侧向量
返回:
解向量 x
"""
n = len(b)
# 前向替换: Ly = b
y = np.zeros(n)
for i in range(n):
y[i] = b[i] - np.dot(L[i, :i], y[:i])
# 回代: Ux = y
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]
return x
# 测试LU分解
L, U = lu_decomposition(A)
print("L矩阵:")
print(L)
print("U矩阵:")
print(U)
print("验证 LU = A:", np.allclose(np.dot(L, U), A))
# 使用LU分解求解
x_lu = solve_lu(L, U, b)
print(f"LU分解解: x = {x_lu}")
4.2 迭代法：雅可比与高斯-赛德尔
对于大型稀疏矩阵，迭代法比直接法更有效。
4.2.1 雅可比迭代法
python
def jacobi_iteration(A, b, x0=None, tol=1e-6, max_iter=1000):
"""
雅可比迭代法求解线性方程组
参数:
A: 系数矩阵
b: 右侧向量
x0: 初始猜测
tol: 容差
max_iter: 最大迭代次数
返回:
解向量 x
"""
n = len(b)
if x0 is None:
x0 = np.zeros(n)
x = x0.copy()
x_new = np.zeros(n)
history = []
for k in range(max_iter):
for i in range(n):
s = 0
for j in range(n):
if j != i:
s += A[i, j] * x[j]
x_new[i] = (b[i] - s) / A[i, i]
error = np.linalg.norm(x_new - x)
history.append((k, error))
if error < tol:
break
x = x_new.copy()
return x_new, history
# 测试雅可比迭代（对角占优矩阵确保收敛）
A_jacobi = np.array([[10, -1, 2, 0],
[-1, 11, -1, 3],
[2, -1, 10, -1],
[0, 3, -1, 8]], dtype=float)
b_jacobi = np.array([6, 25, -11, 15], dtype=float)
x_jacobi, history = jacobi_iteration(A_jacobi, b_jacobi, tol=1e-10)
print("雅可比迭代解:")
print(f"x = {x_jacobi}")
print(f"验证 Ax = b: {np.allclose(np.dot(A_jacobi, x_jacobi), b_jacobi)}")
4.2.2 高斯-赛德尔迭代法
高斯-赛德尔方法使用最新计算的值，通常比雅可比方法收敛更快。
python
def gauss_seidel_iteration(A, b, x0=None, tol=1e-6, max_iter=1000):
"""
高斯-赛德尔迭代法求解线性方程组
参数:
A: 系数矩阵
b: 右侧向量
x0: 初始猜测
tol: 容差
max_iter: 最大迭代次数
返回:
解向量 x
"""
n = len(b)
if x0 is None:
x0 = np.zeros(n)
x = x0.copy()
history = []
for k in range(max_iter):
x_old = x.copy()
for i in range(n):
s1 = np.dot(A[i, :i], x[:i])
s2 = np.dot(A[i, i+1:], x_old[i+1:])
x[i] = (b[i] - s1 - s2) / A[i, i]
error = np.linalg.norm(x - x_old)
history.append((k, error))
if error < tol:
break
return x, history
# 测试高斯-赛德尔迭代
x_gs, history_gs = gauss_seidel_iteration(A_jacobi, b_jacobi, tol=1e-10)
print("高斯-赛德尔迭代解:")
print(f"x = {x_gs}")
print(f"验证 Ax = b: {np.allclose(np.dot(A_jacobi, x_gs), b_jacobi)}")
4.3 矩阵分解高级方法
4.3.1 QR分解
QR分解将矩阵分解为正交矩阵和上三角矩阵的乘积。
python
def householder_qr(A):
"""
使用Householder变换进行QR分解
参数:
A: 矩阵 (m x n)
返回:
Q: 正交矩阵 (m x m)
R: 上三角矩阵 (m x n)
"""
m, n = A.shape
Q = np.eye(m)
R = A.copy().astype(float)
for j in range(min(m, n)):
# 计算Householder向量
x = R[j:, j]
norm_x = np.linalg.norm(x)
sign = 1 if x[0] >= 0 else -1
v = x.copy()
v[0] += sign * norm_x
v = v / np.linalg.norm(v)
# 应用Householder变换
R[j:, j:] -= 2 * np.outer(v, np.dot(v, R[j:, j:]))
# 累积正交变换
Q_j = np.eye(m)
Q_j[j:, j:] -= 2 * np.outer(v, v)
Q = Q @ Q_j.T
return Q, R
# 测试QR分解
A_qr = np.array([[12, -51, 4],
[6, 167, -68],
[-4, 24, -41]], dtype=float)
Q, R = householder_qr(A_qr)
print("Q矩阵:")
print(Q)
print("R矩阵:")
print(R)
print("验证 QR = A:", np.allclose(Q @ R, A_qr))
print("验证 Q是正交矩阵:", np.allclose(Q @ Q.T, np.eye(Q.shape[0])))
4.3.2 特征值计算：幂迭代法
幂迭代法是计算矩阵主特征值的简单方法。
python
def power_iteration(A, num_iterations=100):
"""
幂迭代法计算矩阵的主特征值和特征向量
参数:
A: 方阵
num_iterations: 迭代次数
返回:
主特征值，对应的特征向量
"""
n = A.shape[0]
b_k = np.random.rand(n)
for _ in range(num_iterations):
# 计算矩阵乘向量
b_k1 = np.dot(A, b_k)
# 计算范数
b_k1_norm = np.linalg.norm(b_k1)
# 重新归一化
b_k = b_k1 / b_k1_norm
# 计算特征值 (瑞利商)
eigenvalue = np.dot(b_k.T, np.dot(A, b_k))
return eigenvalue, b_k
# 测试幂迭代法
A_eig = np.array([[2, 1],
[1, 2]], dtype=float)
eigenvalue, eigenvector = power_iteration(A_eig, 100)
print(f"主特征值: {eigenvalue:.6f}")
print(f"特征向量: {eigenvector}")
# 与numpy计算结果比较
eigenvalues_np, eigenvectors_np = np.linalg.eig(A_eig)
print(f"NumPy计算的特征值: {eigenvalues_np}")
print(f"NumPy计算的特征向量:\n{eigenvectors_np}")
4.4 稀疏矩阵处理
对于大型稀疏矩阵，使用专用数据结构可以大幅提高计算效率。
python
from scipy import sparse
def sparse_matrix_demo():
"""稀疏矩阵处理示例"""
# 创建稠密矩阵
n = 1000
density = 0.01 # 1% 非零元素
A_dense = np.random.rand(n, n)
A_dense[A_dense > density] = 0 # 使矩阵稀疏化
# 转换为稀疏矩阵格式
A_csr = sparse.csr_matrix(A_dense)
print(f"稠密矩阵大小: {A_dense.nbytes / 1024 / 1024:.2f} MB")
print(f"稀疏矩阵(CSR)大小: {A_csr.data.nbytes / 1024:.2f} KB")
print(f"压缩比: {A_dense.nbytes / A_csr.data.nbytes:.1f}x")
# 稀疏矩阵运算
x = np.random.rand(n)
# 比较计算时间
import time
# 稠密矩阵向量乘法
start = time.time()
b_dense = A_dense @ x
time_dense = time.time() - start
# 稀疏矩阵向量乘法
start = time.time()
b_sparse = A_csr @ x
time_sparse = time.time() - start
print(f"\n稠密矩阵乘法时间: {time_dense:.6f}秒")
print(f"稀疏矩阵乘法时间: {time_sparse:.6f}秒")
print(f"加速比: {time_dense / time_sparse:.1f}x")
# 验证结果一致性
print(f"结果一致性: {np.allclose(b_dense, b_sparse)}")
sparse_matrix_demo()
5. 插值与多项式逼近
插值是通过已知数据点构造函数的过程，多项式逼近是用多项式函数近似复杂函数。
5.1 拉格朗日插值
拉格朗日插值通过构造基函数实现插值。
python
def lagrange_interpolation(x_points, y_points, x):
"""
拉格朗日插值
参数:
x_points: 已知点的x坐标
y_points: 已知点的y坐标
x: 待插值点的x坐标
返回:
插值结果
"""
n = len(x_points)
result = 0.0
for i in range(n):
# 计算拉格朗日基函数
term = y_points[i]
for j in range(n):
if i != j:
term *= (x - x_points[j]) / (x_points[i] - x_points[j])
result += term
return result
# 测试拉格朗日插值
def runge_function(x):
return 1 / (1 + 25 * x**2)
# 生成插值节点
n_points = 11
x_nodes = np.linspace(-1, 1, n_points)
y_nodes = runge_function(x_nodes)
# 生成密集评估点
x_eval = np.linspace(-1, 1, 1000)
y_true = runge_function(x_eval)
y_interp = np.array([lagrange_interpolation(x_nodes, y_nodes, x) for x in x_eval])
# 计算误差
error = np.abs(y_true - y_interp)
max_error = np.max(error)
mean_error = np.mean(error)
print(f"拉格朗日插值误差:")
print(f"最大误差: {max_error:.6f}")
print(f"平均误差: {mean_error:.6f}")
5.2 牛顿插值
牛顿插值使用差商表，便于添加新节点。
python
def newton_interpolation(x_points, y_points, x):
"""
牛顿插值
参数:
x_points: 已知点的x坐标
y_points: 已知点的y坐标
x: 待插值点的x坐标
返回:
插值结果
"""
n = len(x_points)
# 计算差商表
F = np.zeros((n, n))
F[:, 0] = y_points
for j in range(1, n):
for i in range(n - j):
F[i, j] = (F[i+1, j-1] - F[i, j-1]) / (x_points[i+j] - x_points[i])
# 计算插值结果
result = F[0, 0]
product = 1.0
for i in range(1, n):
product *= (x - x_points[i-1])
result += F[0, i] * product
return result
# 测试牛顿插值
y_newton = np.array([newton_interpolation(x_nodes, y_nodes, x) for x in x_eval])
# 验证与拉格朗日插值结果一致
print(f"牛顿与拉格朗日插值一致性: {np.allclose(y_interp, y_newton, atol=1e-10)}")
5.3 样条插值
样条插值使用分段低次多项式，避免高次插值的Runge现象。
5.3.1 三次样条插值
python
def cubic_spline(x_points, y_points, x_eval):
"""
三次样条插值 (自然边界条件)
参数:
x_points: 节点x坐标
y_points: 节点y坐标
x_eval: 评估点
返回:
插值结果
"""
n = len(x_points)
# 计算步长
h = np.diff(x_points)
# 构建三对角方程组
A = np.zeros((n, n))
b = np.zeros(n)
# 内部点方程
for i in range(1, n-1):
A[i, i-1] = h[i-1]
A[i, i] = 2 * (h[i-1] + h[i])
A[i, i+1] = h[i]
b[i] = 3 * ((y_points[i+1] - y_points[i]) / h[i] -
(y_points[i] - y_points[i-1]) / h[i-1])
# 自然边界条件
A[0, 0] = 1
A[n-1, n-1] = 1
# 求解二阶导数
c = np.linalg.solve(A, b)
# 计算其他系数
a = y_points[:-1]
b_coeff = np.diff(y_points) / h - h * (2 * c[:-1] + c[1:]) / 3
d = np.diff(c) / (3 * h)
# 评估插值
y_eval = np.zeros_like(x_eval)
for i in range(n-1):
mask = (x_points[i] <= x_eval) & (x_eval <= x_points[i+1])
x_local = x_eval[mask] - x_points[i]
y_eval[mask] = (a[i] + b_coeff[i] * x_local +
c[i] * x_local**2 + d[i] * x_local**3)
return y_eval
# 测试三次样条插值
y_spline = cubic_spline(x_nodes, y_nodes, x_eval)
# 计算误差
error_spline = np.abs(y_true - y_spline)
max_error_spline = np.max(error_spline)
mean_error_spline = np.mean(error_spline)
print(f"三次样条插值误差:")
print(f"最大误差: {max_error_spline:.6f}")
print(f"平均误差: {mean_error_spline:.6f}")
# 与高次多项式插值比较
plt.figure(figsize=(10, 6))
plt.plot(x_eval, y_true, 'k-', label='真实函数')
plt.plot(x_eval, y_interp, 'r--', label='拉格朗日插值')
plt.plot(x_eval, y_spline, 'b-.', label='三次样条')
plt.plot(x_nodes, y_nodes, 'go', label='插值节点')
plt.legend()
plt.title('Runge现象: 高次多项式插值 vs 三次样条')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
5.4 最小二乘拟合
当数据有噪声时，最小二乘拟合比插值更合适。
python
def least_squares_fit(x, y, degree):
"""
多项式最小二乘拟合
参数:
x: 数据点x坐标
y: 数据点y坐标
degree: 多项式次数
返回:
拟合多项式系数
"""
# 构建范德蒙矩阵
A = np.vander(x, degree + 1)
# 求解正规方程
coeffs = np.linalg.solve(A.T @ A, A.T @ y)
return coeffs
def polynomial_eval(coeffs, x):
"""
计算多项式值
参数:
coeffs: 多项式系数 [最高次项系数, ..., 常数项]
x: 评估点
返回:
多项式值
"""
n = len(coeffs)
result = np.zeros_like(x)
for i in range(n):
result += coeffs[i] * x**(n-1-i)
return result
# 生成带噪声的数据
np.random.seed(42)
x_data = np.linspace(0, 2*np.pi, 20)
y_true = np.sin(x_data)
y_noisy = y_true + 0.1 * np.random.randn(len(x_data))
# 不同次数的多项式拟合
degrees = [3, 5, 7, 10]
x_fine = np.linspace(0, 2*np.pi, 1000)
plt.figure(figsize=(12, 8))
for i, degree in enumerate(degrees):
coeffs = least_squares_fit(x_data, y_noisy, degree)
y_fit = polynomial_eval(coeffs, x_fine)
plt.subplot(2, 2, i+1)
plt.plot(x_fine, np.sin(x_fine), 'k-', label='真实函数')
plt.plot(x_data, y_noisy, 'ro', label='噪声数据')
plt.plot(x_fine, y_fit, 'b-', label=f'{degree}次拟合')
plt.legend()
plt.title(f'{degree}次多项式拟合')
plt.grid(True)
plt.tight_layout()
plt.show()
# 计算不同次数拟合的误差
errors = []
for degree in range(1, 15):
coeffs = least_squares_fit(x_data, y_noisy, degree)
y_fit = polynomial_eval(coeffs, x_data)
error = np.mean((y_fit - y_true)**2)
errors.append(error)
plt.figure(figsize=(10, 5))
plt.plot(range(1, 15), errors, 'bo-')
plt.xlabel('多项式次数')
plt.ylabel('均方误差')
plt.title('过拟合现象: 误差随多项式次数的变化')
plt.grid(True)
plt.show()
6. 数值微积分
数值微积分解决无法解析求积分或导数的问题。
6.1 数值积分
6.1.1 牛顿-柯特斯公式
python
def newton_cotes_integration(f, a, b, n=100, method='simpson'):
"""
牛顿-柯特斯数值积分
参数:
f: 被积函数
a, b: 积分区间
n: 子区间数量(偶数)
method: 方法('trapezoid', 'simpson')
返回:
积分近似值
"""
if n % 2 != 0 and method == 'simpson':
n += 1 # 辛普森法则需要偶数区间
x = np.linspace(a, b, n + 1)
h = (b - a) / n
y = f(x)
if method == 'trapezoid':
# 梯形法则
integral = h * (0.5 * y[0] + 0.5 * y[-1] + np.sum(y[1:-1]))
elif method == 'simpson':
# 辛普森法则
integral = h / 3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]))
else:
raise ValueError("不支持的积分方法")
return integral
# 测试数值积分
def test_integral(x):
return np.exp(-x**2)
a, b = 0, 2
exact_integral = 0.8820813907624217 # 近似精确值
# 不同方法的积分结果
methods = ['trapezoid', 'simpson']
n_values = [10, 20, 50, 100, 200]
results = {}
for method in methods:
errors = []
for n in n_values:
integral = newton_cotes_integration(test_integral, a, b, n, method)
error = abs(integral - exact_integral)
errors.append(error)
results[method] = errors
# 显示结果
print("数值积分误差比较:")
print("n\t梯形法则误差\t辛普森法则误差")
print("-" * 45)
for i, n in enumerate(n_values):
print(f"{n}\t{results['trapezoid'][i]:.8f}\t{results['simpson'][i]:.8f}")
# 收敛速度分析
plt.figure(figsize=(10, 6))
for method in methods:
plt.loglog(n_values, results[method], 'o-', label=method)
plt.xlabel('子区间数量 n')
plt.ylabel('误差')
plt.title('数值积分方法的收敛速度')
plt.legend()
plt.grid(True)
plt.show()
6.1.2 高斯求积法
高斯求积使用最优节点和权重，达到最高代数精度。
python
def gauss_legendre_integration(f, a, b, n=5):
"""
高斯-勒让德求积
参数:
f: 被积函数
a, b: 积分区间
n: 积分点数量
返回:
积分近似值
"""
# 勒让德多项式的根和权重
if n == 2:
nodes = [-0.5773502692, 0.5773502692]
weights = [1.0, 1.0]
elif n == 3:
nodes = [-0.7745966692, 0.0, 0.7745966692]
weights = [0.5555555556, 0.8888888889, 0.5555555556]
elif n == 5:
nodes = [-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459]
weights = [0.2369268851, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268851]
else:
raise ValueError("不支持的点数")
# 坐标变换到[a, b]
transform = lambda x: (b - a) / 2 * x + (a + b) / 2
nodes_transformed = transform(np.array(nodes))
weights_transformed = np.array(weights) * (b - a) / 2
# 计算积分
integral = np.sum(weights_transformed * f(nodes_transformed))
return integral
# 测试高斯求积
gauss_results = {}
for n in [2, 3, 5]:
integral = gauss_legendre_integration(test_integral, a, b, n)
error = abs(integral - exact_integral)
gauss_results[n] = error
print("高斯求积误差:")
for n, error in gauss_results.items():
print(f"{n}点高斯求积误差: {error:.10f}")
# 与牛顿-柯特斯方法比较
simpson_error = results['simpson'][-1] # n=200的辛普森法则误差
print(f"\n比较: 5点高斯求积误差 {gauss_results[5]:.10f} vs 200点辛普森法则误差 {simpson_error:.10f}")
6.2 数值微分
6.2.1 有限差分法
python
def finite_difference(f, x, h=1e-5, method='central'):
"""
有限差分法数值求导
参数:
f: 函数
x: 求导点
h: 步长
method: 方法('forward', 'backward', 'central')
返回:
导数近似值
"""
if method == 'forward':
return (f(x + h) - f(x)) / h
elif method == 'backward':
return (f(x) - f(x - h)) / h
elif method == 'central':
return (f(x + h) - f(x - h)) / (2 * h)
else:
raise ValueError("不支持的差分方法")
# 测试数值微分
def test_function(x):
return np.sin(x) * np.exp(-x)
def exact_derivative(x):
return np.exp(-x) * (np.cos(x) - np.sin(x))
x0 = 1.0
exact_deriv = exact_derivative(x0)
# 不同步长的误差
h_values = [10**(-i) for i in range(1, 16)]
methods = ['forward', 'backward', 'central']
errors = {method: [] for method in methods}
for h in h_values:
for method in methods:
deriv_approx = finite_difference(test_function, x0, h, method)
error = abs(deriv_approx - exact_deriv)
errors[method].append(error)
# 绘制误差随步长变化
plt.figure(figsize=(10, 6))
for method in methods:
plt.loglog(h_values, errors[method], 'o-', label=method)
plt.xlabel('步长 h')
plt.ylabel('误差')
plt.title('数值微分误差随步长的变化')
plt.legend()
plt.grid(True)
plt.show()
# 最佳步长分析
optimal_h = {}
for method in methods:
min_error = min(errors[method])
optimal_index = errors[method].index(min_error)
optimal_h[method] = h_values[optimal_index]
print(f"{method}差分: 最佳步长 = {optimal_h[method]:.2e}, 最小误差 = {min_error:.2e}")
6.2.2 理查森外推法
理查森外推通过不同步长的组合提高精度。
python
def richardson_extrapolation(f, x, h, n=2):
"""
理查森外推法提高数值微分精度
参数:
f: 函数
x: 求导点
h: 初始步长
n: 外推次数
返回:
高精度导数近似值
"""
# 创建外推表
D = np.zeros((n+1, n+1))
# 第一列: 不同步长的中心差分
for i in range(n+1):
h_i = h / (2**i)
D[i, 0] = finite_difference(f, x, h_i, 'central')
# 外推计算
for j in range(1, n+1):
for i in range(j, n+1):
D[i, j] = D[i, j-1] + (D[i, j-1] - D[i-1, j-1]) / (4**j - 1)
return D[n, n]
# 测试理查森外推
h_initial = 0.1
richardson_deriv = richardson_extrapolation(test_function, x0, h_initial, n=4)
richardson_error = abs(richardson_deriv - exact_deriv)
print(f"理查森外推结果: {richardson_deriv:.10f}")
print(f"精确导数: {exact_deriv:.10f}")
print(f"理查森外推误差: {richardson_error:.2e}")
# 与简单中心差分比较
central_deriv = finite_difference(test_function, x0, h_initial, 'central')
central_error = abs(central_deriv - exact_deriv)
print(f"中心差分误差: {central_error:.2e}")
print(f"精度提升: {central_error/richardson_error:.1f}倍")
7. 常微分方程数值解
常微分方程(ODE)数值解是科学计算中的重要问题。
7.1 初值问题
7.1.1 欧拉方法
python
def euler_method(f, y0, t_range, h):
"""
欧拉方法求解ODE: dy/dt = f(t, y)
参数:
f: 右侧函数
y0: 初始条件
t_range: 时间范围 [t0, tf]
h: 步长
返回:
时间点和解
"""
t0, tf = t_range
t = np.arange(t0, tf + h, h)
n = len(t)
y = np.zeros(n)
y[0] = y0
for i in range(n-1):
y[i+1] = y[i] + h * f(t[i], y[i])
return t, y
# 测试欧拉方法
def ode_example(t, y):
return -2 * t * y # dy/dt = -2ty
def exact_solution(t):
return np.exp(-t**2) # 解析解: y = e^{-t^2}
# 求解ODE
t_range = [0, 2]
y0 = 1
h = 0.1
t_euler, y_euler = euler_method(ode_example, y0, t_range, h)
y_exact = exact_solution(t_euler)
# 计算误差
error = np.abs(y_euler - y_exact)
max_error = np.max(error)
print(f"欧拉方法最大误差: {max_error:.6f}")
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t_euler, y_exact, 'k-', label='解析解')
plt.plot(t_euler, y_euler, 'ro--', label='欧拉方法')
plt.xlabel('t')
plt.ylabel('y')
plt.title('欧拉方法求解ODE')
plt.legend()
plt.grid(True)
plt.show()
7.1.2 龙格-库塔方法
龙格-库塔方法比欧拉方法更精确。
python
def runge_kutta_4(f, y0, t_range, h):
"""
四阶龙格-库塔方法
参数:
f: 右侧函数
y0: 初始条件
t_range: 时间范围 [t0, tf]
h: 步长
返回:
时间点和解
"""
t0, tf = t_range
t = np.arange(t0, tf + h, h)
n = len(t)
y = np.zeros(n)
y[0] = y0
for i in range(n-1):
k1 = h * f(t[i], y[i])
k2 = h * f(t[i] + h/2, y[i] + k1/2)
k3 = h * f(t[i] + h/2, y[i] + k2/2)
k4 = h * f(t[i] + h, y[i] + k3)
y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
return t, y
# 测试龙格-库塔方法
t_rk4, y_rk4 = runge_kutta_4(ode_example, y0, t_range, h)
y_exact_rk4 = exact_solution(t_rk4)
# 计算误差
error_rk4 = np.abs(y_rk4 - y_exact_rk4)
max_error_rk4 = np.max(error_rk4)
print(f"四阶龙格-库塔方法最大误差: {max_error_rk4:.6f}")
# 与欧拉方法比较
print(f"精度提升: {max_error/max_error_rk4:.1f}倍")
# 绘制比较
plt.figure(figsize=(10, 6))
plt.plot(t_euler, y_exact, 'k-', label='解析解')
plt.plot(t_euler, y_euler, 'ro--', label='欧拉方法', markersize=4)
plt.plot(t_rk4, y_rk4, 'bx--', label='四阶RK方法', markersize=4)
plt.xlabel('t')
plt.ylabel('y')
plt.title('ODE求解方法比较')
plt.legend()
plt.grid(True)
plt.show()
7.2 刚性问题
刚性问题需要特殊方法处理，如隐式方法。
python
def implicit_euler(f, jacobian, y0, t_range, h):
"""
隐式欧拉方法求解刚性ODE
参数:
f: 右侧函数
jacobian: 雅可比矩阵函数
y0: 初始条件
t_range: 时间范围
h: 步长
返回:
时间点和解
"""
t0, tf = t_range
t = np.arange(t0, tf + h, h)
n = len(t)
y = np.zeros(n)
y[0] = y0
# 牛顿迭代求解隐式方程
for i in range(n-1):
y_guess = y[i] # 初始猜测
# 牛顿迭代
for _ in range(10): # 最大迭代次数
F = y_guess - y[i] - h * f(t[i+1], y_guess)
if np.linalg.norm(F) < 1e-8:
break
J = np.eye(len(y0)) - h * jacobian(t[i+1], y_guess)
delta = np.linalg.solve(J, -F)
y_guess += delta
y[i+1] = y_guess
return t, y
# 刚性ODE示例
def stiff_ode(t, y):
return -1000 * y + 3000 - 2000 * np.exp(-t)
def stiff_jacobian(t, y):
return np.array([[-1000]])
def stiff_exact(t):
return 3 - 0.998 * np.exp(-1000*t) - 2.002 * np.exp(-t)
# 求解刚性ODE
t_range_stiff = [0, 0.1]
y0_stiff = np.array([0])
h_stiff = 0.0001
# 显式欧拉方法会失效
try:
t_euler_stiff, y_euler_stiff = euler_method(stiff_ode, y0_stiff, t_range_stiff, h_stiff)
except Exception as e:
print(f"显式欧拉方法失败: {e}")
# 隐式欧拉方法
t_implicit, y_implicit = implicit_euler(stiff_ode, stiff_jacobian, y0_stiff, t_range_stiff, h_stiff)
y_exact_stiff = stiff_exact(t_implicit)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t_implicit, y_exact_stiff, 'k-', label='解析解')
plt.plot(t_implicit, y_implicit, 'ro--', label='隐式欧拉', markersize=4)
plt.xlabel('t')
plt.ylabel('y')
plt.title('刚性ODE的数值解')
plt.legend()
plt.grid(True)
plt.show()
8. 数值优化方法
数值优化寻找函数的最小值或最大值。
8.1 无约束优化
8.1.1 梯度下降法
python
def gradient_descent(f, grad_f, x0, learning_rate=0.1, max_iter=1000, tol=1e-6):
"""
梯度下降法
参数:
f: 目标函数
grad_f: 梯度函数
x0: 初始点
learning_rate: 学习率
max_iter: 最大迭代次数
tol: 容差
返回:
最优解和优化历史
"""
x = x0.copy()
history = {'x': [x], 'f': [f(x)]}
for i in range(max_iter):
gradient = grad_f(x)
x_new = x - learning_rate * gradient
history['x'].append(x_new)
history['f'].append(f(x_new))
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
return x, history
# 测试梯度下降
def quadratic_function(x):
return x[0]**2 + 10*x[1]**2 + 5*x[0]*x[1] + 2*x[0] - 3*x[1]
def quadratic_gradient(x):
return np.array([2*x[0] + 5*x[1] + 2, 20*x[1] + 5*x[0] - 3])
# 优化
x0 = np.array([0, 0])
x_opt, history = gradient_descent(quadratic_function, quadratic_gradient, x0, learning_rate=0.05)
print(f"最优解: {x_opt}")
print(f"最优值: {quadratic_function(x_opt)}")
# 绘制优化过程
x_history = np.array(history['x'])
f_history = history['f']
# 创建等高线图
x1 = np.linspace(-2, 2, 100)
x2 = np.linspace(-1, 1, 100)
X1, X2 = np.meshgrid(x1, x2)
Z = quadratic_function([X1, X2])
plt.figure(figsize=(10, 6))
plt.contour(X1, X2, Z, 50, cmap='viridis')
plt.plot(x_history[:, 0], x_history[:, 1], 'ro-', markersize=4)
plt.colorbar()
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('梯度下降优化路径')
plt.grid(True)
plt.show()
# 绘制函数值下降曲线
plt.figure(figsize=(10, 5))
plt.semilogy(f_history)
plt.xlabel('迭代次数')
plt.ylabel('函数值')
plt.title('梯度下降收敛曲线')
plt.grid(True)
plt.show()
8.1.2 牛顿法优化
python
def newton_optimization(f, grad_f, hess_f, x0, max_iter=100, tol=1e-6):
"""
牛顿法优化
参数:
f: 目标函数
grad_f: 梯度函数
hess_f: 海森矩阵函数
x0: 初始点
max_iter: 最大迭代次数
tol: 容差
返回:
最优解和优化历史
"""
x = x0.copy()
history = {'x': [x], 'f': [f(x)]}
for i in range(max_iter):
gradient = grad_f(x)
hessian = hess_f(x)
# 求解牛顿方程
try:
step = np.linalg.solve(hessian, -gradient)
except np.linalg.LinAlgError:
# 海森矩阵奇异，使用梯度下降
step = -gradient / np.linalg.norm(gradient)
x_new = x + step
history['x'].append(x_new)
history['f'].append(f(x_new))
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
return x, history
# 测试牛顿法
def quadratic_hessian(x):
return np.array([[2, 5], [5, 20]])
x_newton, history_newton = newton_optimization(quadratic_function, quadratic_gradient,
quadratic_hessian, x0)
print(f"牛顿法最优解: {x_newton}")
print(f"牛顿法最优值: {quadratic_function(x_newton)}")
# 比较梯度下降和牛顿法
plt.figure(figsize=(10, 5))
plt.semilogy(history['f'], label='梯度下降')
plt.semilogy(history_newton['f'], label='牛顿法')
plt.xlabel('迭代次数')
plt.ylabel('函数值')
plt.title('优化方法收敛速度比较')
plt.legend()
plt.grid(True)
plt.show()
8.2 约束优化
8.2.1 拉格朗日乘子法
python
def lagrangian_method(f, g, grad_f, grad_g, x0, max_iter=100, tol=1e-6):
"""
拉格朗日乘子法求解等式约束优化
参数:
f: 目标函数
g: 约束函数 (返回向量)
grad_f: 目标函数梯度
grad_g: 约束函数雅可比矩阵
x0: 初始点
max_iter: 最大迭代次数
tol: 容差
返回:
最优解和拉格朗日乘子
"""
x = x0.copy()
n = len(x0) # 变量维度
m = len(g(x0)) # 约束个数
# 初始拉格朗日乘子
lambda_val = np.zeros(m)
for i in range(max_iter):
# 构建KKT系统
gradient = grad_f(x) + grad_g(x).T @ lambda_val
constraint = g(x)
# 检查收敛
if np.linalg.norm(gradient) < tol and np.linalg.norm(constraint) < tol:
break
# 构建KKT矩阵
KKT_matrix = np.block([
[grad_g(x).T @ grad_g(x), grad_g(x).T],
[grad_g(x), np.zeros((m, m))]
])
# 构建右侧向量
rhs = np.concatenate([-grad_g(x).T @ grad_f(x), -g(x)])
# 求解KKT系统
solution = np.linalg.solve(KKT_matrix, rhs)
# 更新变量和乘子
x += solution[:n]
lambda_val += solution[n:]
return x, lambda_val
# 测试约束优化
def constraint_optimization_example():
# 目标函数: f(x) = x1^2 + x2^2
def f(x):
return x[0]**2 + x[1]**2
# 约束: g(x) = x1 + x2 - 1 = 0
def g(x):
return np.array([x[0] + x[1] - 1])
def grad_f(x):
return np.array([2*x[0], 2*x[1]])
def grad_g(x):
return np.array([[1, 1]]) # 雅可比矩阵
# 求解
x0 = np.array([0, 0])
x_opt, lambda_opt = lagrangian_method(f, g, grad_f, grad_g, x0)
print(f"约束优化最优解: {x_opt}")
print(f"目标函数值: {f(x_opt)}")
print(f"约束值: {g(x_opt)}")
print(f"拉格朗日乘子: {lambda_opt}")
constraint_optimization_example()
9. 综合实战案例
9.1 物理仿真：弹簧质点系统
python
def spring_mass_system():
"""
弹簧质点系统仿真
"""
# 系统参数
m = 1.0 # 质量
k = 10.0 # 弹簧系数
c = 0.5 # 阻尼系数
# ODE: m * x'' + c * x' + k * x = 0
def spring_ode(t, y):
x, v = y # 位移和速度
dxdt = v
dvdt = (-c * v - k * x) / m
return np.array([dxdt, dvdt])
# 初始条件
y0 = np.array([1.0, 0.0]) # 初始位移1m，初始速度0
t_range = [0, 10] # 仿真时间
h = 0.01 # 步长
# 使用四阶龙格-库塔方法求解
t, y = runge_kutta_4(spring_ode, y0, t_range, h)
x, v = y[:, 0], y[:, 1]
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(t, x, 'b-')
plt.xlabel('时间 (s)')
plt.ylabel('位移 (m)')
plt.title('弹簧质点系统位移响应')
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(t, v, 'r-')
plt.xlabel('时间 (s)')
plt.ylabel('速度 (m/s)')
plt.title('弹簧质点系统速度响应')
plt.grid(True)
plt.tight_layout()
plt.show()
# 相图
plt.figure(figsize=(8, 6))
plt.plot(x, v)
plt.xlabel('位移 (m)')
plt.ylabel('速度 (m/s)')
plt.title('弹簧质点系统相图')
plt.grid(True)
plt.show()
spring_mass_system()
9.2 数据分析：曲线拟合与预测
python
def data_fitting_prediction():
"""
数据曲线拟合与预测示例
"""
# 生成示例数据 (带有噪声的正弦曲线)
np.random.seed(42)
x_data = np.linspace(0, 2*np.pi, 50)
y_true = 2 * np.sin(1.5 * x_data) + 0.5 * np.cos(3 * x_data)
y_noisy = y_true + 0.3 * np.random.randn(len(x_data))
# 使用样条插值拟合数据
from scipy import interpolate
# 三次样条插值
spline = interpolate.CubicSpline(x_data, y_noisy)
# 生成预测点
x_pred = np.linspace(0, 2*np.pi, 500)
y_spline = spline(x_pred)
# 使用傅里叶级数拟合
def fourier_series(x, *coeffs):
result = coeffs[0] * np.ones_like(x)
n_terms = (len(coeffs) - 1) // 2
for i in range(n_terms):
result += coeffs[2*i+1] * np.sin((i+1)*x)
result += coeffs[2*i+2] * np.cos((i+1)*x)
return result
# 非线性最小二乘拟合
from scipy.optimize import curve_fit
# 初始猜测 (常数项 + 3个正弦余弦对)
p0 = [0] * 7
popt, pcov = curve_fit(fourier_series, x_data, y_noisy, p0=p0)
y_fourier = fourier_series(x_pred, *popt)
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(x_data, y_noisy, 'ro', label='噪声数据', markersize=4)
plt.plot(x_pred, y_spline, 'b-', label='样条插值')
plt.plot(x_pred, y_fourier, 'g--', label='傅里叶拟合')
plt.plot(x_pred, 2*np.sin(1.5*x_pred) + 0.5*np.cos(3*x_pred),
'k:', label='真实函数', linewidth=2)
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('曲线拟合比较')
plt.grid(True)
# 预测性能评估
plt.subplot(2, 1, 2)
x_extrapolate = np.linspace(2*np.pi, 3*np.pi, 100)
y_true_ext = 2 * np.sin(1.5 * x_extrapolate) + 0.5 * np.cos(3 * x_extrapolate)
# 外推预测
y_spline_ext = spline(x_extrapolate)
y_fourier_ext = fourier_series(x_extrapolate, *popt)
plt.plot(x_extrapolate, y_true_ext, 'k:', label='真实函数', linewidth=2)
plt.plot(x_extrapolate, y_spline_ext, 'b-', label='样条外推')
plt.plot(x_extrapolate, y_fourier_ext, 'g--', label='傅里叶外推')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('外推预测性能')
plt.grid(True)
plt.tight_layout()
plt.show()
# 计算拟合误差
mse_spline = np.mean((y_spline - (2*np.sin(1.5*x_pred) + 0.5*np.cos(3*x_pred)))**2)
mse_fourier = np.mean((y_fourier - (2*np.sin(1.5*x_pred) + 0.5*np.cos(3*x_pred)))**2)
print(f"样条插值MSE: {mse_spline:.6f}")
print(f"傅里叶拟合MSE: {mse_fourier:.6f}")
data_fitting_prediction()
9.3 金融计算：期权定价
python
def option_pricing():
"""
布莱克-舒尔斯期权定价模型数值解
"""
# 布莱克-舒尔斯偏微分方程
def black_scholes_pde(S, V, sigma, r):
"""
布莱克-舒尔斯PDE: ∂V/∂t + 0.5*σ²S²∂²V/∂S² + rS∂V/∂S - rV = 0
"""
# 使用中心差分近似导数
dS = S[1] - S[0]
# 一阶导数 (中心差分)
dV_dS = np.zeros_like(V)
dV_dS[1:-1] = (V[2:] - V[:-2]) / (2 * dS)
# 边界条件
dV_dS[0] = (V[1] - V[0]) / dS # 前向差分
dV_dS[-1] = (V[-1] - V[-2]) / dS # 后向差分
# 二阶导数
d2V_dS2 = np.zeros_like(V)
d2V_dS2[1:-1] = (V[2:] - 2*V[1:-1] + V[:-2]) / (dS**2)
# 边界条件
d2V_dS2[0] = d2V_dS2[1]
d2V_dS2[-1] = d2V_dS2[-2]
# PDE右侧
return -0.5 * sigma**2 * S**2 * d2V_dS2 - r * S * dV_dS + r * V
# 参数
S_max = 200 # 最大股价
N = 100 # 价格网格点数
M = 1000 # 时间步数
T = 1.0 # 到期时间(年)
K = 100 # 行权价
r = 0.05 # 无风险利率
sigma = 0.2 # 波动率
# 网格设置
S = np.linspace(0, S_max, N)
dt = T / M
# 初始条件 (欧式看涨期权到期 payoff)
V = np.maximum(S - K, 0)
# 时间反向演化 (从到期到当前)
for t in range(M):
# 显式欧拉方法
V = V + dt * black_scholes_pde(S, V, sigma, r)
# 边界条件
V[0] = 0 # S=0时期权价值为0
V[-1] = S_max - K * np.exp(-r * (T - t*dt)) # S很大时的边界条件
# 插值得到特定价格的期权价值
from scipy import interpolate
interp = interpolate.interp1d(S, V)
S0 = 100 # 当前股价
option_value = interp(S0)
print(f"布莱克-舒尔斯期权定价结果:")
print(f"当前股价: {S0}")
print(f"行权价: {K}")
print(f"到期时间: {T}年")
print(f"无风险利率: {r*100}%")
print(f"波动率: {sigma*100}%")
print(f"期权理论价格: {option_value:.4f}")
# 绘制期权价值曲线
plt.figure(figsize=(10, 6))
plt.plot(S, V)
plt.axvline(x=S0, color='r', linestyle='--', label=f'当前股价 = {S0}')
plt.xlabel('股价')
plt.ylabel('期权价值')
plt.title('布莱克-舒尔斯期权定价')
plt.legend()
plt.grid(True)
plt.show()
option_pricing()
10. 总结与展望
10.1 数值分析的关键要点
通过本文的讲解和实例，我们可以总结出数值分析的一些关键要点：
-
误差意识：数值计算必然存在误差，需要理解误差来源并控制误差传播
-
算法选择：不同问题需要选择适当的算法，考虑精度、效率和稳定性
-
问题适应性：了解问题的数学特性（如刚性、病态条件等）对算法选择至关重要
-
实现细节：算法的具体实现方式会显著影响数值结果的精度和计算效率
10.2 Python数值计算的未来发展方向
Python数值计算领域仍在快速发展，主要趋势包括：
| 方向 | 描述 | 相关工具 |
|---|---|---|
| 自动微分 | 自动计算导数，用于机器学习和优化 | JAX, Autograd, PyTorch |
| GPU加速 | 利用GPU进行大规模数值计算 | CuPy, PyTorch, TensorFlow |
| 符号计算 | 符号与数值计算的结合 | SymPy, SageMath |
| 分布式计算 | 大规模并行数值计算 | Dask, MPI for Python |
| 可解释性 | 数值结果的可视化和解释 | Matplotlib, Plotly, Bokeh |
10.3 学习建议与资源
对于想要深入学习数值分析的读者，建议：
-
理论基础：扎实的数学基础是理解数值算法的关键
-
编程实践：通过实现算法加深理解，而不仅仅是调用库函数
-
性能分析：学习使用性能分析工具优化代码
-
社区参与：参与开源项目和学习社区讨论
推荐资源：
-
书籍：《Numerical Recipes》、《Numerical Analysis》
-
在线课程：Coursera、edX上的数值分析课程
-
文档：NumPy、SciPy官方文档
-
社区：Stack Overflow、GitHub、Python科学计算社区
数值分析是一个既深且广的领域，本文仅涵盖了其中的基础内容和常用方法。希望这篇技术博客能为你的数值分析学习之旅提供坚实的起点和实用的参考。
扫一扫
2455
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
