Python 数值计算与数值分析基础

Python 数值计算与数值分析基础

示例演示

当涉及到Python数值计算和数值分析时，下面是20个示例，涵盖了一些常见的用法：

1.数值积分：

在 Python 中，你可以使用 scipy.integrate 模块中的 quad 函数来进行数值积分

import numpy as np
from scipy.integrate import quad

# 定义被积函数
f = lambda x: np.sin(x)

# 计算积分
result, error = quad(f, 0, np.pi)

# 输出结果
print("积分结果:", result)
print("误差估计:", error)

积分结果: 2.0
误差估计: 2.220446049250313e-14

2.数值微分：

import numpy as np

# 定义函数
f = lambda x: np.cos(x)

# 定义微分点和一个很小的h值
x0 = 0.5
h = 1e-5

# 计算数值微分（使用中心差分法）
numerical_derivative = (f(x0 + h) - f(x0 - h)) / (2 * h)

# 输出结果
print("数值导数结果:", numerical_derivative)

数值导数结果: -0.47942553859647846

3.非线性方程求根：

在 Python 中，你可以使用 scipy.optimize 模块的 fsolve 函数来求解方程。

import numpy as np
from scipy.optimize import fsolve

# 定义方程
def equation(x):
    return x**2 - 2

# 使用 fsolve 进行求解
solution = fsolve(equation, 1.5)

# 输出结果
print("方程的解:", solution[0])

方程的解: 1.4142135623730951

4.线性方程组求解：

在 Python 中，使用 NumPy 库可以方便地处理矩阵和线性代数运算。使用 np.linalg.solve(A, b) 函数求解线性方程组Ax=b。

import numpy as np

# 定义矩阵 A 和向量 b
A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])

# 求解线性方程 Ax = b
x = np.linalg.solve(A, b)

# 输出结果
print("解 x:", x)

解 x: [-4.   4.5]

5.曲线拟合：

在 Python 中，你可以使用 NumPy 和 SciPy 库来生成数据，并使用 Matplotlib 进行绘图。对于拟合功能，可以使用 scipy.optimize.curve_fit 来实现。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 定义 sine 函数
def sin_func(x, a, b, c, d):
    return a * np.sin(b * x + c) + d

# 生成数据
x = np.linspace(0, 2 * np.pi, 100)
y = np.sin(x) + np.random.rand(len(x)) * 0.2  # 添加随机噪声

# 拟合数据
initial_guess = [1, 1, 0, 0]  # 初始猜测参数
params, covariance = curve_fit(sin_func, x, y, p0=initial_guess)

# 生成拟合曲线
y_fit = sin_func(x, *params)

# 绘制结果
plt.scatter(x, y, label='Data', color='blue')
plt.plot(x, y_fit, label='Fitted Curve', color='red')
plt.title('Sine Fit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

6.矩阵特征值和特征向量：

在 Python 中，你可以使用 NumPy 库来计算矩阵的特征值和特征向量。使用 np.linalg.eig 函数计算特征值和特征向量。

import numpy as np

# 定义矩阵 A
A = np.array([[1, 2], [3, 4]])

# 计算特征值和特征向量
eig_val, eig_vec = np.linalg.eig(A)

# 输出结果
print("特征值:\n", eig_val)
print("特征向量:\n", eig_vec)

特征值:
 [-0.37228132  5.37228132]
特征向量:
 [[-0.82456484 -0.41597356]
 [ 0.56576746 -0.90937671]]

7.傅里叶变换：

import numpy as np
import matplotlib.pyplot as plt

# 设置采样频率
Fs = 1000
# 生成时间序列
t = np.arange(0, 1, 1/Fs)
# 生成信号
x = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t)
# 进行快速傅里叶变换
y = np.fft.fft(x)
# 计算频率轴
f = np.fft.fftfreq(len(y), d=1/Fs)

# 绘制幅度谱
plt.plot(f, np.abs(y))
plt.title('FFT of the Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, Fs/2)  # 只显示正频率部分
plt.grid()
plt.show()

​

8.信号滤波：

以下代码创建了一个包含 50 Hz 和 120 Hz 正弦波的信号，然后使用带通滤波器对信号进行滤波，并将原始信号和滤波后的信号绘制在同一图中。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

# 设置采样频率
Fs = 1000
# 生成时间序列
t = np.arange(0, 1, 1/Fs)
# 生成信号
x = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t)

# 设计带通滤波器
lowcut = 45  # 低截止频率
highcut = 55  # 高截止频率
b, a = butter(3, [lowcut, highcut], fs=Fs, btype='bandpass')

# 使用 filtfilt 进行零相位滤波
filtered_signal = filtfilt(b, a, x)

# 绘制原始信号和滤波后的信号
plt.figure(figsize=(10, 6))
plt.plot(t, x, label='Original Signal', alpha=0.7)
plt.plot(t, filtered_signal, label='Filtered Signal', linewidth=2)
plt.title('Signal Before and After Bandpass Filtering')
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.legend()
plt.grid()
plt.show()

9.最小二乘拟合：

import numpy as np
import matplotlib.pyplot as plt

# 定义数据点
x = np.array([1, 2, 3, 4, 5])
y = np.array([1, 3, 6, 10, 15])

# 进行二次多项式拟合
p = np.polyfit(x, y, 2)

# 计算拟合值
f = np.polyval(p, x)

# 绘制数据点和拟合曲线
plt.plot(x, y, 'o', label='Data Points')  # 原始数据点
plt.plot(x, f, label='Fitted Curve')      # 拟合曲线
plt.title('Polynomial Fit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

10.数值优化问题求解：

使用 scipy.optimize 库中的 minimize 函数来实现无约束优化。

import numpy as np
from scipy.optimize import minimize

# 定义目标函数
def fun(x):
    return x[0]**2 + x[1]**2

# 初始猜测
x0 = np.array([1, 1])

# 使用 minimize 进行无约束优化
result = minimize(fun, x0)

# 输出优化结果
print("Optimized parameters:", result.x)
print("Objective function value at optimized parameters:", result.fun)

Optimized parameters: [-1.07505143e-08 -1.07505143e-08]
Objective function value at optimized parameters: 2.311471135620994e-16

11.数值积分：

使用 scipy.integrate 模块中的 quad 函数来进行数值积分。

import numpy as np
from scipy.integrate import quad

# 定义被积函数
def func(x):
    return 1 / (1 + x**2)

# 进行定积分
integral_value, error = quad(func, 0, 1)

# 输出积分值和误差
print("Integral value:", integral_value)
print("Estimated error:", error)

Integral value: 0.7853981633974484
Estimated error: 8.719671245021581e-15

12.插值：

使用 numpy 处理数组，使用 matplotlib.pyplot 绘制图形，使用 scipy.interpolate 进行插值。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

# 原始数据点
x = [0, 1, 2, 3]
y = [1, 4, 9, 16]

# 生成插值点
xi = np.arange(0, 3.1, 0.1)

# 创建插值函数
interp_func = interp1d(x, y, kind='linear')

# 计算插值结果
yi = interp_func(xi)

# 绘制原始数据点和插值结果
plt.plot(x, y, 'o', label='Original data')
plt.plot(xi, yi, label='Interpolated data')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Interpolation Example')
plt.legend()
plt.grid()
plt.show()

13.求解常微分方程组：

使用 scipy.integrate 模块中的 solve_ivp 函数来解决常微分方程，并使用 matplotlib 来绘制图形。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定义微分方程
def func(t, y):
    return [y[1], -y[0]]

# 设定时间范围和初始条件
t_span = (0, 10)
y0 = [1, 0]

# 使用 solve_ivp 解决微分方程
sol = solve_ivp(func, t_span, y0, t_eval=np.linspace(0, 10, 100))

# 绘制结果
plt.plot(sol.t, sol.y[0])
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.title('Solution of ODE')
plt.grid()
plt.show()

14.数值统计：

在 Python 中，numpy 库提供了计算均值和标准差的函数。

import numpy as np

# 定义数据
data = [1, 2, 3, 4, 5]

# 计算均值和标准差
mean_value = np.mean(data)
std_value = np.std(data)

# 打印结果
print("Mean:", mean_value)
print("Standard Deviation:", std_value)

Mean: 3.0
Standard Deviation: 1.4142135623730951

15.随机数生成：

在 Python 中，可以使用 numpy 库的 random.rand 函数生成随机数。

import numpy as np

# 生成 1 到 10 之间的随机数
random_num = np.random.rand(1, 10)

# 打印结果
print(random_num)

[[0.08318904 0.2923153  0.1923046  0.96278395 0.93616543 0.48201543
  0.27400168 0.62454063 0.63461248 0.64493564]]

16.多项式求根：

在 Python 中，可以使用 numpy 库的 numpy.roots 函数来计算多项式的根。

import numpy as np

# 定义多项式系数
coefficients = [1, -3, 2]

# 计算多项式的根
roots_of_polynomial = np.roots(coefficients)

# 打印结果
print("Roots of the polynomial:", roots_of_polynomial)

Roots of the polynomial: [2. 1.]

17.矩阵求逆：

在 Python 中，可以使用 numpy 库来创建矩阵并计算其逆。使用 np.linalg.inv(A) 计算矩阵A的逆。

import numpy as np

# 定义矩阵 A
A = np.array([[1, 2], [3, 4]])

# 计算矩阵 A 的逆
inv_A = np.linalg.inv(A)

# 打印结果
print("Inverse of A:\n", inv_A)

Inverse of A:
 [[-2.   1. ]
 [ 1.5 -0.5]]

18.线性插值：

import numpy as np
import matplotlib.pyplot as plt

# 定义数据点
x = [0, 1, 2]
y = [1, 3, 2]

# 定义插值点
xi = np.arange(0, 2.1, 0.1)  # 生成从0到2（包括2）的数，步长为0.1

# 进行线性插值
yi = np.interp(xi, x, y)

# 绘制结果
plt.plot(x, y, 'o', label='Data Points')  # 绘制原始数据点
plt.plot(xi, yi, label='Linear Interpolation')  # 绘制插值曲线
plt.xlabel('x')
plt.ylabel('y')
plt.title('Linear Interpolation Example')
plt.legend()
plt.grid()
plt.show()

19.蒙特卡洛方法：

以下Python 代码，使用 numpy 来生成随机数，并估算圆周率 π 的值。

import numpy as np

# 初始化计数器和样本数
count = 0
n = 100000

# 进行 Monte Carlo 模拟
for _ in range(n):
    x = np.random.rand()  # 生成 [0, 1) 之间的随机数
    y = np.random.rand()  # 生成 [0, 1) 之间的随机数
    if x**2 + y**2 <= 1:  # 检查是否在单位圆内
        count += 1

# 估算 π
pi_estimate = 4 * count / n
print(f"Estimated value of π: {pi_estimate}")

Estimated value of π: 3.14684

20.矩阵求秩：

在 Python 中，可以使用 numpy 库来创建矩阵并计算其秩（rank）。

import numpy as np

# 定义矩阵 A
A = np.array([[1, 2], [3, 4]])

# 计算矩阵 A 的秩
rank_A = np.linalg.matrix_rank(A)

print(f"Rank of matrix A: {rank_A}")

Rank of matrix A: 2

这些示例展示了Python中数值计算和数值分析的一些常见用法。你可以通过这些示例来了解如何使用Python进行数值计算和分析，并可根据具体需求进行进一步修改和调整。在实际应用中，可以根据具体问题选择合适的函数和方法进行数值计算和分析。