哈达玛矩阵
H = hadamard(n) 返回阶次为n的哈达玛矩阵.
n — 矩阵的阶次,非负整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def hadamard_matrix(input_str):
"""
构造Hadamard矩阵，对标MATLAB的hadamard函数
参数:
input_str (str): 输入字符串，应为正整数n（需满足n是2的幂）
返回:
Sympy矩阵或错误信息
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def _recursive_hadamard(n):
"""递归构造2的幂次阶Hadamard矩阵"""
if n == 1:
return sp.Matrix([[1]])
H_prev = _recursive_hadamard(n // 2)
# 矩阵拼接：行拼接和列拼接
top = H_prev.row_join(H_prev)
bottom = H_prev.row_join(-H_prev)
return top.col_join(bottom)
# 验证输入合法性
if expr.is_Integer and expr > 0:
n = int(expr)
# 检查是否为2的幂
if (n & (n - 1)) == 0 or n == 1:
result = _recursive_hadamard(n)
else:
error = True
else:
error = True
if error:
return f"输入错误: {input_str}\n提示: n必须是2的幂的正整数"
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：有效输入 n=4
input1 = '4'
mat1 = hadamard_matrix(input1)
print("示例1结果(n=4):\n", mat1)
# Matrix([[1, 1, 1, 1],
# [1, -1, 1, -1],
# [1, 1, -1, -1],
# [1, -1, -1, 1]])
# 示例2：边界输入 n=1
input2 = '1'
mat2 = hadamard_matrix(input2)
print("示例2结果(n=1):\n", mat2)
# Matrix([[1]])
半波正弦波信号
HalfRectSineWave(A)
A - 数组,或符号表达式
import matplotlib.pyplot as plt
import numpy as np
示例1：模拟电源整流
50Hz交流电整流（市电频率）
t_power, y_power = HalfRectSineWave(t, Fs=1e3, f=50, num_cycles=3)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t_power * 1000, y_power) # 转换为毫秒
plt.title("50Hz半波整流电源输出")
plt.xlabel("时间 (毫秒)")
plt.ylabel("电压 (V)")
plt.grid(True)
添加电容滤波效果（简单模拟）
t_filtered, y_raw = HalfRectSineWave(t, Fs=10e3, f=50, num_cycles=2)
# 简单RC滤波模拟
alpha = 0.05 # 滤波系数
y_filtered = np.zeros_like(y_raw)
y_filtered[0] = y_raw[0]
for i in range(1, len(y_raw)):
y_filtered[i] = alpha * y_raw[i] + (1 - alpha) * y_filtered[i - 1]
plt.subplot(1, 2, 2)
plt.plot(t_filtered * 1000, y_raw, 'b-', alpha=0.5, label='整流后')
plt.plot(t_filtered * 1000, y_filtered, 'r-', linewidth=2, label='滤波后')
plt.title("整流信号与滤波效果")
plt.xlabel("时间 (毫秒)")
plt.ylabel("电压 (V)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
示例2. 音频信号处理
模拟音频信号包络检测
t_audio, y_audio = HalfRectSineWave(t, Fs=44.1e3, f=1000, num_cycles=10)
添加高频载波和低频调制
carrier_freq = 10e3 # 10kHz载波
modulation_freq = 500 # 500Hz调制
AM调制信号
am_signal = (1 + 0.5 * np.sin(2 * np.pi * modulation_freq * t_audio)) * \
np.sin(2 * np.pi * carrier_freq * t_audio) * 100
包络检测（半波整流）
t_env, env_detected = HalfRectSineWave(list(am_signal), Fs=44.1e3, f=carrier_freq)
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t_audio[:2000] * 1000, am_signal[:2000])
plt.title("AM调制信号")
plt.xlabel("时间 (毫秒)")
plt.ylabel("幅度")
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(t_env[:2000] * 1000, env_detected[:2000])
plt.title("包络检测输出")
plt.xlabel("时间 (毫秒)")
plt.ylabel("幅度")
plt.grid(True)
plt.tight_layout()
plt.show()
示例3. 机械振动分析
使用时间变换来模拟非线性响应
t_vib, y_vib = HalfRectSineWave(t * (1 + 0.1*t), Fs=5e3, f=100, num_cycles=8)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t_vib * 1000, y_vib)
plt.title("机械振动冲击响应")
plt.xlabel("时间 (毫秒)")
plt.ylabel("加速度 (m/s²)")
plt.grid(True)
分析不同频率的振动
frequencies = [50, 100, 200]
plt.subplot(1, 2, 2)
for freq in frequencies:
t_temp, y_temp = HalfRectSineWave(t, Fs=5e3, f=freq, num_cycles=4)
plt.plot(t_temp * 1000, y_temp, label=f'{freq}Hz')
plt.title("不同频率振动响应")
plt.xlabel("时间 (毫秒)")
plt.ylabel("加速度 (m/s²)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
示例4. 生物医学信号模拟
使用自定义时间点模拟不规则心跳
irregular_beats = [0, 0.8, 1.9, 3.0, 4.2, 5.1] # 不规则心跳时间点（秒）
生成心电信号中的R波（半波整流正弦波模拟）
t_ecg = np.linspace(0, 6, 6000) # 6秒信号
ecg_signal = np.zeros_like(t_ecg)
for beat_time in irregular_beats:
# 在每个心跳时间生成R波
beat_duration = 0.1 # R波持续时间
beat_samples = int(beat_duration * 1000) # 采样点数
t_beat = np.linspace(0, beat_duration, beat_samples)
# 生成R波（使用半波整流正弦波形状）
_, r_wave = HalfRectSineWave(list(t_beat), f=15, Fs=1000)
# 将R波插入到ECG信号中
start_idx = int(beat_time * 1000)
end_idx = start_idx + len(r_wave)
if end_idx < len(ecg_signal):
ecg_signal[start_idx:end_idx] += r_wave
plt.figure(figsize=(12, 4))
plt.plot(t_ecg, ecg_signal)
plt.title("模拟心电图信号（R波检测）")
plt.xlabel("时间 (秒)")
plt.ylabel("幅度")
plt.grid(True)
标记R波位置
for beat_time in irregular_beats:
plt.axvline(x=beat_time, color='red', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
示例5. 通信系统同步信号
生成时钟同步信号
t_sync, y_sync = HalfRectSineWave(t, Fs=100e3, f=10e3, num_cycles=20)
添加相位调制来模拟时钟抖动
t_jitter, y_jitter = HalfRectSineWave(t + 0.00002*sin(2*pi*1000*t),
Fs=100e3, f=10e3, num_cycles=20)
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t_sync * 1e6, y_sync) # 转换为微秒
plt.title("理想时钟同步信号 (10kHz)")
plt.xlabel("时间 (微秒)")
plt.ylabel("幅度")
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(t_jitter * 1e6, y_jitter)
plt.title("带时钟抖动的同步信号")
plt.xlabel("时间 (微秒)")
plt.ylabel("幅度")
plt.grid(True)
plt.tight_layout()
计算时钟误差
ideal_peaks = np.where(y_sync > 95)[0]
jitter_peaks = np.where(y_jitter > 95)[0]
if len(ideal_peaks) == len(jitter_peaks):
timing_errors = (jitter_peaks - ideal_peaks) / 100e3 * 1e6 # 转换为微秒
print(f"最大时钟抖动: {np.max(np.abs(timing_errors)):.2f} 微秒")
print(f"平均时钟抖动: {np.mean(np.abs(timing_errors)):.2f} 微秒")
plt.show()
示例6. 功率电子器件分析
分析不同导通角的整流效果
conduction_angles = [30, 60, 90, 120] # 度
plt.figure(figsize=(12, 8))
for i, angle_deg in enumerate(conduction_angles):
angle_rad = np.radians(angle_deg)
phase_shift = angle_rad / (2 * np.pi * 1000) # 转换为时间偏移
# 使用相位偏移模拟可控硅的导通角控制
t_scr, y_scr = HalfRectSineWave(t + phase_shift,
Fs=20e3, f=1000, num_cycles=3)
plt.subplot(2, 2, i + 1)
plt.plot(t_scr * 1000, y_scr)
plt.title(f"导通角 {angle_deg}° 的整流输出")
plt.xlabel("时间 (毫秒)")
plt.ylabel("输出电压")
plt.grid(True)
# 计算平均输出电压
avg_voltage = np.mean(y_scr)
plt.text(0.5, 80, f'平均电压: {avg_voltage:.1f}V',
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8))
plt.tight_layout()
plt.show()
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import ast
def half_rect_sine_wave(input_str, Fs=20e3, f=1e3, num_cycles=6):
"""
生成半波整流正弦波。
参数:
input_str (str): 输入表达式或时间列表字符串。
- 若为列表字符串 (如 "[0, 1e-3, 2e-3]"), 表示自定义时间点（秒）。
- 若为SymPy表达式 (如 "t", "t*0.5"), 表示时间变换公式。
Fs (float): 采样率 (Hz)，默认20kHz。
f (float): 信号频率 (Hz)，默认1kHz。
num_cycles (int): 生成信号的周期数（仅当输入为表达式时有效），默认6个周期。
返回:
tuple: (时间数组, 半波整流后的信号数组) 或错误信息字符串。
"""
try:
# 解析输入字符串
try:
# 尝试解析为Python列表
expr = ast.literal_eval(input_str)
if not isinstance(expr, list):
# 不是列表则尝试SymPy解析
expr = sp.sympify(input_str)
except (ValueError, SyntaxError, TypeError):
# 解析失败则作为SymPy表达式处理
expr = sp.sympify(input_str)
# 处理列表输入（自定义时间点）
if isinstance(expr, list):
time_values = np.array(expr, dtype=float)
# 生成正弦波：y = 100*sin(2πf*t)
y = 100 * np.sin(2 * np.pi * f * time_values)
y_rectified = np.maximum(y, 0)
return time_values, y_rectified
# 处理SymPy表达式
elif isinstance(expr, sp.Expr):
# 计算总样本数
samples_per_cycle = int(Fs / f)
total_samples = num_cycles * samples_per_cycle
# 生成时间数组（秒）
time_base = np.arange(total_samples) / Fs
# 将表达式转换为函数
variables = expr.free_symbols
if variables:
func = sp.lambdify(variables, expr, modules='numpy')
# 代入时间计算处理后的时间值
t_processed = func(time_base)
else:
# 无变量则视为常数表达式
t_processed = np.full_like(time_base, float(expr))
# 生成正弦波并整流
y = 100 * np.sin(2 * np.pi * f * t_processed)
y_rectified = np.maximum(y, 0)
return time_base, y_rectified
else:
raise ValueError("输入类型不支持，请输入列表或SymPy表达式。")
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
import matplotlib.pyplot as plt
# 示例1：默认参数生成6个周期的波形
t1, y1 = half_rect_sine_wave("t")
plt.figure()
plt.plot(t1, y1)
plt.title("默认参数半波整流正弦波")
plt.xlabel("时间 (秒)")
plt.ylabel("幅度")
# 示例2：自定义时间点
custom_time = [0, 0.0005, 0.001, 0.0015] # 0, 0.5ms, 1ms, 1.5ms
t2, y2 = half_rect_sine_wave(str(custom_time))
plt.figure()
plt.stem(t2, y2, use_line_collection=True)
plt.title("自定义时间点半波整流")
plt.xlabel("时间 (秒)")
plt.ylabel("幅度")
# 示例3：带相位偏移的波形
t3, y3 = half_rect_sine_wave("t + 0.00025", num_cycles=3) # 0.25ms相位偏移
plt.figure()
plt.plot(t3, y3)
plt.title("带相位偏移的半波整流")
plt.xlabel("时间 (秒)")
plt.ylabel("幅度")
plt.show()
汉明窗
w = hamming(L,sym=1) 返回一个长度为L个点的修正的汉明窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
示例1. 数字滤波器设计
"""设计一个低通FIR滤波器"""
使用Hamming窗设计51阶低通FIR滤波器
N = 51
cutoff_freq = 0.3 # 归一化截止频率
生成理想低通滤波器的冲激响应
h_ideal = []
for n in range(N):
if n == (N - 1) // 2:
h_ideal.append(2 * cutoff_freq)
else:
theta = 2 * np.pi * cutoff_freq * (n - (N - 1) / 2)
h_ideal.append(np.sin(theta) / (np.pi * (n - (N - 1) / 2)))
应用Hamming窗
window = hamming(N)
h_fir = [h * w for h, w in zip(h_ideal, window)]
print(f"低通FIR滤波器系数 (前10个): {h_fir[:10]}")
#低通FIR滤波器系数 (前10个): [5.489886917185572e-18, 0.0010548554750528494, -0.0007683356183186414,
-0.0009550734669726426, 0.0019735003278461144, -3.925899901137517e-18,
-0.003261117480891942, 0.0025651632210178265, 0.0032304059987978613, -0.006511387630530373]
示例2. 频谱分析 - 语音信号处理
"""语音信号的频谱分析"""
模拟语音信号片段
fs = 8000 # 采样率8kHz
t = np.linspace(0, 0.1, 800) # 100ms信号
生成包含基频和泛音的语音-like信号
f0 = 200 # 基频200Hz
signal_data = (np.sin(2 * np.pi * f0 * t) +
0.5 * np.sin(2 * np.pi * 2 * f0 * t) +
0.3 * np.sin(2 * np.pi * 3 * f0 * t) +
0.1 * np.random.normal(0, 0.1, len(t)))
应用Hamming窗进行频谱分析
N = 256
window = hamming(N)
signal_windowed = signal_data[:N] * window
计算FFT
spectrum = np.fft.fft(signal_windowed)
freq = np.fft.fftfreq(N, 1 / fs)
print(f"频谱峰值频率: {freq[np.argmax(np.abs(spectrum[:N // 2]))]:.1f} Hz")
#频谱峰值频率: 187.5 Hz
示例3. 雷达信号处理
"""雷达脉冲压缩（匹配滤波）"""
生成线性调频信号（LFM）
T = 10e-6 # 脉冲宽度10微秒
B = 10e6 # 带宽10MHz
fs = 50e6 # 采样率50MHz
t = np.linspace(-T / 2, T / 2, int(T * fs))
chirp_signal = np.exp(1j * np.pi * (B / T) * t ** 2) # 复数LFM信号
应用Hamming窗减少距离旁瓣
N = len(chirp_signal)
window = hamming(N)
weighted_chirp = chirp_signal * window
匹配滤波器（时域反转共轭）
matched_filter = np.conj(weighted_chirp[::-1])
print(f"加窗后脉冲长度: {len(weighted_chirp)} 点")
print(f"主瓣宽度改善: 使用Hamming窗可降低距离旁瓣约-43dB")
#加窗后脉冲长度: 500 点
#主瓣宽度改善: 使用Hamming窗可降低距离旁瓣约-43dB
示例4. 图像处理 - 频域滤波
"""在频域中使用Hamming窗进行图像滤波"""
创建测试图像（包含不同频率成分）
M, N = 64, 64
x, y = np.meshgrid(np.linspace(-1, 1, N), np.linspace(-1, 1, M))
生成包含高低频的图像
test_image = (np.sin(10 * np.pi * x) + # 高频成分
np.cos(4 * np.pi * y) + # 中频成分
0.5 * np.ones_like(x)) # 直流分量
创建2D Hamming窗
window_1d = hamming(M)
window_2d = np.outer(window_1d, window_1d)
应用窗函数后进行FFT
image_windowed = test_image * window_2d
image_fft = np.fft.fft2(image_windowed)
print(f"加窗前图像能量: {np.sum(test_image ** 2):.3f}")
print(f"加窗后图像能量: {np.sum(image_windowed ** 2):.3f}")
#加窗前图像能量: 5184.000
#加窗后图像能量: 784.075
#Hamming窗有效减少了频谱泄漏
示例5. 通信系统 - OFDM子载波整形
"""OFDM系统中子载波的频谱整形"""
OFDM参数
N_subcarriers = 64 # 子载波数量
CP_length = 16 # 循环前缀长度
symbol_length = N_subcarriers + CP_length
生成随机QAM符号
qam_symbols = np.random.randint(0, 4, N_subcarriers) # 4-QAM
qam_constellation = [1 + 1j, 1 - 1j, -1 + 1j, -1 - 1j]
mapped_symbols = [qam_constellation[s] for s in qam_symbols]
IFFT生成时域信号
ifft_output = np.fft.ifft(mapped_symbols, N_subcarriers)
应用Hamming窗进行时域加窗
window = hamming(symbol_length, 0) # 周期窗
windowed_symbol = np.concatenate([ifft_output[-CP_length:], ifft_output]) * window
print(f"OFDM符号长度: {len(windowed_symbol)}")
#OFDM符号长度: 80
#加窗后PAPR(峰均比)降低，频谱带外衰减更好
示例6. 生物医学信号处理
"""心电信号(ECG)的频谱特征提取"""
模拟ECG信号片段
fs = 360 # 采样率360Hz
t = np.linspace(0, 2, 720) # 2秒信号
模拟ECG特征：P波、QRS复合波、T波
ecg_signal = (0.5 * np.exp(-((t - 0.2) / 0.03) ** 2) + # P波
2.0 * np.exp(-((t - 0.3) / 0.02) ** 2) - # QRS波
0.8 * np.exp(-((t - 0.5) / 0.06) ** 2)) # T波
添加噪声
noise = 0.1 * np.random.normal(0, 1, len(t))
noisy_ecg = ecg_signal + noise
使用Hamming窗进行频谱分析
segment_length = 256
window = hamming(segment_length)
对信号分段加窗分析
segments = []
for i in range(0, len(noisy_ecg) - segment_length, segment_length // 2):
segment = noisy_ecg[i:i + segment_length] * window
segments.append(segment)
print(f"ECG信号分段数: {len(segments)}")
#ECG信号分段数: 4
示例7. 声学测量 - 房间脉冲响应
"""房间声学特性的测量分析"""
模拟房间脉冲响应
N = 512
t = np.linspace(0, 0.1, N) # 100ms响应
模拟包含直达声和早期反射的脉冲响应
impulse_response = np.zeros(N)
impulse_response[10] = 1.0 # 直达声
impulse_response[45] = 0.6 # 第一次反射
impulse_response[78] = 0.4 # 第二次反射
impulse_response[120] = 0.3 # 第三次反射
添加混响尾音（指数衰减）
for i in range(150, N):
impulse_response[i] = 0.8 * np.exp(-0.05 * (i - 150))
应用Hamming窗提取早期反射
window = hamming(200) # 关注前200个采样点
window_full = np.concatenate([window, np.zeros(N - 200)])
windowed_response = impulse_response * window_full
计算混响时间(RT60)等参数
energy = np.cumsum(windowed_response[::-1] ** 2)[::-1]
#Hamming窗减少了频谱泄漏，提高心率检测准确性
#使用Hamming窗分离早期反射和混响尾音
#便于分析房间的早期衰减特性
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.signal as signal
def hamming_window(input_str):
"""
生成Hamming窗，对标MATLAB的hamming函数。
参数:
input_str (str): 输入字符串，支持以下格式：
- 单个数字（如"64"）生成对称窗（'symmetric'）。
- 元组（如"(64, 0)"），第二个元素为0生成周期窗（'periodic'），非零生成对称窗。
返回:
list or str: Hamming窗列表，或错误信息。
"""
try:
expr = sp.sympify(input_str)
result = None
# 处理元组输入（例如 "(m, flag)"）
if isinstance(expr, tuple) and len(expr) == 2:
m, flag = expr
if m.is_number and flag.is_number:
m_int = int(m)
# flag=0 对应周期窗（sym=False），否则对称窗（sym=True）
sym = bool(int(flag) != 0)
window = signal.windows.hamming(m_int, sym)
result = list(window)
else:
return f"输入错误: 参数必须为数字 {input_str}"
# 处理单个数字输入（例如 "m"）
elif expr.is_number:
m_int = int(expr)
window = signal.windows.hamming(m_int, sym=True) # 默认对称窗
result = list(window)
else:
return f"输入格式错误: {input_str}"
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 生成对称窗（64点，MATLAB默认）
print("对称窗（64点）:", hamming_window("64")[:5])
# [0.08000000000000007, 0.08228584333191541, 0.08912065558966992, 0.10043650933837528, 0.11612094253961253]
# 生成周期窗（64点，第二个参数为0）
print("周期窗（64点）:", hamming_window("64, 0")[:5])
# [0.08000000000000007, 0.08221502573078954, 0.0888387710145141, 0.099807445563184, 0.11501541504480817]
汉宁窗
w = hann(L,sym=1) 返回一个长度为L个点的修正的汉明窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
示例1. 音频信号处理 - 音乐频谱分析
模拟音乐信号（包含和弦）
fs = 44100 # CD音质采样率
duration = 0.1 # 100ms分析窗口
t = np.linspace(0, duration, int(fs * duration))
C大调和弦：C(261.63Hz), E(329.63Hz), G(392.00Hz)
chord_signal = (0.6 * np.sin(2 * np.pi * 261.63 * t) +
0.4 * np.sin(2 * np.pi * 329.63 * t) +
0.3 * np.sin(2 * np.pi * 392.00 * t) +
0.1 * np.random.normal(0, 0.05, len(t)))
应用Hann窗进行频谱分析
N = len(t)
hann_win = hann(N)
windowed_signal = chord_signal * hann_win
FFT分析
spectrum = np.fft.fft(windowed_signal)
freqs = np.fft.fftfreq(N, 1 / fs)
查找峰值频率
magnitude = 20 * np.log10(np.abs(spectrum[:N // 2]) + 1e-10)
peak_freqs = []
for i in range(1, len(magnitude) - 1):
if magnitude[i] > magnitude[i - 1] and magnitude[i] > magnitude[i + 1]:
if magnitude[i] > -40: # 高于-40dB的峰值
peak_freqs.append(freqs[i])
检测到的和弦频率成分:
for freq in peak_freqs[:6]: # 显示前6个峰值
print(f" {freq:.1f} Hz")
#50.0 Hz
#120.0 Hz
#140.0 Hz
#160.0 Hz
#190.0 Hz
#260.0 Hz
示例2. 结构健康监测 - 振动模态分析
模拟桥梁振动信号
fs = 1000 # 采样率1kHz
t = np.linspace(0, 10, 10000) # 10秒数据
结构固有频率：1.5Hz, 4.2Hz, 8.7Hz (典型桥梁频率)
vibration = (2.0 * np.sin(2 * np.pi * 1.5 * t) * np.exp(-0.01 * t) + # 一阶模态，衰减
1.2 * np.sin(2 * np.pi * 4.2 * t) * np.exp(-0.02 * t) + # 二阶模态
0.8 * np.sin(2 * np.pi * 8.7 * t) * np.exp(-0.03 * t) + # 三阶模态
0.3 * np.random.normal(0, 0.2, len(t))) # 环境噪声
分段Hann窗分析（短时傅里叶变换）
segment_length = 2048
overlap = 1024
hann_win = hann(segment_length)
spectrogram = []
time_points = []
for i in range(0, len(vibration) - segment_length, segment_length - overlap):
segment = vibration[i:i + segment_length] * hann_win
seg_fft = np.fft.fft(segment)
spectrogram.append(20 * np.log10(np.abs(seg_fft[:segment_length // 2]) + 1e-10))
time_points.append(t[i])
结构振动模态分析完成:
print(f"使用{segment_length}点Hann窗，重叠{overlap}点")
print("Hann窗提供优良的频率分辨率和较低的频谱泄漏")
#使用2048点Hann窗，重叠1024点
#Hann窗提供优良的频率分辨率和较低的频谱泄漏
示例3. 雷达信号处理 - 多目标分辨
雷达参数
fs = 100e6 # 100MHz采样率
T = 20e-6 # 20微秒脉冲宽度
t = np.linspace(0, T, int(fs * T))
模拟三个不同距离的目标回波
target_delays = [5e-6, 8e-6, 8.5e-6] # 微秒延迟
target_amplitudes = [1.0, 0.8, 0.6] # 目标强度
radar_echo = np.zeros(len(t))
for delay, amp in zip(target_delays, target_amplitudes):
idx = int(delay * fs)
if idx < len(t):
radar_echo[idx] = amp
添加噪声和展宽效应
radar_echo = np.convolve(radar_echo, np.exp(-(t - 1e-6) ** 2 / (2e-6) ** 2), 'same')
radar_echo += 0.1 * np.random.normal(0, 0.05, len(t))
应用Hann窗提高距离分辨力
hann_win = hann(len(t))
windowed_echo = radar_echo * hann_win
距离谱分析
range_spectrum = np.fft.fft(windowed_echo)
range_bins = np.fft.fftfreq(len(t), 1 / fs) * 3e8 / 2 # 转换为距离
目标距离检测结果:
peaks, _ = signal.find_peaks(np.abs(range_spectrum[:len(t) // 2]), height=0.1)
for peak in peaks[:3]:
distance = range_bins[peak]
if distance > 0:
print(f" 目标距离: {distance:.1f} 米")
# 目标距离: 75000000000000.0 米
# 目标距离: 97500000000000.0 米
# 目标距离: 120000000000000.0 米
示例4. 医学影像 - 超声成像
模拟超声回波信号（A模式）
fs = 50e6 # 50MHz超声采样率
depth = 0.15 # 15cm探测深度
t = np.linspace(0, depth * 2 / 1540, int(fs * depth * 2 / 1540)) # 声速1540m/s
模拟不同组织的反射界面
ultrasound_signal = np.zeros(len(t))
# 组织界面位置（米）
interfaces = [0.02, 0.05, 0.08, 0.12]
for interface in interfaces:
idx = int(interface * 2 / 1540 * fs)
if idx < len(t):
ultrasound_signal[idx] = 0.8 * np.exp(-interface * 10) # 深度衰减
添加超声系统的冲激响应
impulse_response = np.exp(-(np.arange(-10, 11) * 1e-7) ** 2 / (2e-7) ** 2)
ultrasound_signal = np.convolve(ultrasound_signal, impulse_response, 'same')
加噪声
ultrasound_signal += 0.1 * np.random.normal(0, 0.05, len(t))
应用Hann窗进行包络检测
hann_win = hann(len(t))
windowed_signal = ultrasound_signal * hann_win
Hilbert变换提取包络
analytic_signal = signal.hilbert(windowed_signal)
envelope = np.abs(analytic_signal)
超声成像处理完成:
print("Hann窗有效抑制旁瓣，提高纵向分辨力")
print(f"检测到 {len(interfaces)} 个组织界面")
#Hann窗有效抑制旁瓣，提高纵向分辨力
#检测到 4 个组织界面
示例5. 通信系统 - 5G NR波形生成
5G NR参数
scs = 30e3 # 30kHz子载波间隔
fft_size = 4096
num_symbols = 14 # 一个时隙的OFDM符号数
生成5G NR上行参考信号
reference_signal = np.random.randn(fft_size) + 1j * np.random.randn(fft_size)
reference_signal = reference_signal / np.sqrt(2) # 功率归一化
时域OFDM符号生成
ofdm_symbols = []
for i in range(num_symbols):
# IFFT
time_domain = np.fft.ifft(reference_signal, fft_size)
# 应用Hann窗进行频谱整形（降低带外辐射）
hann_win = hann(len(time_domain), 1) # 对称窗
windowed_symbol = time_domain * hann_win
ofdm_symbols.append(windowed_symbol)
计算功率谱密度
full_signal = np.concatenate(ofdm_symbols)
f, psd = signal.welch(full_signal, fs=fft_size * scs, window='hann',
nperseg=1024, return_onesided=False)
5G NR波形生成完成:
print(f"使用{fft_size}点FFT，{scs / 1000}kHz子载波间隔")
print("Hann窗有效降低带外辐射，满足3GPP频谱模板要求")
#使用4096点FFT，30.0kHz子载波间隔
#Hann窗有效降低带外辐射，满足3GPP频谱模板要求
示例6. 地球物理勘探 - 地震信号处理
模拟地震检波器信号
fs = 500 # 500Hz采样率
duration = 2.0 # 2秒记录
t = np.linspace(0, duration, int(fs * duration))
模拟地震事件：P波、S波、面波
seismic_signal = np.zeros(len(t))
P波（纵波，速度快）
p_wave_time = 0.8
p_wave_idx = int(p_wave_time * fs)
if p_wave_idx < len(t):
seismic_signal[p_wave_idx:p_wave_idx + 100] += 1.0 * np.exp(-(np.arange(100) - 50) ** 2 / 100)
S波（横波，速度慢）
s_wave_time = 1.2
s_wave_idx = int(s_wave_time * fs)
if s_wave_idx < len(t):
seismic_signal[s_wave_idx:s_wave_idx + 150] += 0.7 * np.exp(-(np.arange(150) - 75) ** 2 / 200)
面波（低频，幅度大）
surface_wave_time = 1.5
surface_idx = int(surface_wave_time * fs)
if surface_idx < len(t):
surface_wave = 2.0 * np.sin(2 * np.pi * 2 * (t[surface_idx:] - surface_wave_time)) * \
np.exp(-5 * (t[surface_idx:] - surface_wave_time))
seismic_signal[surface_idx:surface_idx + len(surface_wave)] += surface_wave
添加环境噪声
seismic_signal += 0.3 * np.random.normal(0, 0.1, len(t))
使用Hann窗进行时频分析
window_length = 256
hann_win = hann(window_length)
计算频谱图
f, t_spec, Sxx = signal.spectrogram(seismic_signal, fs, window=hann_win,
nperseg=window_length, noverlap=window_length // 2)
地震信号分析完成:
print("检测到地震波序列: P波 → S波 → 面波")
print("Hann窗提供良好的时频局部化特性")
#检测到地震波序列: P波 → S波 → 面波
#Hann窗提供良好的时频局部化特性
示例7. 电力系统 - 谐波分析
模拟电网电压信号（50Hz基波 + 谐波）
fs = 3200 # 3.2kHz采样率
t = np.linspace(0, 0.2, 640) # 10个周期（200ms）
理想50Hz正弦波 + 谐波污染
fundamental = 230 * np.sqrt(2) * np.sin(2 * np.pi * 50 * t) # 230V RMS
谐波成分：3次(2%), 5次(5%), 7次(1%), 11次(0.5%)
harmonics = (0.02 * 230 * np.sqrt(2) * np.sin(2 * np.pi * 150 * t) + # 3次谐波
0.05 * 230 * np.sqrt(2) * np.sin(2 * np.pi * 250 * t) + # 5次谐波
0.01 * 230 * np.sqrt(2) * np.sin(2 * np.pi * 350 * t) + # 7次谐波
0.005 * 230 * np.sqrt(2) * np.sin(2 * np.pi * 550 * t)) # 11次谐波
voltage_signal = fundamental + harmonics + 5 * np.random.normal(0, 1, len(t))
应用Hann窗进行精确的谐波分析
N = len(voltage_signal)
hann_win = hann(N)
windowed_voltage = voltage_signal * hann_win
FFT分析
spectrum = np.fft.fft(windowed_voltage)
freqs = np.fft.fftfreq(N, 1 / fs)
计算各次谐波的THD（总谐波失真）
fundamental_power = 0
harmonic_power = 0
for i, freq in enumerate(freqs[:N // 2]):
if 45 <= abs(freq) <= 55: # 基波
fundamental_power += np.abs(spectrum[i]) ** 2
elif abs(freq) > 55: # 谐波
harmonic_power += np.abs(spectrum[i]) ** 2
thd = np.sqrt(harmonic_power / fundamental_power) * 100
print("电能质量分析结果:")
print(f" 基波频率: 50 Hz")
print(f" 总谐波失真(THD): {thd:.2f}%")
print(f" Hann窗有效抑制频谱泄漏，提高谐波测量精度")
#电能质量分析结果:
#基波频率: 50 Hz
#总谐波失真(THD): 5.94%
#Hann窗有效抑制频谱泄漏，提高谐波测量精度
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.signal as signal
def hann_window(input_str):
"""
生成Hann窗（汉宁窗），对标MATLAB的hann函数。
参数:
input_str (str): 输入字符串，支持以下格式：
- 单个数字（如"64"）生成对称窗（'symmetric'）。
- 元组（如"(64, 0)"），第二个元素为0生成周期窗（'periodic'），非0生成对称窗。
返回:
list or str: Hann窗列表，或错误信息。
"""
try:
expr = sp.sympify(input_str)
result = None
# 处理元组输入（例如 "(m, flag)"）
if isinstance(expr, tuple) and len(expr) == 2:
m, flag = expr
if m.is_number and flag.is_number:
m_int = int(m)
# flag=0 对应周期窗（sym=False），否则对称窗（sym=True）
sym = bool(int(flag) != 0)
window = signal.windows.hann(m_int, sym)
result = list(window)
else:
return f"输入错误: 参数必须为数字 {input_str}"
# 处理单个数字输入（例如 "m"）
elif expr.is_number:
m_int = int(expr)
window = signal.windows.hann(m_int, sym=True) # 默认对称窗
result = list(window)
else:
return f"输入格式错误: {input_str}"
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 生成对称窗（64点，MATLAB默认）
print("对称窗（64点）:", hann_window("64")[:5])
# [0.0, 0.002484612317299295, 0.009913756075728086, 0.022213597106929606, 0.03926189406479613]
# 生成周期窗（64点，第二个参数为0）
print("周期窗（64点）:", hann_window("(64, 0)")[:5])
# [0.0, 0.002407636663901591, 0.009607359798384785, 0.021529832133895588, 0.03806023374435663]
其特征值位于复平面中的垂直线上的矩阵
H = hanowa(n,alpha) 返回 2×2 分块矩阵, 其中包含四个 n/2×n/2 分块，形式如下.
[alpha*eye(n/2), -diag(1:n/2)
diag(1:n/2), alpha*eye(n/2)]
参量 n 必须为偶数整数. alpha 的默认值为 -1. A 具有复数特征值，形式为 alpha ± k*i (1 <= k <= n/2).
n, alpha — 输入,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def hanowa_matrix(input_str):
"""
构造Hanowa矩阵，对标MATLAB的hanowa函数
参数:
input_str (str): 输入字符串，格式为：
- 单个偶数n: 生成n×n矩阵，默认d=-1
- 元组(n, d): n为偶数，d为标量参数
返回:
Sympy矩阵或错误信息
Hanowa矩阵结构：
[[d*I, -D],
[ D, d*I]]
其中D = diag(1,2,...,m), m = n/2, I为m×m单位矩阵
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def build_hanowa(n, d=-1):
"""核心构造函数"""
# 验证n是否为偶数
if n % 2 != 0 or n <= 0:
raise ValueError("n必须是正偶数")
m = n // 2
I = sp.eye(m) # 单位矩阵
D = sp.diag(*range(1, m + 1)) # 对角矩阵1~m
# 分块矩阵构造
top_row = sp.Matrix.hstack(d * I, -D)
bottom_row = sp.Matrix.hstack(D, d * I)
return sp.Matrix.vstack(top_row, bottom_row)
# 输入解析逻辑
if isinstance(expr, tuple) and len(expr) == 2:
# 处理(n, d)元组输入
if expr[0].is_Integer and expr[1].is_number:
n = int(expr[0])
d = float(expr[1]) if expr[1].is_Float else int(expr[1])
result = build_hanowa(n, d)
else:
error = True
elif expr.is_Integer:
# 处理单个整数输入
n = int(expr)
result = build_hanowa(n)
else:
error = True
if error:
return f"输入错误: {input_str}\n格式应为'(n, d)'或'n'，其中n为正偶数"
return result
except ValueError as ve:
return f"参数错误: {ve}"
except Exception as e:
return f"运行时错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：有效输入 (4, -1)
input1 = '(4, -1)'
mat1 = hanowa_matrix(input1)
print("示例1结果(n=4, d=-1):\n", mat1)
# Matrix([[-1, 0, -1, 0],
# [0, -1, 0, -2],
# [1, 0, -1, 0],
# [0, 2, 0, -1]])
# 示例2：默认d值测试
input2 = '6'
mat2 = hanowa_matrix(input2)
print("示例2结果(n=6, d=-1):\n", mat2[:, :2]) # 显示前两列
# Matrix([[-1, 0],
# [0, -1],
# [0, 0],
# [1, 0],
# [0, 2],
# [0, 0]])
# 示例3：自定义d值
input3 = '(2, 5)'
mat3 = hanowa_matrix(input3)
print("示例3结果(n=2, d=5):\n", mat3)
# Matrix([[5, -1],
# [1, 5]])
汉克尔矩阵
H = hankel(c) 返回正方形汉克尔矩阵,其中c定义矩阵的第一列,主反对角线以下的元素为零.
H = hankel(c,r) 返回汉克尔矩阵,第一列为c,最后一行为r.如果c的最后一个元素不同于r的第一个元素,则hankel会发出警告,并对反对角线使用c的最后一个元素
c — 汉克尔矩阵的第一列,标量,向量
r — 汉克尔矩阵的最后一行,标量,向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def hankel_matrix(input_str):
"""
构造Hankel矩阵，对标MATLAB的hankel函数。
参数:
input_str (str): 输入字符串，格式为：
- 两个列表的元组：'([c1, c2, ...], [r1, r2, ...])'，生成由列向量c和行向量r定义的Hankel矩阵。
- 单个列表：'[a1, a2, ..., an]'，生成n×n的Hankel矩阵，右上三角填充0。
返回:
Sympy矩阵或错误信息。
"""
try:
expr = sp.sympify(input_str)
H = None
error = False
# 处理两个参数的情况：列向量c和行向量r
if isinstance(expr, tuple) and len(expr) == 2 and \
isinstance(expr[0], list) and isinstance(expr[1], list):
c, r = expr[0], expr[1]
rows, cols = len(c), len(r)
H = sp.Matrix.zeros(rows, cols)
for i in range(rows):
for j in range(cols):
if i + j < rows:
H[i, j] = c[i + j] # 填充c的元素
else:
idx = i + j - rows # 计算r的索引
H[i, j] = r[idx] if idx < len(r) else 0 # 超出部分填充0
# 处理单个参数的情况，生成方阵
elif isinstance(expr, list):
c = expr
n = len(c)
H = sp.Matrix.zeros(n, n)
for i in range(n):
for j in range(n):
H[i, j] = c[i + j] if i + j < n else 0 # 右上三角填充0
else:
error = True
return H if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：两个参数生成Hankel矩阵
input1 = '([1, 2, 3], [3, 4, 5])'
result1 = hankel_matrix(input1)
print("示例1结果（两参数）:\n", result1)
# Matrix([[1, 2, 3],
# [2, 3, 3],
# [3, 3, 4]])
# 示例2：单个参数生成方阵
input2 = '[1, 2, 3]'
result2 = hankel_matrix(input2)
print("示例2结果（单参数）:\n", result2)
# Matrix([[1, 2, 3],
# [2, 3, 0],
# [3, 0, 0]])
# 示例3：r长度不足时填充0
input3 = '([1, 2, 3], [4, 5])'
result3 = hankel_matrix(input3)
print("示例3结果（r长度不足）:\n", result3)
# Matrix([[1, 2],
# [2, 3],
# [3, 4]])
第一类汉克尔函数
HankelHone(n,z)在复数z平面上有一个从@inf到0的分支切割不连续点.
n -- 输入, 数值,符号变量,向量,矩阵,多维数组
z -- 输入, 数值,符号变量,向量,矩阵,多维数组
示例1. 电磁波传播 - 圆柱形天线的辐射场
天线参数
frequency = 2.4e9 # 2.4GHz (WiFi频段)
wavelength = 3e8 / frequency
k = 2 * np.pi / wavelength # 波数
观察点位置（极坐标）
rho = np.linspace(0.1, 5, 100) * wavelength # 径向距离
phi = np.pi / 4 # 方位角
计算不同阶数的汉克尔函数（表示不同模态）
orders = [0, 1, 2] # 0阶、1阶、2阶模态
H_fields = []
for n in orders:
# 第一类汉克尔函数 H1_n(kρ) 表示外向波
Hn = [HankelHone(n, k * r) for r in rho]
H_fields.append(Hn)
计算电场强度（简化模型）
E_field = np.abs(H_fields[0]) # 主要考虑0阶模态
print(f"圆柱天线辐射场分析 (频率: {frequency / 1e9} GHz)")
print(f"波数 k = {k:.3f} rad/m")
print(f"电场随距离衰减特性计算完成")
#圆柱天线辐射场分析 (频率: 2.4 GHz)
#波数 k = 50.265 rad/m
#电场随距离衰减特性计算完成
示例2. 声学散射 - 圆柱体声波散射
声波参数
sound_speed = 343 # 空气中声速 m/s
frequency = 1000 # 1kHz
k = 2 * np.pi * frequency / sound_speed # 声波波数
圆柱体半径
a = 0.5 # 0.5米半径
入射平面波角度
theta_incident = 0 # 沿x轴方向入射
计算散射场系数（使用汉克尔函数）
max_order = 10 # 最大展开阶数
scattering_coeffs = []
for n in range(max_order + 1):
# 第一类汉克尔函数 H1_n(ka)
Hn_ka = HankelHone(n, k * a)
# 散射系数（简化模型）
if n == 0:
coeff = - Hn_ka / (Hn_ka + 1j) # 0阶散射系数
else:
coeff = - (2 * Hn_ka) / (Hn_ka + 1j) # 高阶散射系数
scattering_coeffs.append(coeff)
计算远场散射模式
theta = np.linspace(0, 2 * np.pi, 360)
scattering_pattern = np.zeros_like(theta, dtype=complex)
for i, angle in enumerate(theta):
for n in range(max_order + 1):
scattering_pattern[i] += scattering_coeffs[n] * np.cos(n * (angle - theta_incident))
print("圆柱体声波散射分析完成")
print(f"圆柱半径: {a} m, 频率: {frequency} Hz")
print(f"使用{max_order}阶汉克尔函数展开")
#圆柱体声波散射分析完成
#圆柱半径: 0.5 m, 频率: 1000 Hz
#使用10阶汉克尔函数展开
示例3. 量子力学 - 球对称势散射
粒子参数（中子散射）
hbar = 1.0545718e-34 # 约化普朗克常数
mass = 1.6749275e-27 # 中子质量 kg
energy = 1.0e-19 # 粒子能量 J
计算波数
k = np.sqrt(2 * mass * energy) / hbar
势垒半径
R = 1e-15 # 1飞米（原子核尺度）
计算分波展开的相移（使用汉克尔函数）
l_values = [0, 1, 2] # 角动量量子数
phase_shifts = []
scattering_amplitudes = []
for l in l_values:
# 球汉克尔函数 h1_l(kr) = sqrt(pi/2kr) * H1_{l+1/2}(kr)
kr = k * R
H_l_half = HankelHone(l + 0.5, kr)
spherical_hankel = np.sqrt(np.pi / (2 * kr)) * H_l_half
# 简化模型计算相移（假设硬球势）
j_l = np.sin(kr - l * np.pi / 2) / (kr) # 球贝塞尔函数近似
n_l = -np.cos(kr - l * np.pi / 2) / (kr) # 球诺依曼函数近似
tan_delta = j_l / n_l
delta_l = np.arctan(tan_delta)
phase_shifts.append(delta_l)
# 散射振幅
f_l = (2 * l + 1) * np.exp(1j * delta_l) * np.sin(delta_l) / k
scattering_amplitudes.append(f_l)
print("量子力学散射分析:")
print(f"粒子能量: {energy:.2e} J, 波数: {k:.2e} 1/m")
print(f"角动量分波: l={l_values}")
print(f"相移: {[f'{delta:.4f}' for delta in phase_shifts]}")
#量子力学散射分析:
#粒子能量: 1.00e-19 J, 波数: 1.74e+11 1/m
#角动量分波: l=[0, 1, 2]
#相移: ['-0.0002', '1.5706', '-0.0002']
示例4. 地震波传播 - 层状介质中的波场
地层参数
layer_thickness = [100, 200, np.inf] # 各层厚度 (m)，最后一层半无限
wave_speeds = [800, 1200, 1800] # 各层波速 (m/s)
densities = [1800, 2000, 2200] # 各层密度 (kg/m³)
frequency = 10 # 地震波频率 10Hz
omega = 2 * np.pi * frequency
计算各层的波数
wave_numbers = [omega / v for v in wave_speeds]
使用汉克尔函数表示柱面波传播
r_values = np.linspace(10, 1000, 50) # 距离震源的距离
wave_fields = []
for i, (k, rho) in enumerate(zip(wave_numbers, densities)):
layer_field = []
for r in r_values:
# 第一类汉克尔函数 H1_0(kr) 表示外向柱面波
H0 = HankelHone(0, k * r)
# 波场幅度（考虑几何扩散和介质特性）
amplitude = np.abs(H0) / np.sqrt(r) * rho
layer_field.append(amplitude)
wave_fields.append(layer_field)
print("层状介质地震波传播建模:")
print(f"频率: {frequency} Hz, 层数: {len(layer_thickness)}")
print("使用0阶汉克尔函数表示柱面波传播")
#层状介质地震波传播建模:
#频率: 10 Hz, 层数: 3
#使用0阶汉克尔函数表示柱面波传播
示例5. 微波工程 - 圆柱波导模式分析
波导参数
radius = 0.05 # 5cm半径
frequency = 10e9 # 10GHz
omega = 2 * np.pi * frequency
c = 3e8 # 光速
计算截止波数（TM模式）
modes_tm = [] # TM_{mn} 模式
modes_te = [] # TE_{mn} 模式
TM模式的特征方程根（贝塞尔函数的根）
bessel_roots_tm = {
(0, 1): 2.4048, (1, 1): 3.8317, (2, 1): 5.1356,
(0, 2): 5.5201, (1, 2): 7.0156, (2, 2): 8.4172
}
TE模式的特征方程根（贝塞尔函数导数的根）
bessel_roots_te = {
(0, 1): 3.8317, (1, 1): 1.8412, (2, 1): 3.0542,
(0, 2): 7.0156, (1, 2): 5.3314, (2, 2): 6.7061
}
计算TM模式
for (m, n), root in bessel_roots_tm.items():
kc = root / radius # 截止波数
cutoff_freq = c * kc / (2 * np.pi)
if frequency > cutoff_freq:
# 传播常数
beta = np.sqrt((omega / c) ** 2 - kc ** 2)
# 场分布使用汉克尔函数
r_test = radius / 2
H_field = HankelHone(m, kc * r_test)
modes_tm.append({
'mode': f'TM_{m}{n}',
'cutoff_freq': cutoff_freq,
'propagation_constant': beta,
'field_sample': H_field
})
print("圆柱波导模式分析:")
print(f"波导半径: {radius}m, 工作频率: {frequency / 1e9}GHz")
print(f"支持的TM模式: {[mode['mode'] for mode in modes_tm]}")
#圆柱波导模式分析:
#波导半径: 0.05m, 工作频率: 10.0GHz
#支持的TM模式: ['TM_01', 'TM_11', 'TM_21', 'TM_02', 'TM_12', 'TM_22']
示例6. 水下声学 - 海洋波导中的声传播
海洋环境参数
sound_speed_water = 1500 # 水中声速 m/s
sound_speed_bottom = 1600 # 海底声速 m/s
depth = 100 # 海水深度 m
frequency = 500 # 500Hz声源
k_water = 2 * np.pi * frequency / sound_speed_water
距离范围
ranges = np.linspace(100, 5000, 100) # 100m到5km
使用汉克尔函数表示柱面波（简正波理论）
propagation_loss = []
for r in ranges:
# 简正波求和（简化模型，只考虑前几阶）
total_field = 0
num_modes = 5
for n in range(1, num_modes + 1):
# 第n阶简正波的横向波数
k_n = k_water * np.sin(n * np.pi / depth)
# 使用汉克尔函数表示径向传播
H0 = HankelHone(0, k_n * r)
# 简正波幅度（深度函数简化）
mode_amplitude = np.sin(n * np.pi * 0.5) # 在50m深度
total_field += mode_amplitude * H0
propagation_loss.append(20 * np.log10(np.abs(total_field) + 1e-10))
print("海洋声波传播建模:")
print(f"频率: {frequency}Hz, 水深: {depth}m")
print(f"使用{num_modes}阶简正波和0阶汉克尔函数")
#海洋声波传播建模:
#频率: 500Hz, 水深: 100m
#使用5阶简正波和0阶汉克尔函数
示例7. 光学 - 圆柱形微粒的光散射
光学参数
wavelength = 0.6328e-6 # He-Ne激光波长 632.8nm
k0 = 2 * np.pi / wavelength # 真空波数
圆柱微粒参数
radius = 1e-6 # 1微米半径
n_particle = 1.59 # 颗粒折射率（玻璃）
n_medium = 1.33 # 介质折射率（水）
相对折射率
m = n_particle / n_medium
k = k0 * n_medium # 介质中的波数
x = k * radius # 尺寸参数
计算散射系数（米氏理论）
max_terms = 20
scattering_coeffs_an = [] # an系数
scattering_coeffs_bn = [] # bn系数
for n in range(1, max_terms + 1):
# 计算汉克尔函数及其导数
Hn_x = HankelHone(n, x)
Hn_mx = HankelHone(n, m * x)
# 数值导数（中心差分）
dx = 1e-8
Hn_x_plus = HankelHone(n, x + dx)
Hn_x_minus = HankelHone(n, x - dx)
Hn_prime_x = (Hn_x_plus - Hn_x_minus) / (2 * dx)
Hn_mx_plus = HankelHone(n, m * x + dx)
Hn_mx_minus = HankelHone(n, m * x - dx)
Hn_prime_mx = (Hn_mx_plus - Hn_mx_minus) / (2 * dx)
# 散射系数（简化计算）
an = (m * Hn_prime_mx * Hn_x - Hn_mx * Hn_prime_x) / \
(m * Hn_prime_mx * Hn_x - Hn_mx * Hn_prime_x + 1j * (Hn_mx * Hn_x - m * Hn_prime_mx * Hn_prime_x))
bn = (Hn_prime_mx * Hn_x - m * Hn_mx * Hn_prime_x) / \
(Hn_prime_mx * Hn_x - m * Hn_mx * Hn_prime_x + 1j * (m * Hn_mx * Hn_x - Hn_prime_mx * Hn_prime_x))
scattering_coeffs_an.append(an)
scattering_coeffs_bn.append(bn)
# 计算散射效率
Q_sca = 0
for n in range(1, max_terms + 1):
Q_sca += (2 * n + 1) * (np.abs(scattering_coeffs_an[n - 1]) ** 2 +
np.abs(scattering_coeffs_bn[n - 1]) ** 2)
Q_sca *= 2 / (x ** 2)
print("圆柱微粒米氏散射分析:")
print(f"波长: {wavelength * 1e9:.1f}nm, 颗粒半径: {radius * 1e6:.1f}μm")
print(f"散射效率 Q_sca = {Q_sca:.4f}")
print(f"使用{max_terms}项汉克尔函数展开")
#圆柱微粒米氏散射分析:
#波长: 632.8nm, 颗粒半径: 1.0μm
#散射效率 Q_sca = 1.9334
#使用20项汉克尔函数展开
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.special as sci
def hankel1_bessel_representation(input_str):
"""
计算第一类汉克尔函数H1_v(x)，对标Wolfram的HankelH1[v,x]。
参数:
input_str (str): 输入字符串，格式为 "(v, x)"，其中：
- v 可以是标量或矩阵
- x 可以是标量或矩阵
- 若v和x均为矩阵，则形状必须相同
返回:
SymPy矩阵/数值 或 str: 计算结果或错误信息
"""
try:
expr = sp.sympify(input_str)
# 检查输入是否为长度为2的SymPy元组（例如 "(v, x)"）
if not (isinstance(expr, tuple) and len(expr) == 2):
return f"输入错误: 需要两个参数组成的元组，例如 '(v, x)'"
if all(e.is_number for e in expr):
v = float(expr[0])
x = complex(expr[1])
result = sci.hankel1(v, x)
elif any(e.free_symbols for e in expr):
result = sp.hankel1(*expr)
else:
error = True
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 标量计算
print("H1_2(5) =", hankel1_bessel_representation("2, 5+2j"))
# (0.02047643629410032+0.050063353169436886j)
# 示例2: v为标量，x为矩阵
x_matrix = "1"
print("\nH1_1(x_matrix):")
print(hankel1_bessel_representation(f"(1, {x_matrix})"))
# (0.44005058574493355-0.7812128213002889j)
v_matrix = "2"
print("\nH1_v_matrix(x_matrix):")
print(hankel1_bessel_representation(f"({v_matrix}, {x_matrix})"))
# (0.11490348493190057-1.6506826068162548j)
第二类汉克尔函数
HankelHtwo(n,z)在复数z平面上有一个从@inf到0的分支切割不连续点.
n -- 输入, 数值,符号变量,向量,矩阵,多维数组
z -- 输入, 数值,符号变量,向量,矩阵,多维数组
示例1. 电磁波吸收 - 完美匹配层(PML)理论
PML参数
frequency = 10e9 # 10GHz
wavelength = 3e8 / frequency
k0 = 2 * np.pi / wavelength # 自由空间波数
PML层中的复坐标拉伸
sigma_max = 10 # 最大电导率
pml_thickness = 0.1 # PML厚度 (m)
在PML中的位置（复数坐标）
z_values = np.linspace(0, pml_thickness, 50)
complex_z = z_values + 1j * sigma_max * (z_values / pml_thickness) ** 2 / (2 * k0)
使用第二类汉克尔函数表示被吸收的波
wave_amplitude = []
for z in complex_z:
# H2_0(k0z) 表示内向衰减波
H2 = HankelHtwo(0, k0 * z)
wave_amplitude.append(np.abs(H2))
print("完美匹配层(PML)波吸收分析:")
print(f"频率: {frequency / 1e9}GHz, PML厚度: {pml_thickness}m")
print("使用第二类汉克尔函数模拟波在PML中的衰减")
#完美匹配层(PML)波吸收分析:
#频率: 10.0GHz, PML厚度: 0.1m
#使用第二类汉克尔函数模拟波在PML中的衰减
示例2. 声学成像 - 时间反转镜技术
声学参数
sound_speed = 1500 # 水中声速 m/s
frequency = 100e3 # 100kHz超声
k = 2 * np.pi * frequency / sound_speed
传感器阵列位置（圆形阵列）
num_sensors = 32
array_radius = 0.1 # 10cm半径
sensor_angles = np.linspace(0, 2 * np.pi, num_sensors, endpoint=False)
目标点位置
target_r = 0.05 # 5cm距离
target_theta = np.pi / 3 # 60度方向
时间反转过程：接收和重新发射波场
time_reversed_field = np.zeros(len(sensor_angles), dtype=complex)
for i, sensor_angle in enumerate(sensor_angles):
# 计算从目标到传感器的距离
distance = np.sqrt(array_radius ** 2 + target_r ** 2 -
2 * array_radius * target_r * np.cos(sensor_angle - target_theta))
# 接收过程使用H2（内向波）
H2_receive = HankelHtwo(f"(0, {k * distance})")
# 时间反转后重新发射（相位共轭）
time_reversed_field[i] = np.conj(H2_receive)
重构图像场
image_plane = np.zeros((50, 50), dtype=complex)
x_vals = np.linspace(-0.1, 0.1, 50)
y_vals = np.linspace(-0.1, 0.1, 50)
for i, x in enumerate(x_vals):
for j, y in enumerate(y_vals):
r = np.sqrt(x ** 2 + y ** 2)
if r > 0:
theta = np.arctan2(y, x)
# 合成焦点场使用H1（外向波）
for sensor_idx, sensor_angle in enumerate(sensor_angles):
sensor_x = array_radius * np.cos(sensor_angle)
sensor_y = array_radius * np.sin(sensor_angle)
dist_to_sensor = np.sqrt((x - sensor_x) ** 2 + (y - sensor_y) ** 2)
H1_transmit = HankelHtwo(0, k * dist_to_sensor)
image_plane[i, j] += time_reversed_field[sensor_idx] * H1_transmit
print("时间反转声学成像模拟:")
print(f"阵列: {num_sensors}个传感器, 半径: {array_radius}m")
print(f"目标位置: r={target_r}m, θ={np.degrees(target_theta):.1f}°")
#时间反转声学成像模拟:
#阵列: 32个传感器, 半径: 0.1m
#目标位置: r=0.05m, θ=60.0°
示例3. 量子力学 - 束缚态问题的渐近行为
量子系统参数（氢原子类似势）
hbar = 1.0545718e-34
mass = 9.1093837e-31 # 电子质量
energy = -13.6 * 1.602e-19 # 束缚态能量 -13.6eV
有效波数（虚数，表示指数衰减）
k_imag = np.sqrt(-2 * mass * energy) / hbar
径向距离
r_values = np.linspace(1e-11, 5e-10, 100) # 0.1Å 到 5Å
计算s波（l=0）束缚态波函数的渐近形式
对于束缚态，r→∞时波函数~ H2_0(iκr) ~ e^{-κr}
bound_state_wavefunction = []
for r in r_values:
# 使用第二类汉克尔函数表示指数衰减
# H2_0(iκr) 在 r→∞ 时表现为 e^{-κr}
H2 = HankelHtwo(0, 1j * k_imag * r)
bound_state_wavefunction.append(np.real(H2))
验证指数衰减特性
theoretical_decay = np.exp(-k_imag * r_values)
print("量子束缚态波函数分析:")
print(f"束缚能: {energy / 1.602e-19:.2f} eV")
print(f"衰减常数: κ = {k_imag:.2e} m⁻¹")
print("使用第二类汉克尔函数描述束缚态渐近行为")
#量子束缚态波函数分析:
#束缚能: -13.60 eV
#衰减常数: κ = 1.89e+10 m⁻¹
#使用第二类汉克尔函数描述束缚态渐近行为
示例4. 地震学 - 地下反射波分析
地震参数
dominant_freq = 30 # 30Hz主频
sound_speed = 2000 # 地层速度 m/s
k = 2 * np.pi * dominant_freq / sound_speed
地下反射界面
reflector_depth = 1000 # 1km深度
geophone_spacing = 50 # 50m道间距
num_geophones = 41 # 21个检波器
地面检波器位置
geophone_positions = np.linspace(-1000, 1000, num_geophones)
计算下行波（震源到反射点）和上行波（反射点到接收点）
使用H2表示下行波（从震源向下传播）
down_going_waves = []
up_going_waves = []
for x in geophone_positions:
# 震源在(0,0)，反射点在(x/2, reflector_depth)
down_distance = np.sqrt((x / 2) ** 2 + reflector_depth ** 2)
up_distance = np.sqrt((x / 2) ** 2 + reflector_depth ** 2)
total_distance = down_distance + up_distance
# 下行波（使用H2表示内向传播特性）
H2_down = HankelHtwo(0, k * down_distance)
# 上行波（使用H1表示外向传播特性）
H1_up = HankelHtwo(f"(0, {k * up_distance})")
down_going_waves.append(H2_down)
up_going_waves.append(H1_up)
总波场（反射响应）
reflection_response = [down * up for down, up in zip(down_going_waves, up_going_waves)]
print("地下反射波分析:")
print(f"反射深度: {reflector_depth}m, 道数: {num_geophones}")
print("使用H2表示下行波，H1表示上行波")
#地下反射波分析:
#反射深度: 1000m, 道数: 41
#使用H2表示下行波，H1表示上行波
示例5. 光学谐振腔 - 内向传播模式
谐振腔参数
wavelength = 1.55e-6 # 通信波长 1550nm
k0 = 2 * np.pi / wavelength
resonator_radius = 10e-6 # 10微米半径
材料参数
n_core = 3.45 # 硅波导折射率
n_clad = 1.44 # SiO2包层折射率
谐振条件：波程 = m * 波长
m_values = np.arange(50, 100) # 模式阶数
resonant_frequencies = []
quality_factors = []
for m in m_values:
# 谐振条件：2πR * n_eff = m * λ
circumference = 2 * np.pi * resonator_radius
n_eff_required = m * wavelength / circumference
# 使用第二类汉克尔函数分析内向耦合
# 对于环形谐振腔，内向波有助于形成谐振
k_eff = k0 * n_eff_required
# 计算谐振场的径向依赖
r_test = resonator_radius
H2_field = HankelHtwo(m, k_eff * r_test)
# 估计品质因子（简化模型）
# Q ≈ 2π × 存储能量 / 每周期损耗能量
field_magnitude = np.abs(H2_field)
Q_approx = 1e4 * field_magnitude # 简化估计
resonant_frequencies.append(n_eff_required)
quality_factors.append(Q_approx)
print("光学环形谐振腔模式分析:")
print(f"谐振腔半径: {resonator_radius * 1e6:.1f}μm")
print(f"分析模式阶数: {m_values[0]} 到 {m_values[-1]}")
print("使用第二类汉克尔函数分析内向传播谐振模式")
#光学环形谐振腔模式分析:
#谐振腔半径: 10.0μm
#分析模式阶数: 50 到 99
#使用第二类汉克尔函数分析内向传播谐振模式
示例6. 水下通信 - 收敛区传播
海洋环境参数（深海声道）
sound_speed_profile = [1500, 1480, 1520] # 声速剖面 m/s
layer_depths = [0, 500, 1000] # 层深度 m
frequency = 500 # 500Hz
计算会聚区位置（简化的射线理论）
source_depth = 100 # 声源深度 m
max_range = 50e3 # 50km最大距离
会聚区大约出现在~50-60km间隔
convergence_ranges = [50e3, 110e3] # 第一、第二会聚区
在会聚区使用第二类汉克尔函数表示聚焦波
range_samples = np.linspace(45e3, 55e3, 100) # 第一会聚区附近
convergence_gain = []
for r in range_samples:
# 有效波数（考虑平均声速）
c_avg = np.mean(sound_speed_profile)
k = 2 * np.pi * frequency / c_avg
# 在会聚区，波前曲率导致能量集中
# 使用高阶汉克尔函数表示聚焦效应
H2_focus = HankelHtwo(10, k * r) # 高阶表示聚焦
# 会聚增益（相对于球面扩展）
spherical_spread = 1 / r
convergence_gain.append(np.abs(H2_focus) / spherical_spread)
print("水下声学会聚区分析:")
print(f"频率: {frequency}Hz, 声源深度: {source_depth}m")
print(f"第一会聚区范围: {convergence_ranges[0] / 1000:.0f}km")
print("使用高阶第二类汉克尔函数模拟会聚效应")
#水下声学会聚区分析:
#频率: 500Hz, 声源深度: 100m
#第一会聚区范围: 50km
#使用高阶第二类汉克尔函数模拟会聚效应
示例7. 电磁兼容性 - 屏蔽效能分析
屏蔽体参数（圆柱形屏蔽壳）
frequency = 1e9 # 1GHz
k0 = 2 * np.pi * frequency / 3e8
shield_radius = 0.1 # 10cm半径
shield_thickness = 0.001 # 1mm厚度
材料参数（铜）
conductivity = 5.8e7 # 铜电导率 S/m
mu0 = 4 * np.pi * 1e-7
epsilon0 = 8.854e-12
金属中的复波数
k_metal = np.sqrt(1j * 2 * np.pi * frequency * mu0 * conductivity)
计算屏蔽效能（内外场比）
observation_angles = np.linspace(0, 2 * np.pi, 36)
shielding_effectiveness = []
for phi in observation_angles:
# 外部入射场（平面波近似）
# 内部场使用第二类汉克尔函数表示衰减波
r_inside = shield_radius - shield_thickness / 2
# 对于良好导体，内部场急剧衰减
H2_inside = HankelHtwo(f"(0, {k_metal * r_inside})")
# 外部参考场（无屏蔽）
H2_outside = HankelHtwo(0, k0 * shield_radius)
# 屏蔽效能 = 20log10(外部场/内部场)
SE = 20 * np.log10(np.abs(H2_outside) / (np.abs(H2_inside) + 1e-10))
shielding_effectiveness.append(SE)
print("电磁屏蔽效能分析:")
print(f"频率: {frequency / 1e9}GHz, 屏蔽体半径: {shield_radius}m")
print(f"材料: 铜 (σ={conductivity:.1e} S/m)")
print("使用第二类汉克尔函数分析屏蔽体内的场衰减")
#电磁屏蔽效能分析:
#频率: 1.0GHz, 屏蔽体半径: 0.1m
#材料: 铜 (σ=5.8e+07 S/m)
#使用第二类汉克尔函数分析屏蔽体内的场衰减
示例8. 地球物理学 - 地幔柱上升流模拟
地幔柱参数
plume_radius = 100e3 # 100km半径
thermal_diffusivity = 1e-6 # 热扩散系数 m²/s
upwelling_velocity = 0.1 # 上升速度 m/year → m/s
upwelling_velocity_ms = upwelling_velocity / (365 * 24 * 3600)
特征时间和长度尺度
characteristic_time = plume_radius ** 2 / thermal_diffusivity
peclet_number = upwelling_velocity_ms * plume_radius / thermal_diffusivity
柱坐标系中的温度场（稳态）
r_values = np.linspace(0, 2 * plume_radius, 50)
z_values = np.linspace(0, 500e3, 50) # 500km深度范围
temperature_field = np.zeros((len(r_values), len(z_values)))
for i, r in enumerate(r_values):
for j, z in enumerate(z_values):
if r > 0:
# 使用第二类汉克尔函数表示热柱的径向温度分布
# 对于上升流，温度场满足修正的贝塞尔方程
k_thermal = np.sqrt(peclet_number / plume_radius ** 2)
# 简化模型：温度分布 ~ H2_0(αr) * exp(-βz)
H2_radial = HankelHtwo(0, k_thermal * r)
vertical_decay = np.exp(-z / (2 * plume_radius))
temperature_field[i, j] = np.real(H2_radial) * vertical_decay
print("地幔柱热上升流模拟:")
print(f"地幔柱半径: {plume_radius / 1000:.0f}km")
print(f"佩克莱数: {peclet_number:.2e}")
print("使用第二类汉克尔函数描述柱对称热场分布")
#地幔柱热上升流模拟:
#地幔柱半径: 100km
#佩克莱数: 3.17e+02
#使用第二类汉克尔函数描述柱对称热场分布
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.special as sci
def hankel2_bessel_representation(input_str):
"""
计算第一类汉克尔函数H1_v(x)，对标Wolfram的HankelH1[v,x]。
参数:
input_str (str): 输入字符串，格式为 "(v, x)"，其中：
- v 可以是标量或矩阵
- x 可以是标量或矩阵
- 若v和x均为矩阵，则形状必须相同
返回:
SymPy矩阵/数值 或 str: 计算结果或错误信息
"""
try:
expr = sp.sympify(input_str)
# 检查输入是否为长度为2的SymPy元组（例如 "(v, x)"）
if not (isinstance(expr, tuple) and len(expr) == 2):
return f"输入错误: 需要两个参数组成的元组，例如 '(v, x)'"
if all(e.is_number for e in expr):
v = float(expr[0])
x = complex(expr[1])
result = sci.hankel2(v, x)
elif any(e.free_symbols for e in expr):
result = sp.hankel2(*expr)
else:
error = True
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 标量计算
print("H1_2(5) =", hankel2_bessel_representation("2, 5+2j"))
# (-0.0858294011199319-2.289177410281067j)
# 示例2: v为标量，x为矩阵
x_matrix = "1"
print("\nH1_1(x_matrix):")
print(hankel2_bessel_representation(f"(1, {x_matrix})"))
# (0.44005058574493355+0.7812128213002889j)
v_matrix = "2"
print("\nH1_v_matrix(x_matrix):")
print(hankel2_bessel_representation(f"({v_matrix}, {x_matrix})"))
# (0.11490348493190057+1.6506826068162548j)
调和数
Harmonic（x）返回x的调和函数。对于x的整数值，Harmonic（x）生成调和数。
x —— 输入,数字|向量|矩阵|多维数组|符号变量|符号表达式
示例1. 算法分析 - 快速排序的平均比较次数
快速排序的平均比较次数 = 2n * H(n) - 3n
n_values = [10, 100, 1000]
快速排序算法分析:
print("平均比较次数 = 2n * H(n) - 3n")
#平均比较次数 = 2n * H(n) - 3n
for n in n_values:
H_n = Harmonic(n)
avg_comparisons = 2 * n * float(H_n) - 3 * n
print(f"n = {n:6d}: H(n) = {float(H_n):.4f}, 平均比较次数 = {avg_comparisons:,.0f}")
#n = 10: H(n) = 2.9290, 平均比较次数 = 29
#n = 100: H(n) = 5.1874, 平均比较次数 = 737
#n = 1000: H(n) = 7.4855, 平均比较次数 = 11,971
渐近分析
n_large = 1000
H_n_large = Harmonic(n_large)
asymptotic_complexity = 2 * n_large * float(H_n_large) - 3 * n_large
print(f"\n渐近分析 (n = 1,000,000):")
print(f"理论复杂度: O(n log n)")
print(f"实际计算值: {asymptotic_complexity:,.0f} 次比较")
#渐近分析 (n = 1,000,000):
#理论复杂度: O(n log n)
#实际计算值: 11,971 次比较
示例2. 数论 - 调和级数的发散性研究
和级数部分和 H(n) = 1 + 1/2 + ... + 1/n
milestones = [10, 100, 1000]
print("调和级数发散性研究:")
print("H(n) = 1 + 1/2 + 1/3 + ... + 1/n")
#调和级数发散性研究:
#H(n) = 1 + 1/2 + 1/3 + ... + 1/n
results = []
for n in milestones:
H_n = Harmonic(n)
results.append((n, float(H_n)))
print(f"H({n:6d}) = {float(H_n):.6f}")
#H( 10) = 2.928968
#H( 100) = 5.187378
#H( 1000) = 7.485471
与对数的比较
for n, H_n in results:
ln_n = np.log(n)
difference = H_n - ln_n
print(f"H({n}) - ln({n}) = {difference:.6f}")
#H(10) - ln(10) = 0.626383
#H(100) - ln(100) = 0.582207
#H(1000) - ln(1000) = 0.577716
欧拉-马歇罗尼常数近似
gamma_approx = Harmonic(1000) - np.log(1000)
print(f"\n欧拉-马歇罗尼常数 γ ≈ H(1000000) - ln(1000000) = {gamma_approx:.10f}")
#欧拉-马歇罗尼常数 γ ≈ H(1000000) - ln(1000000) = 0.5777155816
示例3. 概率论 - 优惠券收集问题
收集n种不同优惠券的期望次数 = n * H(n)
coupon_types = [5, 10, 20, 50, 100]
优惠券收集问题分析:
print("收集全部n种优惠券的期望次数 = n × H(n)")
#收集全部n种优惠券的期望次数 = n × H(n)
expected_trials = []
for n in coupon_types:
H_n = Harmonic(n)
expected = n * float(H_n)
expected_trials.append(expected)
print(f"n = {n:3d}: 期望次数 = {expected:7.2f}")
#n = 5: 期望次数 = 11.42
#n = 10: 期望次数 = 29.29
#n = 20: 期望次数 = 71.95
#n = 50: 期望次数 = 224.96
#n = 100: 期望次数 = 518.74
方差分析
for n in [10, 20, 50]:
H_n = Harmonic(n)
H_n_sq = Harmonic(n ** 2)
# 方差近似公式
variance_approx = (np.pi ** 2 * n ** 2) / 6 - n * float(H_n) - n ** 2 * float(H_n_sq)
std_dev = np.sqrt(max(variance_approx, 0))
print(f"n = {n}: 标准差 ≈ {std_dev:.2f}")
#n = 10: 标准差 ≈ 0.00
#n = 20: 标准差 ≈ 0.00
#n = 50: 标准差 ≈ 0.00
示例4. 信息论 - 最优编码长度分析
对于n个符号，最优编码的平均长度接近 H(概率分布)
symbol_counts = [8, 16, 32, 64] # 不同符号数量
信息论编码分析:
print("均匀分布下最优编码的平均长度 ≈ H(n)")
#均匀分布下最优编码的平均长度 ≈ H(n)
for n in symbol_counts:
H_n = Harmonic(n)
# 均匀分布下的最优编码平均长度
optimal_length = float(H_n) / np.log(2) # 转换为比特
# 固定长度编码的比较
fixed_length = np.ceil(np.log2(n))
print(f"符号数 {n:2d}: 最优编码 {optimal_length:5.2f} 比特, 固定编码 {fixed_length:2.0f} 比特")
#符号数 8: 最优编码 3.92 比特, 固定编码 3 比特
#符号数 16: 最优编码 4.88 比特, 固定编码 4 比特
#符号数 32: 最优编码 5.86 比特, 固定编码 5 比特
#符号数 64: 最优编码 6.84 比特, 固定编码 6 比特
非均匀分布示例
p = 0.5 # 几何分布参数
theoretical_entropy = -p * np.log2(p) - (1 - p) * np.log2(1 - p)
使用调和数近似计算
n_terms = 100
geometric_entropy_approx = 0
for k in range(1, n_terms + 1):
prob = p * (1 - p) ** (k - 1)
if prob > 1e-10:
geometric_entropy_approx -= prob * np.log2(prob)
print(f"几何分布(p={p})的理论熵: {theoretical_entropy:.4f} 比特")
print(f"数值计算熵: {geometric_entropy_approx:.4f} 比特")
#几何分布(p=0.5)的理论熵: 1.0000 比特
#数值计算熵: 2.0000 比特
示例5. 物理学 - 理想链的熵计算
对于N个链节的理想链，熵与调和数有关
chain_lengths = [10, 50, 100, 500, 1000]
高分子理想链构象熵分析:
print("熵 S ≈ k_B × [ln(2) + H(N-1)]")
#熵 S ≈ k_B × [ln(2) + H(N-1)]
k_B = 1.380649e-23 # 玻尔兹曼常数
entropies = []
for N in chain_lengths:
if N > 1:
H_N_minus_1 = Harmonic(N - 1)
entropy = k_B * (np.log(2) + float(H_N_minus_1))
entropies.append(entropy)
print(f"N = {N:4d}: 熵 = {entropy:.2e} J/K")
#N = 10: 熵 = 4.86e-23 J/K
#N = 50: 熵 = 7.14e-23 J/K
#N = 100: 熵 = 8.11e-23 J/K
#N = 500: 熵 = 1.03e-22 J/K
#N = 1000: 熵 = 1.13e-22 J/K
与随机游走的关系
N_large = 1000
H_large = Harmonic(N_large - 1)
rms_distance = np.sqrt(N_large) # 均方根距离
print(f"N = {N_large} 的随机游走:")
print(f" 均方根距离: {rms_distance:.1f} 步长")
print(f" 构象熵: {k_B * (np.log(2) + float(H_large)):.2e} J/K")
#N = 1000 的随机游走:
# 均方根距离: 31.6 步长
# 构象熵: 1.13e-22 J/K
示例6. 金融数学 - 债券久期计算
麦考利久期计算涉及调和数的加权平均
bond_years = [5, 10, 20, 30] # 债券期限
coupon_rate = 0.05 # 票面利率5%
yield_rate = 0.04 # 市场收益率4%
债券久期计算:
print("麦考利久期 = Σ [t × PV(CF_t)] / 债券价格")
#麦考利久期 = Σ [t × PV(CF_t)] / 债券价格
for years in bond_years:
# 计算久期（简化模型，使用调和数近似）
H_years = Harmonic(years)
# 麦考利久期近似公式
macaulay_duration = (1 + yield_rate) / yield_rate - \
(1 + yield_rate + years * (coupon_rate - yield_rate)) / \
(coupon_rate * ((1 + yield_rate) ** years - 1) + yield_rate)
# 使用调和数的另一种近似
harmonic_approx = float(H_years) * coupon_rate / yield_rate
print(f"{years:2d}年期债券: 理论久期 = {macaulay_duration:5.2f} 年, " +
f"调和近似 = {harmonic_approx:5.2f} 年")
# 5年期债券: 理论久期 = 4.56 年, 调和近似 = 2.85 年
#10年期债券: 理论久期 = 8.19 年, 调和近似 = 3.66 年
#20年期债券: 理论久期 = 13.54 年, 调和近似 = 4.50 年
#30年期债券: 理论久期 = 17.19 年, 调和近似 = 4.99 年
永续债券的特殊情况
perpetuity_duration = (1 + yield_rate) / yield_rate
print(f"永续债券久期 = (1 + r) / r = {perpetuity_duration:.2f} 年")
#永续债券久期 = (1 + r) / r = 26.00 年
示例7. 计算机科学 - 哈希表冲突分析
将m个球投入n个箱子，空箱子数量的期望 = n × (1 - 1/n)^m ≈ n × exp(-m/n)
使用调和数可以得到更精确的结果
hash_table_sizes = [100, 1000, 10000] # 哈希表大小
load_factors = [0.5, 1.0, 2.0] # 负载因子 α = m/n
哈希表冲突分析:
print("负载因子 α = m/n")
#负载因子 α = m/n
results = []
for n in hash_table_sizes:
print(f"\n哈希表大小 n = {n}:")
for alpha in load_factors:
m = int(alpha * n) # 元素数量
# 空槽位数量的精确期望值
expected_empty = n * (1 - 1 / n) ** m
# 使用调和数的近似
if m <= n:
harmonic_approx = n * float(Harmonic(n)) - \
n * float(Harmonic(n - m))
else:
harmonic_approx = 0
# 冲突次数的期望
expected_collisions = m - n + expected_empty
print(f" α = {alpha:.1f}: 期望冲突次数 = {expected_collisions:6.1f}, " +
f"空槽位 = {expected_empty:6.1f}")
results.append((n, alpha, expected_collisions))
#哈希表大小 n = 100:
# α = 0.5: 期望冲突次数 = 10.5, 空槽位 = 60.5
# α = 1.0: 期望冲突次数 = 36.6, 空槽位 = 36.6
# α = 2.0: 期望冲突次数 = 113.4, 空槽位 = 13.4
#哈希表大小 n = 1000:
# α = 0.5: 期望冲突次数 = 106.4, 空槽位 = 606.4
# α = 1.0: 期望冲突次数 = 367.7, 空槽位 = 367.7
# α = 2.0: 期望冲突次数 = 1135.2, 空槽位 = 135.2
#哈希表大小 n = 10000:
# α = 0.5: 期望冲突次数 = 1065.2, 空槽位 = 6065.2
# α = 1.0: 期望冲突次数 = 3678.6, 空槽位 = 3678.6
# α = 2.0: 期望冲突次数 = 11353.2, 空槽位 = 1353.2
生日悖论分析
birthday_probabilities = []
for n_people in [23, 50, 100]:
prob_no_collision = 1.0
for i in range(1, n_people):
prob_no_collision *= (365 - i) / 365
prob_collision = 1 - prob_no_collision
# 使用调和数近似
H_365 = harmonic_number("365")
birthday_approx = 1 - np.exp(-n_people * (n_people - 1) / (2 * 365))
print(f"{n_people} 人: 理论概率 = {prob_collision:.4f}, " +
f"调和近似 = {birthday_approx:.4f}")
#23 人: 理论概率 = 0.5073, 调和近似 = 0.5000
#50 人: 理论概率 = 0.9704, 调和近似 = 0.9651
#100 人: 理论概率 = 1.0000, 调和近似 = 1.0000
示例8. 统计学 - 极端值理论
在极值理论中，n个独立同分布随机变量的最大值的期望与调和数有关
sample_sizes = [10, 50, 100, 1000]
distribution_types = ['指数分布', '均匀分布', '正态分布']
极端值统计分析:
for dist_type in distribution_types:
print(f"\n{dist_type}:")
for n in sample_sizes:
if dist_type == '指数分布':
# 指数分布 λ=1 的最大值期望 = H(n)
expected_max = Harmonic(n)
elif dist_type == '均匀分布':
# 均匀分布 [0,1] 的最大值期望 = n/(n+1)
expected_max = n / (n + 1)
elif dist_type == '正态分布':
# 标准正态分布的最大值期望近似
# 使用极值理论的近似公式
gamma = 0.5772156649 # 欧拉常数
expected_max = np.sqrt(2 * np.log(n)) - \
(gamma + np.log(np.log(n)) + np.log(4 * np.pi)) / \
(2 * np.sqrt(2 * np.log(n)))
print(f" n = {n:4d}: E[max] = {expected_max:.4f}")
#指数分布:
# n = 10: E[max] = 2.9290
# n = 50: E[max] = 4.4992
# n = 100: E[max] = 5.1874
# n = 1000: E[max] = 7.4855
#均匀分布:
# n = 10: E[max] = 0.9091
# n = 50: E[max] = 0.9804
# n = 100: E[max] = 0.9901
# n = 1000: E[max] = 0.9990
#正态分布:
# n = 10: E[max] = 1.2274
# n = 50: E[max] = 1.9977
# n = 100: E[max] = 2.2712
# n = 1000: E[max] = 3.0388
记录值分析:
print("在n个独立同分布的观测中，记录值的期望数量 = H(n)")
#在n个独立同分布的观测中，记录值的期望数量 = H(n)
for n in [10, 100, 1000]:
expected_records = float(Harmonic(n)
print(f"n = {n}: 期望记录值数量 = {expected_records:.4f}")
#n = 10: 期望记录值数量 = 2.9290
#n = 100: 期望记录值数量 = 5.1874
#n = 1000: 期望记录值数量 = 7.4855
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def harmonic_number(input_str):
"""
计算调和数H(n) = 1 + 1/2 + ... + 1/n，对标MATLAB的harmonic函数。
参数:
input_str (str): 输入字符串，可以是：
- 正整数（如"5"）
- 符号表达式（如"n"）
- 矩阵（如"[[1, 2], [3, 4]]"）
返回:
SymPy数值/矩阵/符号表达式 或 str: 计算结果或错误信息
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if expr.is_number or expr.free_symbols:
# 标量或符号变量
# 只有sympy harmonic支持复数, 负数参数，numpy和scipy均不支持。
result = sp.harmonic(expr).evalf()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 标量输入
print("H(5) =", harmonic_number("5+2j"))
# 2.34499020414819 + 0.347995611554697*I
# 示例2: 符号输入
n = sp.symbols('n')
print("\n符号计算 H(n):", harmonic_number("n"))
# harmonic(n)
# 示例3: 非整数输入
print("\n非整数测试:", harmonic_number("2.5"))
# 1.68037230554678
# 示例4: 负数输入
print("\n负数测试:", harmonic_number("-3"))
# zoo
豪斯多夫距离是一种数学构造，用于度量度量空间子集的两组点的“接近度”。
这种度量可用于为两个轨迹、数据云或任何点集之间的相似性分配标量分数。
此函数将返回两组点之间的豪斯多夫距离。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import shapely.geometry as sh
def hausdorff_distance_shapely(input_str):
"""
使用坐标点计算两个几何对象之间的豪斯多夫距离。
参数:
input_str (str): 一个字符串，代表包含两个几何对象的列表。每个对象是一个坐标点列表（元组或列表）。
第一个对象可以是点（1个点）、线（2个点）或多边形（3个及以上点）。第二个对象的定义类似。
返回:
float: 两个几何对象之间的豪斯多夫距离。
示例:
>>> hausdorff_distance_shapely('[[(0,0)], [(1,1)]]')
1.4142135623730951
>>> hausdorff_distance_shapely('[[(0,0), (1,0)], [(2,0), (3,0), (4,0)]]')
3.0
"""
# 将输入字符串解析为Python列表结构
try:
expr = ast.literal_eval(input_str)
except (SyntaxError, ValueError) as e:
raise ValueError("输入格式无效。输入必须是列表的字符串表示形式。") from e
# 检查输入是否恰好包含两个几何对象
if not isinstance(expr, (list, tuple)) or len(expr) != 2:
raise ValueError("输入必须恰好包含两个几何对象。")
a_points, b_points = expr[0], expr[1]
# 从点列表创建Shapely对象的函数
def create_shapely_object(points):
if not isinstance(points, list) or len(points) < 1:
raise ValueError("无效的几何对象点。必须是一个非空的坐标列表。")
# 将所有点转换为元组（以防它们是列表）
points = [tuple(point) for point in points]
if len(points) == 1:
return sh.Point(points[0])
elif len(points) == 2:
return sh.LineString(points)
else:
return sh.Polygon(points)
# 为两个几何对象创建Shapely对象
try:
obj_a = create_shapely_object(a_points)
obj_b = create_shapely_object(b_points)
except Exception as e:
raise ValueError(f"创建几何对象时出错: {e}") from e
# 计算并返回豪斯多夫距离
return obj_a.hausdorff_distance(obj_b)
# 示范代码
if __name__ == "__main__":
# 定义两个几何对象的坐标，这里是一个点和一条线
input_str = "[[(0, 0)], [(1, 1), (2, 2)]]"
# 调用函数计算豪斯多夫距离
result = hausdorff_distance_shapely(input_str)
print(f"两个几何对象之间的豪斯多夫距离为: {result}")
# 2.8284271247461903
海维赛阶跃函数
heaviside(x)在x<0时返回值0,在x>0时返回1,在x＝0时返回1/2.
x —— 标量,向量，矩阵, 多维数组,符号表达式
示例1. 电路分析 - 开关电路响应
在 t=0 时刻施加5V电压
voltage_response = heaviside(5 * Heaviside(t))
print(f"阶跃电压响应: {voltage_response}")
#阶跃电压响应: Heaviside(5*Heaviside(t))
RC电路的阶跃响应
RC_circuit = heaviside(5 * (1 - exp(-t/0.1)) * Heaviside(t))
print(f"RC电路响应: {RC_circuit}\n")
#RC电路响应: Heaviside((5 - 5*exp(-10.0*t))*Heaviside(t))
示例2. 控制系统 - 单位阶跃响应
控制系统传递函数的阶跃响应
control_system = heaviside(1 - exp(-2*t)*cos(3*t))
print(f"控制系统响应: {control_system}")
#控制系统响应: Heaviside(1 - exp(-2*t)*cos(3*t))
示例3. 信号处理 - 矩形脉冲信号
从 t=1 到 t=3 的矩形脉冲
rectangular_pulse = heaviside(Heaviside(t-1) - Heaviside(t-3))
print(f"矩形脉冲信号: {rectangular_pulse}")
#矩形脉冲信号: Heaviside(-Heaviside(t - 3) + Heaviside(t - 1))
示例4. 机械系统 - 冲击载荷
在 t=2 时刻施加冲击力
impact_force = heaviside(100 * Heaviside(t-2))
print(f"冲击载荷: {impact_force}")
#冲击载荷: Heaviside(100*Heaviside(t - 2))
示例5. 经济学 - 政策冲击
政策在特定时间点实施的影响
policy_effect = heaviside(0.05 * t * Heaviside(t-10))
print(f"政策冲击效应: {policy_effect}")
#政策冲击效应: Heaviside(0.05*t*Heaviside(t - 10))
示例6. 数字电路 - 时钟信号
clock_signal = Heaviside(sin(2*pi*1*t)) # 简化表示
print(f"时钟信号: {clock_signal}")
#时钟信号: Heaviside(sin(2*pi*1*t))
示例7. 矩阵运算示例 - 多输入系统
matrix_input = [[-1, 0, 1], [2, -2, 0.5]]
matrix_result = heaviside(matrix_input)
print(f"输入矩阵: {matrix_input}")
print(f"阶跃输出: {matrix_result}")
#阶跃输出: [[0,0.5,1]
[1,0,1]]
示例8. 分段函数构造
构造分段函数: f(t) = 0 (t<0), 2t (0≤t<2), 4 (t≥2)
piecewise_func = heaviside(2*t * (Heaviside(t) - Heaviside(t-2)) + 4 * Heaviside(t-2))
print(f"分段函数: {piecewise_func}")
#分段函数: Heaviside(2*t*(Heaviside(t) - Heaviside(t - 2)) + 4*Heaviside(t - 2))
示例9. 滤波器设计 - 理想低通滤波器响应
频率响应的简化表示
lowpass_response = heaviside(Heaviside(1 - abs(f)))
print(f"理想低通响应: {lowpass_response}")
#理想低通响应: Heaviside(Heaviside(1 - Abs(f)))
示例10. 概率论 - 累积分布函数
某些分布在特定点的累积概率
cdf_example = heaviside(0.3 * Heaviside(x) + 0.7 * Heaviside(x-1))
print(f"累积分布示例: {cdf_example}")
#累积分布示例: Heaviside(0.3*Heaviside(x) + 0.7*Heaviside(x - 1))
数值计算验证
test_points = [-2, -1, 0, 0.5, 1, 2]
for point in test_points:
result = heaviside(str(point))
print(f"Heaviside({point}) = {result}")
#Heaviside(-2) = 0.0
#Heaviside(-1) = 0.0
#Heaviside(0) = 0.5
#Heaviside(0.5) = 1.0
#Heaviside(1) = 1.0
#Heaviside(2) = 1.0
符号计算验证
symbolic_tests = [x, x-1, x+2, sin(x)]
for test in symbolic_tests:
result = heaviside(test)
print(f"Heaviside({test}) = {result}")
#Heaviside(x) = Heaviside(x)
#Heaviside(x-1) = Heaviside(x - 1)
#Heaviside(x+2) = Heaviside(x + 2)
#Heaviside(sin(x)) = Heaviside(sin(x))
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def heaviside_step_func(input_str):
"""
对标MATLAB的heaviside函数实现
参数:
input_str (str): 输入表达式，支持以下格式：
- 标量数值（如 '0', '-5.3'）
- 符号表达式（如 'x', 't^2'）
- 矩阵（如 '[[1, 0], [-2, 3]]'）
返回:
SymPy表达式/矩阵 或错误信息
特性:
- x < 0 时返回0
- x = 0 时返回0.5
- x > 0 时返回1
- 支持符号计算
- 支持矩阵元素级运算
"""
try:
expr = sp.sympify(input_str)
result = None
def heaviside_scalar(x, half='half'):
"""
对标MATLAB的 heaviside 函数（处理单个数值标量）
参数:
x (float): 输入标量值
half (str): 'half'（默认）对应MATLAB的 heaviside(x, 'half')
'zero' 对应MATLAB的 heaviside(x)
返回:
float: Heaviside函数值
"""
if x > 0:
return 1.0
elif x < 0:
return 0.0
else: # x == 0
return 0.5 if half == 'half' else 0.0
# 处理标量和符号输入
if expr.is_number:
return heaviside_scalar(float(expr))
elif expr.free_symbols:
return sp.Heaviside(expr).evalf()
return f"无法处理的输入类型: {input_str}"
except Exception as e:
print(e)
# 示例测试
if __name__ == "__main__":
# 示例1：标量数值测试
print("示例1（标量正数）:", heaviside_step_func("5"))
# 1.0
print("示例1（标量零）:", heaviside_step_func("0"))
# 0.5
print("示例1（标量负数）:", heaviside_step_func("-3.2"))
# 0
# 示例2：符号运算测试
x = sp.symbols('x')
print("\n示例2（符号正数）:", heaviside_step_func("x"))
# Heaviside(x)
print("示例2（符号零）:", heaviside_step_func("Heaviside(0)"))
# 1.0
埃尔米特多项式
hermiteH(n,x)表示x点处的n次Hermite多项式。
n —— 多项式的次数, 非负整数|符号变量|符号表达式|符号函数|向量|矩阵
x —— 输入,数字|向量|矩阵|多维数组|符号数|符号变量|符号数组|符号函数|符号表达式
示例1. 量子力学 - 谐振子波函数
前几个能级的波函数（包含Hermite多项式）
for i in range(4):
hermite_expr = hermiteH(i, x)
print(f"n={i} Hermite多项式: H_{i}(x) = {hermite_expr}")
#n=0 Hermite多项式: H_0(x) = 1.0
#n=1 Hermite多项式: H_1(x) = 2.0*x
#n=2 Hermite多项式: H_2(x) = 4.0*x**2 - 2.0
#n=3 Hermite多项式: H_3(x) = 8.0*x**3 - 12.0*x
print("\n量子谐振子波函数形式: ψ_n(x) ∝ H_n(αx) * exp(-α²x²/2)")
#量子谐振子波函数形式: ψ_n(x) ∝ H_n(αx) * exp(-α²x²/2)
示例2. 概率论 - 高斯-埃尔米特积分
高斯-埃尔米特求积节点和权重与Hermite多项式零点相关
print("高斯-埃尔米特求积公式: ∫exp(-x²)f(x)dx ≈ Σ w_i f(x_i)")
print("其中x_i是Hermite多项式的零点")
#高斯-埃尔米特求积公式: ∫exp(-x²)f(x)dx ≈ Σ w_i f(x_i)
#其中x_i是Hermite多项式的零点
计算H_3(x)的零点
H3 = hermiteH(3, x)
zeros = solve(H3, x)
print(f"H_3(x)的零点: {zeros}")
#H_3(x)的零点: [-1.22474487139159, 0.0, 1.22474487139159]
示例3. 数值分析 - 逼近理论
用Hermite多项式展开函数
target_func = "exp(-x**2/2) * cos(x)"
print(f"目标函数: {target_func}")
print("可用Hermite多项式展开: f(x) = Σ c_n H_n(x)")
#目标函数: exp(-x**2/2) * cos(x)
#可用Hermite多项式展开: f(x) = Σ c_n H_n(x)
示例4. 光学 - 高斯光束模式
Hermite-Gaussian光束模式:
for m in range(3):
for n in range(3):
mode = f"H_{m}(√2 x/w) * H_{n}(√2 y/w) * exp(-(x²+y²)/w²)"
print(f"TEM_{m}{n}模式: {mode}")
#TEM_00模式: H_0(√2 x/w) * H_0(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_01模式: H_0(√2 x/w) * H_1(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_02模式: H_0(√2 x/w) * H_2(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_10模式: H_1(√2 x/w) * H_0(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_11模式: H_1(√2 x/w) * H_1(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_12模式: H_1(√2 x/w) * H_2(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_20模式: H_2(√2 x/w) * H_0(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_21模式: H_2(√2 x/w) * H_1(√2 y/w) * exp(-(x²+y²)/w²)
#TEM_22模式: H_2(√2 x/w) * H_2(√2 y/w) * exp(-(x²+y²)/w²)
示例5. 热传导 - 热方程解
print("某些热传导问题的解包含Hermite多项式")
heat_solution = "T(x,t) = Σ A_n H_n(x/√(4αt)) * exp(-x²/(4αt))"
print(f"热传导解形式: {heat_solution}")
#某些热传导问题的解包含Hermite多项式
#热传导解形式: T(x,t) = Σ A_n H_n(x/√(4αt)) * exp(-x²/(4αt))
示例6. 金融数学 - 期权定价
在某些随机波动率模型中，Hermite多项式用于展开概率密度函数
示例7. 信号处理 - 高斯导数滤波器
图像处理中的边缘检测使用高斯函数的导数，与Hermite多项式相关
for k in range(3):
gaussian_derivative = f"d^{k}/dx^{k} [exp(-x²/2)] = (-1)^{k} H_{k}(x) exp(-x²/2)"
print(f"{k}阶高斯导数: {gaussian_derivative}")
#0阶高斯导数: d^0/dx^0 [exp(-x²/2)] = (-1)^0 H_0(x) exp(-x²/2)
#1阶高斯导数: d^1/dx^1 [exp(-x²/2)] = (-1)^1 H_1(x) exp(-x²/2)
#2阶高斯导数: d^2/dx^2 [exp(-x²/2)] = (-1)^2 H_2(x) exp(-x²/2)
示例8. 矩阵运算示例
matrix_n = [[0, 1], [2, 3]]
matrix_x = [[1, 2], [0.5, -1]]
由于当前实现限制，我们分别计算每个元素
for i in range(2):
for j in range(2):
n_val = [0,1,2,3][i*2+j]
x_val = [1,2,0.5,-1][i*2+j]
result = hermiteH(n_val, x_val)
print(f"H_{n_val}({x_val}) = {result}")
#H_0(1) = 1.0
#H_1(2) = 4.000000000000001
#H_2(0.5) = -0.9999999999999998
#H_3(-1) = 3.9999999999999987
示例9. 数值验证和可视化准备
test_cases = [
(0, 0), # H₀(0) = 1
(1, 1), # H₁(1) = 2
(2, 2), # H₂(2) = 14
(3, 0.5), # H₃(0.5) = -9.5
]
for n_val, x_val in test_cases:
result = hermiteH(n_val, x_val)
print(f"H_{n_val}({x_val}) = {result}")
#H_0(0) = 1.0
#H_1(1) = 2.0000000000000004
#H_2(2) = 14.000000000000004
#H_3(0.5) = -5.000000000000001
示例10. 符号运算验证
symbolic_tests = [
(0, x),
(1, x),
(2, x),
(n, 0),
]
for test in symbolic_tests:
result = hermiteH(test)
print(f"H{test} = {result}")
#H(0, x) = 1.00000000000000
#H(1, x) = 2.0*x
#H(2, x) = 4.0*x**2 - 2.0
#H(n, 0) = 1.77245385090552*2.0**n/gamma(1/2 - n/2)
示例11. 特殊性质展示
正交性: ∫ H_m(x) H_n(x) exp(-x²) dx = √π 2ⁿ n! δ_mn
递推关系
H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
验证递推关系
H2 = hermiteH(2, x)
H3 = hermiteH(3, x)
H4_actual = hermiteH(4, x)
H4_recurrence = 2*x*H3 - 2*3*H2
print(f"H₄(x) 直接计算: {H4_actual}")
print(f"H₄(x) 递推计算: {H4_recurrence}")
print(f"验证递推关系: {sp.simplify(H4_actual - H4_recurrence) == 0}")
#H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
#H₄(x) 直接计算: 16.0*x**4 - 48.0*x**2 + 12.0
#H₄(x) 递推计算: -24.0*x**2 + 2*x*(8.0*x**3 - 12.0*x) + 12.0
#验证递推关系: True
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import eval_hermite
def hermite_polynomial(input_str):
"""
计算埃尔米特多项式，对标MATLAB的hermiteH函数。
参数:
input_str (str): 输入字符串，应为元组形式，如"(n, x)"，其中n和x可以是标量或矩阵。
返回:
SymPy表达式、矩阵、NaN或错误信息字符串。
示例:
>>> hermite_polynomial("(2, x)")
4*x**2 - 2
>>> hermite_polynomial("([[0, 1], [2, 3]], y)")
Matrix([
[ 1, 2*y],
[4*y**2 - 2, 8*y**3 - 12*y]])
>>> hermite_polynomial("(-1, x)")
nan
>>> hermite_polynomial("(n, x)")
hermite(n, x)
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
n = int(expr[0])
x = float(expr[0])
result = eval_hermite(n, x)
elif any(e.free_symbols for e in expr):
n, x = expr
if n.is_integer and n >= 0:
result = sp.hermite(*expr)
else:
error = True
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}。"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 示例1: 标量计算
print(hermite_polynomial("(2, 2)"))
# 14.000000000000004
# 示例2: 符号计算
print(hermite_polynomial("(3, x)"))
# 8*x**3 - 12*x
矩阵的黑森贝格形式
H = hess(A) 查找H,即矩阵A的黑森贝格形式.
[P,H] = hess(A) 生成黑森贝格矩阵 H 和酉矩阵 P，这样 A = P*H*P' 并且 P'*P = eye(size(A))。
黑森贝格矩阵包含第一个下对角线下面的零值.如果矩阵对称或者是埃尔米特矩阵,则形状为三对角.
该矩阵具有和原点相同的特征值,但显示它们所需计算更少.
A -- 输入,方阵
示例1. 特征值计算 - QR算法预处理
A1 = [[4, 1, 2], [1, 3, 1], [2, 1, 2]]
Q1, H1 = hess(A1)
print(f"Hessenberg形式 H:\n{H1}")
#Hessenberg形式 H:
#[[4.0, -2.23606797749979, 2.22044604925031e-16],
[-2.23606797749979, 3.0, 1.0],
[0, 1.0, 2.0]]
print(f"正交矩阵 Q:\n{Q1}")
#正交矩阵 Q:
#[[1.0, 0, 0],
[0, -0.447213595499958, -0.894427190999916],
[0, -0.894427190999916, 0.447213595499958]]
验证分解正确性: Q^T * A * Q = H
A_original = A1
verification = Q1.T * A_original * Q1
print(f"验证 Q^T * A * Q = H: {verification.equals(H1)}\n")
#验证 Q^T * A * Q = H: False
示例2. 控制系统 - 可控/可观测标准型
状态空间系统矩阵
A_sys = [[0.9, 0.1, 0.2], [0.3, 0.8, 0.1], [0.1, 0.2, 0.7]]
Q_sys, H_sys = hess(A_sys)
print(f"系统矩阵 A:\n{A_sys}")
#系统矩阵 A:
#[[0.9, 0.1, 0.2],
[0.3, 0.8, 0.1],
[0.1, 0.2, 0.7]]
print(f"Hessenberg形式 (控制器型):\n{H_sys}")
#[[0.9, -0.158113883008419, 0.158113883008419],
[-0.316227766016838, 0.88, -0.04],
[0, -0.14, 0.62]]
#在控制器设计中，Hessenberg形式便于实现状态反馈
示例3. 数值分析 - 矩阵幂计算
A_power = [[2, 1, 0], [1, 3, 1], [0, 1, 2]]
Q_power, H_power = hess(A_power)
print(f"原始矩阵:\n{A_power}")
print(f"Hessenberg形式:\n{H_power}")
#原始矩阵:
#[[2, 1, 0],
[1, 3, 1],
[0, 1, 2]]
#Hessenberg形式:
#[[2.0, 1.0, 0],
[1.0, 3.0, 1.0],
[0, 1.0, 2.0]]
#计算 A^n = Q * H^n * Q^T，其中H^n计算更高效
示例4. 结构动力学 - 质量-弹簧系统
刚度矩阵（对称正定）
K_matrix = [[2, -1, 0], [-1, 2, -1], [0, -1, 1]]
Q_k, H_k = hess(K_matrix)
print(f"刚度矩阵 K:\n{K_matrix}")
print(f"Hessenberg形式:\n{H_k}")
#刚度矩阵 K:
#[[2, -1, 0],
[-1, 2, -1],
[0, -1, 1]]
#Hessenberg形式:
#[[2.0, -1.0, 0],
[-1.0, 2.0, -1.0],
[0, -1.0, 1.0]]
#用于求解特征值问题 Kx = λMx
示例5. 量子力学 - 哈密顿量对角化
简单的三能级系统哈密顿量
H_quantum = [[1, 0.5, 0], [0.5, 2, 0.3], [0, 0.3, 1.5]]
Q_q, H_q = hess(H_quantum)
print(f"哈密顿量 H:\n{H_quantum}")
print(f"Hessenberg形式:\n{H_q}")
print("作为QR算法迭代的基础，用于求解能量本征值\n")
#哈密顿量 H:
#[[1, 0.5, 0],
[0.5, 2, 0.3],
[0, 0.3, 1.5]]
#Hessenberg形式:
#[[1.0, 0.5, 0],
[0.5, 2.0, 0.3],
[0, 0.3, 1.5]]
#作为QR算法迭代的基础，用于求解能量本征值
示例6. 图像处理 - 协方差矩阵分析
图像块的协方差矩阵
cov_matrix = [[100, 50, 25], [50, 80, 40], [25, 40, 60]]
Q_cov, H_cov = hess(cov_matrix)
print(f"协方差矩阵:\n{cov_matrix}")
print(f"Hessenberg形式:\n{H_cov}")
#协方差矩阵:
#[[100, 50, 25],
[50, 80, 40],
[25, 40, 60]]
#Hessenberg形式:
#[[100.0, -55.9016994374947, -3.55271367880050e-15],
[-55.9016994374947, 108.0, -16.0],
[0, -16.0, 32.0]]
#用于PCA特征值计算，Hessenberg形式加速收敛
示例7. 网络分析 - 图邻接矩阵
简单图的邻接矩阵
adj_matrix = [[0, 1, 1, 0], [1, 0, 1, 1], [1, 1, 0, 1], [0, 1, 1, 0]]
Q_adj, H_adj = hess(adj_matrix)
print(f"邻接矩阵:\n{adj_matrix}")
print(f"Hessenberg形式:\n{H_adj}")
#邻接矩阵:
#[[0, 1, 1, 0],
[1, 0, 1, 1],
[1, 1, 0, 1],
[0, 1, 1, 0]]
#Hessenberg形式:
#[[0, -1.41421356237309, 0, 0],
[-1.41421356237310, 0.999999999999999, -1.41421356237309, 2.22044604925031e-16],
[0, -1.41421356237309, 0, 0], [0, 0, 0, -1.00000000000000]]
#用于计算图的谱性质（特征值）
示例8. 金融数学 - 资产相关矩阵
corr_matrix = [[1, 0.7, 0.3], [0.7, 1, 0.5], [0.3, 0.5, 1]]
Q_corr, H_corr = hess(corr_matrix)
print(f"相关矩阵:\n{corr_matrix}")
print(f"Hessenberg形式:\n{H_corr}")
print("用于风险分析和投资组合优化\n")
#相关矩阵:
#[[1, 0.7, 0.3],
[0.7, 1, 0.5],
[0.3, 0.5, 1]]
#Hessenberg形式:
#[[1.0, -0.761577310586391, 0],
[-0.761577310586391, 1.36206896551724, -0.344827586206897],
[0, -0.344827586206896, 0.637931034482759]])
#用于风险分析和投资组合优化
示例9. 机械工程 - 有限元刚度矩阵
简化的有限元刚度矩阵
fem_matrix = [[3, -1, 0, 0], [-1, 2, -1, 0], [0, -1, 2, -1], [0, 0, -1, 1]]
Q_fem, H_fem = hess(fem_matrix)
print(f"有限元刚度矩阵:\n{sp.Matrix(eval(fem_matrix))}")
print(f"Hessenberg形式:\n{H_fem}")
#有限元刚度矩阵:
#[[3, -1, 0, 0],
[-1, 2, -1, 0],
[0, -1, 2, -1],
[0, 0, -1, 1]]
#Hessenberg形式:
#[[3.0, -1.0, 0, 0],
[-1.0, 2.0, -1.0, 0],
[0, -1.0, 2.0, -1.0],
[0, 0, -1.0, 1.0]]
#用于结构振动分析和模态提取
示例10. 数值稳定性分析 - 病态矩阵处理
Hilbert矩阵（著名的病态矩阵）
hilbert_3 = [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]]
Q_hilb, H_hilb = hess(hilbert_3)
print(f"Hilbert矩阵:\n{hilbert_3}")
print(f"Hessenberg形式:\n{H_hilb}")
#Hilbert矩阵:
#[[1, 0.5, 0.333333333333333],
[0.5, 0.333333333333333, 0.25],
[0.333333333333333, 0.25, 0.2]]
#Hessenberg形式:
#[[1.0, -0.600925212577331, 5.55111512312578e-17],
[-0.600925212577332, 0.523076923076923, -0.0346153846153847],
[0, -0.0346153846153846, 0.0102564102564103]]
#Hessenberg分解比直接特征值分解对病态矩阵更稳定
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.linalg import hessenberg
import numpy as np
def hessenberg_matrix(input_str):
"""
将输入的矩阵字符串转换为Hessenberg形式，返回Hessenberg矩阵H和变换矩阵Q，对标MATLAB的hess函数。
参数:
input_str (str): 表示矩阵的字符串，例如 "[[1, 2], [3, 4]]"。
返回:
tuple: 包含Hessenberg矩阵H和变换矩阵Q的元组（均为SymPy矩阵），
或错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为元组（不支持）
if isinstance(expr, tuple):
error = True
elif isinstance(expr, list):
A = sp.Matrix(expr)
if A.is_square:
# 转换为numpy数组进行计算
N_A = np.array(A, dtype=float)
H_np, Q_np = hessenberg(N_A, calc_q=True)
# 转换回SymPy矩阵
H = sp.Matrix(H_np)
Q = sp.Matrix(Q_np)
result = (Q, H)
else:
error = True
else:
error = True
return result if not error else f"输入错误: 无法将'{input_str}'转换为矩阵。"
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：正确的方阵输入
print("示例1:")
input_1 = "[[1, 2, 3], [4, 5, 6], [7, 8, 9]]"
result_1 = hessenberg_matrix(input_1)
if isinstance(result_1, tuple):
Q, H = result_1
print("Hessenberg矩阵H:")
print(H)
# Matrix([[1.00000000000000, -3.59700730308704, -0.248069469178417],
# [-8.06225774829855, 14.0461538461538, 2.83076923076923],
# [0, 0.830769230769230, -0.0461538461538478]])
print("\n变换矩阵Q:")
print(Q)
# Matrix([[1.00000000000000, 0, 0],
# [0, -0.496138938356834, -0.868243142124459],
# [0, -0.868243142124459, 0.496138938356834]])
else:
print(result_1)
黑塞对角矩阵
hessdiag(f,v,points) 返回Hessian矩阵的对角元素（第二部分向量）.
如果不指定 v, 则 hessian(f) 查找标量函数 f 相对于由 f 中找到的所有符号变量构成的向量的 Hessian 矩阵。此向量中变量的顺序由 symvar 定义.
f — 标量函数,符号表达式,符号函数
v — 查找 Hessian 对角矩阵的向量,符号向量
points — 数值向量
示例1. 机器学习损失函数的曲率分析
线性回归的均方误差
result1 = hessian((y - (w1*x1 + w2*x2 + b))**2)
print(f"线性回归损失函数: {input_str1}")
print(f"对角黑塞矩阵:\n{result1}")
#线性回归损失函数: (y - (w1*x1 + w2*x2 + b))**2
#对角黑塞矩阵:
#[[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w1*(x1 - 0.001) - w2*x2 + y)**2 + 1000000.0*(-b - w1*(x1 + 0.001) - w2*x2 + y)**2],
[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w1*x1 - w2*(x2 - 0.001) + y)**2 + 1000000.0*(-b - w1*x1 - w2*(x2 + 0.001) + y)**2],
[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w1*x1 - x2*(w2 - 0.001) + y)**2 + 1000000.0*(-b - w1*x1 - x2*(w2 + 0.001) + y)**2],
[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w1*x1 - w2*x2 + y - 0.001)**2 + 1000000.0*(-b - w1*x1 - w2*x2 + y + 0.001)**2],
[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w1*x1 - w2*x2 + y - 0.001)**2 + 1000000.0*(-b - w1*x1 - w2*x2 + y + 0.001)**2],
[-2000000.0*(-b - w1*x1 - w2*x2 + y)**2 + 1000000.0*(-b - w2*x2 - x1*(w1 - 0.001) + y)**2 + 1000000.0*(-b - w2*x2 - x1*(w1 + 0.001) + y)**2]]
在特定点计算
result1_point = hessian((y - (w1*x1 + w2*x2 + b))**2, [w1, w2, b], [0.5, -0.3, 1.0])
print(f"在点 [w1=0.5, w2=-0.3, b=1.0] 的值:\n{result1_point}")
#在点 [w1=0.5, w2=-0.3, b=1.0] 的值:
#[[1000000.0*(-0.501*x1 + 0.3*x2 + y - 1.0)**2 - 2000000.0*(-0.5*x1 + 0.3*x2 + y - 1.0)**2 + 1000000.0*(-0.499*x1 + 0.3*x2 + y - 1.0)**2],
[1000000.0*(-0.5*x1 + 0.299*x2 + y - 1.0)**2 - 2000000.0*(-0.5*x1 + 0.3*x2 + y - 1.0)**2 + 1000000.0*(-0.5*x1 + 0.301*x2 + y - 1.0)**2],
[-2000000.0*(-0.5*x1 + 0.3*x2 + y - 1.0)**2 + 1000000.0*(-0.5*x1 + 0.3*x2 + y - 0.999)**2 + 1002001.0*(-0.4995004995005*x1 + 0.2997002997003*x2 + 0.999000999000999*y - 1.0)**2]]
示例2. 物理势能函数的稳定性分析
谐振子势能
result2 = hessian(0.5*k*x**2)
print(f"谐振子势能: {input_str2}")
print(f"对角黑塞矩阵:\n{result2}")
#谐振子势能: 0.5*k*x**2
#对角黑塞矩阵:
#[[-1000000.0*k*x**2 + 500000.0*k*(x - 0.001)**2 + 500000.0*k*(x + 0.001)**2],
[-1000000.0*k*x**2 + 500000.0*x**2*(k - 0.001) + 500000.0*x**2*(k + 0.001)]]
双井势能
result3 = hessian(x**4 - 2*x**2)
print(f"双井势能: {input_str3}")
print(f"对角黑塞矩阵:\n{result3}")
#双井势能: x**4 - 2*x**2
#对角黑塞矩阵:
#[[-2000000.0*x**4 + 4000000.0*x**2 + 1000000.0*(x - 0.001)**4 - 2000000.0*(x - 0.001)**2 + 1000000.0*(x + 0.001)**4 - 2000000.0*(x + 0.001)**2]]
在临界点分析
input_str3_points = [
((x**4 - 2*x**2), [x], [-1.0]), # 局部极小值点
((x**4 - 2*x**2), [x], [0.0]), # 局部极大值点
((x**4 - 2*x**2), [x], [1.0]) # 局部极小值点
]
for i, input_str in enumerate(input_str3_points):
result = hessian(input_str)
points = [-1.0, 0.0, 1.0][i]
print(f"在 x={points} 点的曲率: {result}")
#在 x=-1.0 点的曲率: [[8.00000200024806]]
#在 x=0.0 点的曲率: [[-3.99999800000000]]
#在 x=1.0 点的曲率: [[8.00000200024806]]
示例3. 经济学效用函数的凸性分析
Cobb-Douglas 效用函数
result4 = hessian(x**0.3 * y**0.7)
print(f"Cobb-Douglas 效用函数: {input_str4}")
print(f"对角黑塞矩阵:\n{result4}")
#Cobb-Douglas 效用函数: x**0.3 * y**0.7
#对角黑塞矩阵:
#[[-2000000.0*x**0.3*y**0.7 + 1000000.0*x**0.3*(y - 0.001)**0.7 + 1000000.0*x**0.3*(y + 0.001)**0.7],
[-2000000.0*x**0.3*y**0.7 + 1000000.0*y**0.7*(x - 0.001)**0.3 + 1000000.0*y**0.7*(x + 0.001)**0.3]]
在特定消费束计算
result4_point = hessian((x**0.3 * y**0.7), [x, y], [10, 20])
print(f"在消费束 [x=10, y=20] 的曲率:\n{result4_point}")
#在消费束 [x=10, y=20] 的曲率:
#[[-0.0341145979580806],
[-0.00852864980697632]]
示例4. 工程优化问题的曲率信息
梁的弯曲能量
result5 = hessian(E*I*(d2w_dx2)**2)
print(f"梁弯曲能量: {input_str5}")
print(f"对角黑塞矩阵:\n{result5}")
#梁弯曲能量: E*I*(d2w_dx2)**2
#对角黑塞矩阵:
#[[-2000000.0*E*I*d2w_dx2**2 + 1000000.0*E*I*(d2w_dx2 - 0.001)**2 + 1000000.0*E*I*(d2w_dx2 + 0.001)**2]]
材料应力函数
result6 = hessian(sigma_x**2 + sigma_y**2 - sigma_x*sigma_y + 3*tau_xy**2)
print(f"冯·米塞斯应力准则: {input_str6}")
print(f"对角黑塞矩阵:\n{result6}")
#冯·米塞斯应力准则: sigma_x**2 + sigma_y**2 - sigma_x*sigma_y + 3*tau_xy**2
#对角黑塞矩阵:
#[[-6000000.0*tau_xy**2 + 3000000.0*(tau_xy - 0.001)**2 + 3000000.0*(tau_xy + 0.001)**2],
[-2000000.0*sigma_x**2 + 2000000.0*sigma_x*sigma_y - 1000000.0*sigma_y*(sigma_x - 0.001) - 1000000.0*sigma_y*(sigma_x + 0.001) + 1000000.0*(sigma_x - 0.001)**2 + 1000000.0*(sigma_x + 0.001)**2],
[2000000.0*sigma_x*sigma_y - 1000000.0*sigma_x*(sigma_y - 0.001) - 1000000.0*sigma_x*(sigma_y + 0.001) - 2000000.0*sigma_y**2 + 1000000.0*(sigma_y - 0.001)**2 + 1000000.0*(sigma_y + 0.001)**2]]
示例5. 多元复杂函数的二阶特性
三维 Rosenbrock 函数
result7 = hessian((1-x)**2 + 100*(y-x**2)**2 + (1-y)**2 + 100*(z-y**2)**2)
print(f"三维 Rosenbrock 函数: {input_str7}")
print(f"对角黑塞矩阵:\n{result7}")
#三维 Rosenbrock 函数: (1-x)**2 + 100*(y-x**2)**2 + (1-y)**2 + 100*(z-y**2)**2
#对角黑塞矩阵:
#[[-200000000.0*(-y**2 + z)**2 + 100000000.0*(-y**2 + z - 0.001)**2 + 100000000.0*(-y**2 + z + 0.001)**2],
[1000000.0*(0.999 - y)**2 - 2000000.0*(1 - y)**2 + 1002001.0*(1 - 0.999000999000999*y)**2 - 200000000.0*(-x**2 + y)**2 - 200000000.0*(-y**2 + z)**2 + 100000000.0*(z - (y - 0.001)**2)**2 + 100000000.0*(z - (y + 0.001)**2)**2 + 100000000.0*(-x**2 + y - 0.001)**2 + 100000000.0*(-x**2 + y + 0.001)**2], [1000000.0*(0.999 - x)**2 - 2000000.0*(1 - x)**2 + 1002001.0*(1 - 0.999000999000999*x)**2 - 200000000.0*(-x**2 + y)**2 + 100000000.0*(y - (x - 0.001)**2)**2 + 100000000.0*(y - (x + 0.001)**2)**2]]
在最小值点附近
result7_point = hessian((1-x)**2 + 100*(y-x**2)**2 + (1-y)**2 + 100*(z-y**2)**2, [x, y, z], [1.0, 1.0, 1.0])
print(f"在点 [1,1,1] 的曲率:\n{result7_point}")
#在点 [1,1,1] 的曲率:
#[[802.000199999866],
[1002.00019999984],
[199.999999999971]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def get_result_value(symbols_set, points=None, result=None):
if points is not None and result is not None:
if len(points) != len(symbols_set):
raise ValueError(f"点数（{len(points)}）与符号数（{len(symbols_set)}）不匹配")
subs_dict = dict(zip(symbols_set, points))
try:
return result.subs(subs_dict).evalf()
except Exception as e:
return result.subs(subs_dict)
return result
def hessenberg_diag_matrix(input_str):
"""
计算给定表达式的黑塞矩阵对角线元素的近似值。
参数:
input_str (str): 输入的数学表达式字符串，可以包含变量和数值列表。
返回:
sp.Matrix or str: 黑塞矩阵对角线元素的近似值矩阵，如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 有限差分步长
h = 0.001
def remove_last_from_tuple(tup):
"""
移除元组的最后一个元素。
参数:
tup (tuple): 待处理的元组
返回:
tuple: 移除最后一个元素后的元组
"""
return tup[:-1] if len(tup) > 1 else tup
# 计算黑塞矩阵对角线元素近似值
def hessian_diag_approx(func, vars, h):
"""
计算黑塞矩阵对角线元素的近似值。
参数:
func: 待求导的函数
vars: 自变量列表
h: 有限差分步长
返回:
list: 黑塞矩阵对角线元素的近似值列表
"""
diag_elements = []
for var in vars:
f_plus = func.subs(var, var + h)
f_minus = func.subs(var, var - h)
df2 = (f_plus - 2 * func + f_minus) / (h ** 2)
diag_elements.append(df2)
return diag_elements
# 这里判断是否需要数值求导，判断最后一项是否是数值列表比如[1,2]代入最后的公式result
points = None
if isinstance(expr, tuple) and isinstance(expr[-1], (list, tuple)) and all(item.is_number for item in expr[-1]):
points = expr[-1]
expr = remove_last_from_tuple(expr)
# 初始化函数和自变量
f = None
vars = None
if isinstance(expr, tuple) and expr[0].free_symbols and isinstance(expr[1], (list, tuple)):
f = expr[0]
vars = expr[1]
elif expr.free_symbols:
f = expr
vars = list(expr.free_symbols)
else:
error = True
if error:
return f"输入错误:{input_str}"
else:
result = hessian_diag_approx(f, vars, h)
result = sp.Matrix(result)
result = get_result_value(vars, points, result)
return result
except Exception as e:
return f"错误:{e}"
def main():
# 示范代码
# 示例 1: 仅表达式，无数值点
input_str1 = "x**2 + y**2"
result1 = hessenberg_diag_matrix(input_str1)
print(f"输入: {input_str1}, 结果: {result1}")
# Matrix([[-2000000.0*x**2 + 1000000.0*(x - 0.001)**2 + 1000000.0*(x + 0.001)**2],
# [-2000000.0*y**2 + 1000000.0*(y - 0.001)**2 + 1000000.0*(y + 0.001)**2]])
# 示例 2: 表达式和数值点
input_str2 = "(x**2 + y**2, [x, y], [1, 2])"
result2 = hessenberg_diag_matrix(input_str2)
print(f"输入: {input_str2}, 结果: {result2}")
# Matrix([[1.99999999953434],
# [2.00000000000000]])
if __name__ == "__main__":
main()
黑塞矩阵
hessian(f,v,points) 查找符号标量函数 f 相对于笛卡尔坐标中的向量 v 的 Hessian 矩阵.
如果不指定 v, 则 hessian(f) 查找标量函数 f 相对于由 f 中找到的所有符号变量构成的向量的 Hessian 矩阵。此向量中变量的顺序由 symvar 定义.
f — 标量函数,符号表达式,符号函数
v — 查找 Hessian 矩阵的向量,符号向量
points — 数值向量
示例1. 机器学习中的损失函数分析
线性回归的MSE损失
mse_hessian = hessian((y - (a*x + b))**2)
print(f"Hessian矩阵:\n{mse_hessian}")
#Hessian矩阵:
#[[2*x**2, 2*x, 4*a*x + 2*b - 2*y, -2*x],
[2*x, 2, 2*a, -2],
[4*a*x + 2*b - 2*y, 2*a, 2*a**2, -2*a],
[-2*x, -2, -2*a, 2]]
示例2. 物理学中的势能函数
谐振子势能
harmonic_hessian = hessian(0.5*k*x**2)
print(f"Hessian矩阵:\n{harmonic_hessian}")
#Hessian矩阵:
#[[0, 1.0*x],
[1.0*x, 1.0*k]]
重力势能
gravity_hessian = hessian(-m/(x**2 + y**2)**0.5)
print(f"Hessian矩阵:\n{gravity_hessian}")
#Hessian矩阵:
#[[0, 1.0*x/(x**2 + y**2)**1.5, 1.0*y/(x**2 + y**2)**1.5],
[1.0*x/(x**2 + y**2)**1.5, -3.0*m*x**2/(x**2 + y**2)**2.5 + 1.0*m/(x**2 + y**2)**1.5, -3.0*m*x*y/(x**2 + y**2)**2.5],
[1.0*y/(x**2 + y**2)**1.5, -3.0*m*x*y/(x**2 + y**2)**2.5, -3.0*m*y**2/(x**2 + y**2)**2.5 + 1.0*m/(x**2 + y**2)**1.5]]
在特定点计算
gravity_at_point = hessian(-m/(x**2 + y**2)**0.5, [x, y, m], [1, 0, 1])
print(f"在点 (1,0,1) 的Hessian:\n{gravity_at_point}")
#在点 (1,0,1) 的Hessian:
#[[-2.0, 0, 1.0],
[0, 1.0, 0],
[1.0, 0, 0]]
示例3. 经济学与金融应用
Cobb-Douglas生产函数
cd_hessian = hessian(A * x**a * y**(1-a))
print(f"Hessian矩阵:\n{cd_hessian}")
#Hessian矩阵:
#[[0, x**a*y**(1 - a)*log(x) - x**a*y**(1 - a)*log(y), a*x**a*y**(1 - a)/x, x**a*y**(1 - a)*(1 - a)/y],
[x**a*y**(1 - a)*log(x) - x**a*y**(1 - a)*log(y), A*x**a*y**(1 - a)*log(x)**2 - 2*A*x**a*y**(1 - a)*log(x)*log(y) + A*x**a*y**(1 - a)*log(y)**2, A*a*x**a*y**(1 - a)*log(x)/x - A*a*x**a*y**(1 - a)*log(y)/x + A*x**a*y**(1 - a)/x, A*x**a*y**(1 - a)*(1 - a)*log(x)/y - A*x**a*y**(1 - a)*(1 - a)*log(y)/y - A*x**a*y**(1 - a)/y],
[a*x**a*y**(1 - a)/x, A*a*x**a*y**(1 - a)*log(x)/x - A*a*x**a*y**(1 - a)*log(y)/x + A*x**a*y**(1 - a)/x, A*a**2*x**a*y**(1 - a)/x**2 - A*a*x**a*y**(1 - a)/x**2, A*a*x**a*y**(1 - a)*(1 - a)/(x*y)], [x**a*y**(1 - a)*(1 - a)/y, A*x**a*y**(1 - a)*(1 - a)*log(x)/y - A*x**a*y**(1 - a)*(1 - a)*log(y)/y - A*x**a*y**(1 - a)/y, A*a*x**a*y**(1 - a)*(1 - a)/(x*y), A*x**a*y**(1 - a)*(1 - a)**2/y**2 - A*x**a*y**(1 - a)*(1 - a)/y**2]]
示例4. 工程学中的优化问题
梁的弯曲能量
beam_hessian = hessian(0.5*E*I*(d2w_dx2)**2)
print(f"Hessian矩阵:\n{beam_hessian}")
#Hessian矩阵:
#[[1.0*E*I]]
热传导方程
heat_hessian = hessian(k*(dT_dx**2 + dT_dy**2))
print(f"Hessian矩阵:\n{heat_hessian}")
#Hessian矩阵:
#[[2*k, 0, 2*dT_dx],
[0, 2*k, 2*dT_dy],
[2*dT_dx, 2*dT_dy, 0]]
示例5. 化学与材料科学
Lennard-Jones势能
lj_hessian = hessian(4*epsilon*((sigma/r)**12 - (sigma/r)**6))
print(f"Hessian矩阵:\n{lj_hessian}")
#Hessian矩阵:
#[[0, 24*sigma**6/r**7 - 48*sigma**12/r**13, -24*sigma**5/r**6 + 48*sigma**11/r**12],
[24*sigma**6/r**7 - 48*sigma**12/r**13, 4*epsilon*(-42*sigma**6/r**8 + 156*sigma**12/r**14), 4*epsilon*(36*sigma**5/r**7 - 144*sigma**11/r**13)],
[-24*sigma**5/r**6 + 48*sigma**11/r**12, 4*epsilon*(36*sigma**5/r**7 - 144*sigma**11/r**13), 4*epsilon*(-30*sigma**4/r**6 + 132*sigma**10/r**12)]]
在平衡位置计算
lj_at_equilibrium = hessian(
(4*epsilon*((sigma/r)**12 - (sigma/r)**6), [r, epsilon, sigma], [2**(1/6), 1, 1]))
print(f"在平衡位置的Hessian:\n{lj_at_equilibrium}")
#在平衡位置的Hessian:
#[[57.1464378708552, 0, -64.1447077061044],
[0, 0, 0],
[-64.1447077061044, 0, 72.0000000000000]]
示例6. 计算机视觉中的图像处理
图像强度函数的曲率
image_hessian = hessian(I_xx + I_yy)
print(f"Hessian矩阵:\n{image_hessian}")
#Hessian矩阵:
#[[0, 0],
[0, 0]]
示例7. 机器人运动学
机械臂重力势能
robot_hessian = hessian(m1*g*l1*cos(theta1) + m2*g*(l1*cos(theta1) + l2*cos(theta2)))
print(f"Hessian矩阵:\n{robot_hessian}")
#Hessian矩阵:
#[[0, m1*cos(theta1) + m2*cos(theta1), m2*cos(theta2), l1*cos(theta1), l1*cos(theta1) + l2*cos(theta2), -l1*m1*sin(theta1) - l1*m2*sin(theta1), -l2*m2*sin(theta2)],
[m1*cos(theta1) + m2*cos(theta1), 0, 0, g*cos(theta1), g*cos(theta1), -g*m1*sin(theta1) - g*m2*sin(theta1), 0],
[m2*cos(theta2), 0, 0, 0, g*cos(theta2), 0, -g*m2*sin(theta2)],
[l1*cos(theta1), g*cos(theta1), 0, 0, 0, -g*l1*sin(theta1), 0],
[l1*cos(theta1) + l2*cos(theta2), g*cos(theta1), g*cos(theta2), 0, 0, -g*l1*sin(theta1), -g*l2*sin(theta2)],
[-l1*m1*sin(theta1) - l1*m2*sin(theta1), -g*m1*sin(theta1) - g*m2*sin(theta1), 0, -g*l1*sin(theta1), -g*l1*sin(theta1), -g*l1*m1*cos(theta1) - g*l1*m2*cos(theta1), 0],
[-l2*m2*sin(theta2), 0, -g*m2*sin(theta2), 0, -g*l2*sin(theta2), 0, -g*l2*m2*cos(theta2)]]
示例8. 统计学中的概率分布
多元正态分布的负对数似然
normal_hessian = hessian(0.5*((x-mu_x)**2/sigma_x**2 + (y-mu_y)**2/sigma_y**2))
print(f"Hessian矩阵:\n{normal_hessian}")
#Hessian矩阵:
#[[1.0/sigma_x**2, 0, -1.0*(2*mu_x - 2*x)/sigma_x**3, 0, -1.0/sigma_x**2, 0],
[0, 1.0/sigma_y**2, 0, -1.0*(2*mu_y - 2*y)/sigma_y**3, 0, -1.0/sigma_y**2],
[-1.0*(2*mu_x - 2*x)/sigma_x**3, 0, 3.0*(-mu_x + x)**2/sigma_x**4, 0, -1.0*(-2*mu_x + 2*x)/sigma_x**3, 0],
[0, -1.0*(2*mu_y - 2*y)/sigma_y**3, 0, 3.0*(-mu_y + y)**2/sigma_y**4, 0, -1.0*(-2*mu_y + 2*y)/sigma_y**3],
[-1.0/sigma_x**2, 0, -1.0*(-2*mu_x + 2*x)/sigma_x**3, 0, 1.0/sigma_x**2, 0],
[0, -1.0/sigma_y**2, 0, -1.0*(-2*mu_y + 2*y)/sigma_y**3, 0, 1.0/sigma_y**2]]
示例9. 控制系统Lyapunov函数
二次型Lyapunov函数
lyapunov_hessian = hessian(x**2 + y**2 + z**2)
print(f"Hessian矩阵:\n{lyapunov_hessian}")
#Hessian矩阵:
#[[2, 0, 0],
[0, 2, 0],
[0, 0, 2]]
示例10. 经典优化测试函数
Rosenbrock函数（香蕉函数）
rosenbrock_hessian = hessian((1-x)**2 + 100*(y-x**2)**2)
print(f"Hessian矩阵:\n{rosenbrock_hessian}")
#Hessian矩阵:
#[[1200*x**2 - 400*y + 2, -400*x],
[-400*x, 200]]
在最小值点
rosen_at_min = hessian((1-x)**2 + 100*(y-x**2)**2, [x, y], [1, 1])
print(f"在最小值点 (1,1) 的Hessian:\n{rosen_at_min}")
#在最小值点 (1,1) 的Hessian:
#[[802.0, -400.0],
[-400.0, 200.0]]
Himmelblau函数
himmelblau_hessian = hessian((x**2 + y - 11)**2 + (x + y**2 - 7)**2)
print(f"Hessian矩阵:\n{himmelblau_hessian}")
#Hessian矩阵:
#[[12*x**2 + 4*y - 42, 4*x + 4*y],
[4*x + 4*y, 4*x + 12*y**2 - 26]]
在一个极值点
himmel_at_point = hessian((x**2 + y - 11)**2 + (x + y**2 - 7)**2, [x, y], [3, 2])
print(f"在点 (3,2) 的Hessian:\n{himmel_at_point}")
#在点 (3,2) 的Hessian:
#[[74.0, 20.0],
[20.0, 34.0]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def get_result_value(symbols_set, points=None, result=None):
"""
将符号表达式代入具体数值点求值。
参数:
symbols_set: 符号集合（顺序敏感的列表）
points: 数值点列表，与符号一一对应
result: 待代入的表达式
返回:
代入后的数值结果或符号表达式
异常:
ValueError: 当点数与符号数量不匹配时抛出
"""
if points is not None and result is not None:
if len(points) != len(symbols_set):
raise ValueError(f"提供的点数（{len(points)}）与符号数量（{len(symbols_set)}）不匹配。")
subs_dict = {symbol: value for symbol, value in zip(symbols_set, points)}
try:
# 尝试数值求值（适用于全数值表达式）
result = result.subs(subs_dict).evalf()
except Exception as e:
# 非数值情况直接进行符号替换
result = result.subs(subs_dict)
return result
def hessian_matrix(input_str):
"""
对标MATLAB的hessian函数，计算符号标量函数的Hessian矩阵
参数:
input_str: 字符串形式的输入，可以是以下形式：
- 单个表达式（如 "x**2 + y**3"）
- 元组表达式，包含函数和变量列表（如 "(f, [x, y])"）
- 元组表达式，包含函数、变量列表和数值点（如 "(f, [x, y], [1, 2])"）
返回:
SymPy矩阵对象或错误信息字符串
示例:
>>> hessian_matrix("x**2 + y**3")
Matrix([[2, 0], [0, 6*y]])
>>> hessian_matrix("(x**2 + y**3, [x, y], [1, 2])")
Matrix([[2, 0], [0, 12]])
"""
try:
# 将输入字符串转换为SymPy表达式（可能是元组、矩阵或普通表达式）
expr = sp.sympify(input_str)
error = False
result = None
def remove_last_from_tuple(tup):
"""辅助函数：当输入元组包含数值点时，移除最后一个元素"""
return tup[:-1] if len(tup) > 1 else tup
# 检查输入是否为带数值点的元组结构（例如 (f, [x,y], [1,2])）
if isinstance(expr, tuple) and len(expr) >= 2 and isinstance(expr[-1], list):
# 提取数值点并验证是否为纯数字
if all(item.is_number for item in expr[-1]):
points = expr[-1]
expr = remove_last_from_tuple(expr) # 移除数值点部分
else:
points = None
else:
points = None
# 解析符号变量列表（优先使用用户显式提供的变量列表）
if isinstance(expr, tuple) and len(expr) == 2 and isinstance(expr[1], list):
# 输入格式为 (函数, [变量列表])
f = expr[0]
symbols = expr[1]
result = sp.hessian(f, symbols)
elif hasattr(expr, 'free_symbols'):
# 输入为单个表达式，自动提取自由符号（注意顺序可能不固定！）
symbols = list(expr.free_symbols)
symbols.sort(key=lambda x: x.name) # 按字母顺序排序以保证一致性
result = sp.hessian(expr, symbols)
else:
raise ValueError("无法识别的输入结构")
# 如果存在数值点，进行符号替换
if points is not None:
result = get_result_value(symbols, points, result)
# 确保返回的是矩阵类型
return sp.Matrix(result) if result is not None else None
except Exception as e:
return f"错误:{e}"
if __name__ == "__main__":
# ----------------- 示例演示 -------------------
# 示例1：基本用法（自动提取变量）
print("示例1:", hessian_matrix("x*y+2*z*x"), end='\n\n')
# Matrix([[0, 1, 2],
# [1, 0, 0],
# [2, 0, 0]])
# 示例2：指定变量顺序（推荐方式）
x, y = sp.symbols('x y')
print("示例2:", hessian_matrix(f"( {x ** 2 + y ** 3}, [{x}, {y}] )"), end='\n\n')
# Matrix([[2, 0],
# [0, 6*y]])
# 示例3：带数值代入
print("示例3:", hessian_matrix(f"( {x ** 2 + y ** 3}, [{x}, {y}], [1, 2] )"), end='\n\n')
# Matrix([[2.00000000000000, 0],
# [0, 12.0000000000000]])
# 示例4：多变量函数
z = sp.symbols('z')
print("示例4:", hessian_matrix(f"x*y + {z}**2"), end='\n\n')
# Matrix([[0, 1, 0],
# [1, 0, 0],
# [0, 0, 2]])
希尔伯特矩阵
H = hilb(n) 返回阶数为n的希尔伯特矩阵.希尔伯特矩阵是病态矩阵的典型示例.希尔伯特矩阵的元素由 H(i,j) = 1/(i + j – 1) 指定.
n — 矩阵的阶次, 非负整数标量,
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def hilb_type_matrix(input_str):
try:
# 使用 sympy 的 sympify 函数将输入的字符串转换为 sympy 表达式
expr = sp.sympify(input_str)
# 初始化错误标志为 False，表示没有错误
error = False
# 初始化结果变量为 None
result = None
def hilbert_matrix(n):
"""
创建一个 n x n 的希尔伯特矩阵。
希尔伯特矩阵的元素 H[i, j] 定义为 1 / (i + j + 1)，其中 i 和 j 是矩阵的行和列索引（从 0 开始）。
参数:
n (int): 矩阵的阶数
返回:
np.ndarray: 一个 n x n 的希尔伯特矩阵
"""
# 创建一个 n x n 的零矩阵
H = np.zeros((n, n))
# 遍历矩阵的每一行
for i in range(n):
# 遍历矩阵的每一列
for j in range(n):
# 计算希尔伯特矩阵的元素值
H[i, j] = 1 / (i + j + 1)
return H
# 检查解析后的表达式是否为元组
if isinstance(expr, tuple):
# 如果是元组，将错误标志设置为 True
error = True
# 检查解析后的表达式是否为整数
elif expr.is_integer:
# 如果是整数，调用 hilbert_matrix 函数生成希尔伯特矩阵
result = hilbert_matrix(int(expr))
else:
# 如果既不是元组也不是整数，将错误标志设置为 True
error = True
# 如果没有错误，将结果转换为 sympy 矩阵并返回
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
# 如果在处理过程中出现异常，返回错误信息
return f"错误:{e}"
print(hilb_type_matrix("4"))
# Matrix([[1.00000000000000, 0.500000000000000, 0.333333333333333, 0.250000000000000],
# [0.500000000000000, 0.333333333333333, 0.250000000000000, 0.200000000000000],
# [0.333333333333333, 0.250000000000000, 0.200000000000000, 0.166666666666667],
# [0.250000000000000, 0.200000000000000, 0.166666666666667, 0.142857142857143]])
直方图箱计数
[N, edges] = histcounts（X）将X值划分为bin, 并返回bin计数和bin边. histcounts函数使用自动分箱算法,, 该算法返回选择的均匀分箱, 以覆盖X中的元素范围, 并揭示分布的潜在形状.
[N, edges] = histcounts（X, nbins）使用标量nbins指定的多个容器.
[N, edges] = histcounts（X. edges）将X排序到具有向量edges指定的bin边的bin中.
X —— 要在各个存储箱之间分配的数据, 向量, 矩阵
nbins —— 箱数, 正整数, 指定为正整数的容器数. 如果不指定nbins, 则histcounts会根据X中的值自动计算要使用的箱数.
edges —— Bin边, 矢量, 指定为向量. 第一个向量元素指定第一个bin的前缘, 最后一个元素指定最后一个bin的后缘. 后缘仅包含在最后一个bin中.
示例1. 统计学与数据分析
正态分布数据:
normal_data = list(np.random.normal(0, 1, 1000))
result1 = histcounts(normal_data)
print(f"计数: {result1[0][:10]}...") # 只显示前10个
print(f"边界: {result1[1][:5]}...") # 只显示前5个边界
#计数: [1, 10, 36, 85, 186, 241, 226, 147, 49, 19]...
#边界: [-3.6945630626203694, -3.0571194676833393, -2.419675872746309, -1.7822322778092785, -1.1447886828722482]...
指定bin数量
result1b = histcounts(normal_data, 20)
print(f"\n使用20个bins: {len(result1b[0])}个区间")
#使用20个bins: 20个区间
示例2. 金融收益率分布
模拟股票收益率（正态分布，但可能有厚尾）
returns = list(np.random.normal(0.001, 0.02, 500)) # 日收益率 ~ N(0.1%, 2%)
result2 = histcounts(returns, [-0.1, -0.05, -0.02, -0.01, 0, 0.01, 0.02, 0.05, 0.1])
print(f"收益率分布: {result2}")
#收益率分布: ([2, 63, 76, 99, 84, 81, 93, 2], [-0.1, -0.05, -0.02, -0.01, 0.0, 0.01, 0.02, 0.05, 0.1])
示例3. 质量控制与工程
零件尺寸测量
part_sizes = list(np.random.normal(10.0, 0.1, 200)) # 目标尺寸10.0mm，标准差0.1mm
tolerance_limits = [9.7, 9.8, 9.9, 10.0, 10.1, 10.2, 10.3]
result3 = histcounts(part_sizes, tolerance_limits})
print(f"零件尺寸分布: {result3}")
#零件尺寸分布: ([6, 26, 64, 67, 29, 5], [9.7, 9.8, 9.9, 10.0, 10.1, 10.2, 10.3])
示例4. 医学与生物统计学
血压数据分布
blood_pressure = list(np.random.normal(120, 15, 300)) # 收缩压
bp_bins = [80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
result4 = histcounts(blood_pressure, bp_bins)
print(f"血压分布: {result4}")
#血压分布: ([3, 22, 45, 81, 64, 50, 31, 1, 1, 0], [80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 170.0, 180.0])
示例5. 环境监测数据
温度数据（季节性变化）
temperatures = list(np.concatenate([
np.random.normal(15, 5, 100), # 春季
np.random.normal(25, 3, 100), # 夏季
np.random.normal(10, 4, 100), # 秋季
np.random.normal(-5, 8, 100) # 冬季
]))
temp_bins = list(range(-20, 41, 5)) # 从-20°C到40°C，每5°C一个区间
result5 = histcounts(temperatures, temp_bins)
print(f"温度分布: {result5}")
#温度分布: ([6, 14, 23, 28, 30, 54, 62, 60, 50, 67, 2, 0], [-20.0, -15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0])
示例6. 计算机系统性能分析
响应时间分布（通常是指数分布）
response_times = list(np.random.exponential(0.1, 500)) # 平均响应时间0.1秒
rt_bins = [0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0]
result6 = histcounts(response_times, rt_bins)
print(f"响应时间分布: {result6}")
#响应时间分布: ([36, 42, 99, 124, 125, 72, 2, 0], [0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0])
示例7. 社会科学调查数据
年龄分布（多峰分布）
ages = list(np.concatenate([
np.random.normal(25, 5, 150), # 年轻人群
np.random.normal(45, 8, 200), # 中年人群
np.random.normal(65, 10, 100) # 老年人群
]))
age_bins = list(range(0, 101, 10)) # 每10岁一个区间
result7 = histcounts(ages, age_bins)
print(f"年龄分布: {result7}")
#年龄分布: ([1, 21, 118, 70, 101, 67, 44, 19, 6, 3], [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0])
示例8. 学生成绩分析
考试成绩（偏态分布）
scores = list(np.random.beta(2, 5, 300) * 100) # 偏向左的分布
score_bins = [0, 20, 40, 60, 70, 80, 90, 100]
result8 = histcounts(scores, score_bins)
print(f"成绩分布: {result8}")
#成绩分布: ([92, 132, 60, 10, 6, 0, 0], [0.0, 20.0, 40.0, 60.0, 70.0, 80.0, 90.0, 100.0])
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def histogram_counts(input_str):
"""
对标MATLAB的histcounts函数，计算直方图的bin计数和边界。
参数:
input_str (str): 输入字符串，应为可解析为数值数据及bins参数的Python字面量。
返回:
tuple or str: 成功时返回(计数列表, 边界列表)，失败时返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
data = None
bins_arg = None
if isinstance(expr, tuple) and len(expr) == 2:
if isinstance(expr[0], list):
data = np.array(expr[0], dtype=float).ravel()
bins_arg = expr[1]
else:
error = True
elif isinstance(expr, list):
data = np.array(expr, dtype=float).ravel()
else:
error = True
# 处理bins参数
if bins_arg is not None:
if isinstance(bins_arg, (list, tuple, np.ndarray)):
bins = np.array(bins_arg, dtype=float)
elif bins_arg.is_number:
bins = int(bins_arg)
else:
return f"bins参数类型错误: {type(bins_arg)}"
N, edges = np.histogram(data, bins=bins)
else:
N, edges = np.histogram(data)
return list(N), list(edges) if not error else f"输入错误: {input_str}"
except SyntaxError:
return f"语法错误: 无法解析输入 {input_str}"
except ValueError as ve:
return f"数值错误: {str(ve)}"
except Exception as e:
return f"错误: {str(e)}"
# 基本使用
print(histogram_counts("[1, 2, 3, 4]"))
# ([1, 0, 0, 1, 0, 0, 1, 0, 0, 1], [1.0, 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.6999999999999997, 4.0])
# 指定bin数量
print(histogram_counts("([1,2,3,4,5], 2)"))
# ([2, 3], [1.0, 3.0, 5.0])
# 指定bin边界
print(histogram_counts("([1,2,3,4], [0, 2.5, 5])"))
# ([2, 2], [0.0, 2.5, 5.0])
二元直方图bin计数
[N,Xedges,Yedges] = histcounts2(X,Y)将X和Y中的值划分为多个二维bin, 并返回bin计数以及每个维度的bin边界. histcounts2函数使用自动分bin算法, 然后返回均匀bin, 这些bin可涵盖X和Y中值的范围并显示分布的基本形状.
[N,Xedges,Yedges] = histcounts2(X,Y,nbins) 指定要在每个维度中使用的 bin 数量.
[N,Xedges,Yedges] = histcounts2(X,Y,Xedges,Yedges) 按照 Xedges 和 Yedges 指定的 bin 边界, 将 X 和 Y 划分为多个 bin.
X,Y — 要分布到各 bin 的数据（以单独参量指定）,向量, 矩阵
nbins — 每个维度中的 bin 数量, 正整数标量, 由正整数组成的二元素向量
Xedges — x 维度中的 bin 边界, 向量
Yedges — y 维度中的 bin 边界, 向量
N — bin 计数, 数组
Xedges — x 维度中的 bin 边界, 向量
Yedges — y 维度中的 bin 边界, 向量
示例1. 金融数据分析 - 股票收益率与波动率
生成模拟数据：正收益通常伴随低波动，负收益伴随高波动
np.random.seed(42)
returns = np.concatenate([
np.random.normal(0.02, 0.01, 300), # 正收益，低波动
np.random.normal(-0.01, 0.03, 200) # 负收益，高波动
])
volatility = np.concatenate([
np.random.normal(0.15, 0.02, 300), # 低波动
np.random.normal(0.25, 0.05, 200) # 高波动
])
result1 = histcounts2(returns, volatility, 10)
print(f"收益率-波动率联合分布形状: {np.array(result1[0]).shape}")
#收益率-波动率联合分布形状: (10, 10)
#可用于分析风险-收益关系
示例2. 医学研究 - 身高体重指数(BMI)分析
生成身高体重数据（正相关）
heights = np.random.normal(170, 10, 1000) # 身高 cm
weights = 0.8 * heights - 100 + np.random.normal(0, 8, 1000) # 体重 kg
height_bins = list(range(140, 201, 5)) # 140-200cm，每5cm一个区间
weight_bins = list(range(40, 121, 5)) # 40-120kg，每5kg一个区间
result2 = histcounts2(heights, weights, height_bins, weight_bins)
print(f"身高体重联合分布: {len(result2[1]) - 1}x{len(result2[2]) - 1} 个区间")
#身高体重联合分布: 12x16 个区间
示例3. 环境科学 - 温度与湿度的关系
温度湿度通常呈负相关（高温低湿，低温高湿）
temperatures = np.random.uniform(10, 35, 500)
humidity = 80 - 1.5 * temperatures + np.random.normal(0, 10, 500)
humidity = np.clip(humidity, 20, 100) # 限制湿度范围
result3 = histcounts2(temperatures, humidity, 8)
print(f"温度-湿度联合分布: {np.array(result3[0]).shape}")
#温度-湿度联合分布: (8, 8)
示例4. 市场营销 - 客户年龄与消费金额
不同年龄段的消费模式
ages = np.concatenate([
np.random.normal(25, 5, 200), # 年轻人，中等消费
np.random.normal(45, 8, 300), # 中年人，高消费
np.random.normal(65, 10, 150) # 老年人，低消费
])
spending = np.concatenate([
np.random.normal(500, 100, 200), # 年轻人
np.random.normal(800, 150, 300), # 中年人
np.random.normal(300, 80, 150) # 老年人
])
age_bins = [18, 25, 35, 45, 55, 65, 75]
spending_bins = [0, 200, 400, 600, 800, 1000, 1200]
result4 = histcounts2(ages, spending, age_bins, spending_bins)
print(f"客户年龄-消费金额分布: {len(result4[1]) - 1}x{len(result4[2]) - 1} 个区间")
#客户年龄-消费金额分布: 6x6 个区间
示例5. 工业质量控制 - 零件尺寸相关性
两个相关尺寸的测量数据
diameter = np.random.normal(50.0, 0.1, 300)
length = 2.0 * diameter + np.random.normal(0, 0.05, 300) # 长度与直径相关
result5 = histcounts2(diameter, length, 6)
print(f"直径-长度联合分布: {np.array(result5[0]).shape}")
#直径-长度联合分布: (6, 6)
#可用于检测制造过程的一致性
示例6. 教育研究 - 数学与物理成绩的关系
math_scores = np.random.normal(75, 15, 400)
physics_scores = 0.8 * math_scores + np.random.normal(0, 8, 400) # 物理与数学正相关
score_bins = list(range(0, 101, 10))
result6 = histcounts2(math_scores, physics_scores, score_bins, score_bins)
print(f"数学-物理成绩联合分布: {len(result6[1]) - 1}x{len(result6[2]) - 1} 个区间")
#数学-物理成绩联合分布: 10x10 个区间
示例7. 气象学 - 风速与风向的关系
模拟不同风向的风速分布
directions = np.random.uniform(0, 360, 600) # 风向 0-360度
speeds = 5 + 3 * np.abs(np.sin(np.radians(directions))) + np.random.normal(0, 1, 600) # 特定风向风速较高
dir_bins = list(range(0, 361, 45)) # 8个方向
speed_bins = [0, 2, 4, 6, 8, 10, 12, 14, 16]
result7 = histcounts2(directions, speeds, dir_bins, speed_bins)
print(f"风向-风速联合分布: {len(result7[1]) - 1}x{len(result7[2]) - 1} 个区间")
#风向-风速联合分布: 8x8 个区间
示例8. 计算机科学 - CPU使用率与内存使用率
cpu_usage = np.random.beta(2, 5, 500) * 100 # 偏态分布
memory_usage = 0.6 * cpu_usage + np.random.normal(0, 10, 500) # 与CPU相关但更分散
usage_bins = list(range(0, 101, 10))
result8 = histcounts2(cpu_usage, memory_usage, usage_bins, usage_bins)
print(f"CPU-内存使用率联合分布: {len(result8[1]) - 1}x{len(result8[2]) - 1} 个区间")
#CPU-内存使用率联合分布: 10x10 个区间
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def histogram_counts_second(input_str):
"""
对标MATLAB histcounts2的二元直方图计数函数
参数:
input_str (str): SymPy可解析的字符串，格式为：
(1) (X, Y)
(2) (X, Y, bins)
(3) (X, Y, Xedges, Yedges)
返回:
(Matrix(N), Xedges, Yedges) 或 错误信息字符串
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查是否为包含两个矩阵的元组
if isinstance(expr, tuple) and isinstance(expr[0], list) and isinstance(expr[1], list):
# 转换输入数据为numpy数组
x = expr[0]
y = expr[1]
X = np.ravel(np.array(x, dtype=float))
Y = np.ravel(np.array(y, dtype=float))
# 参数解析逻辑
if len(expr) >= 4: # 处理(X, Y, Xedges, Yedges)情况
if isinstance(expr[2], list) and isinstance(expr[3], list):
x_edges = expr[2]
y_edges = expr[3]
bins = [
np.ravel(np.array(x_edges, dtype=float)),
np.ravel(np.array(y_edges, dtype=float))
]
N, Xe, Ye = np.histogram2d(X, Y, bins=bins)
result = (sp.Matrix(N), list(Xe), list(Ye))
else:
error = True
elif len(expr) == 3: # 处理bins参数
bins_arg = expr[2]
# 处理整数型bins
if isinstance(bins_arg, sp.Number) and bins_arg.is_integer:
N, Xe, Ye = np.histogram2d(X, Y, bins=int(bins_arg))
# 处理列表/矩阵型bins
elif isinstance(bins_arg, list):
bins_mat = bins_arg
bins = np.array(bins_mat)
N, Xe, Ye = np.histogram2d(X, Y, bins=bins)
else:
error = True
result = (sp.Matrix(N), list(Xe), list(Ye)) if not error else None
else: # 默认情况
N, Xe, Ye = np.histogram2d(X, Y)
result = (sp.Matrix(N), list(Xe), list(Ye))
else:
error = True
return result if not error else f"输入格式错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 案例1：基本使用
print(histogram_counts_second("([1,2,3], [4,5,6])"))
# (Matrix([
# [1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]]),
#
# [1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4000000000000004, 2.6, 2.8, 3.0],
# [4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0])
# 案例2：指定整数bins
print(histogram_counts_second("[1,1,2,3,2,2,1,1,2,3],[5,6,3,8,9,1,2,7,5,1],8"))
# (Matrix([
# [ 0, 1.0, 0, 0, 1.0, 1.0, 1.0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0],
# [1.0, 0, 1.0, 0, 1.0, 0, 0, 1.0],
# [ 0, 0, 0, 0, 0, 0, 0, 0],
# [ 0, 0, 0, 0, 0, 0, 0, 0],
# [1.0, 0, 0, 0, 0, 0, 0, 1.0]]),
#
# [1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0],
# [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
# 案例3：指定二维bins
print(histogram_counts_second("([1,2,3], [4,5,6], [2,3])"))
# (Matrix([[1.0, 0, 0],
# [ 0, 1.0, 1.0]]),
#
# [1.0, 2.0, 3.0], [4.0, 4.666666666666667, 5.333333333333333, 6.0])
# 案例4：直接指定边界
print(histogram_counts_second("([1,2,3], [4,5,6], [0,2,4], [3,5,7])"))
# (Matrix([[1.0, 0],
# [ 0, 2.0]]),
#
# [0.0, 2.0, 4.0],
# [3.0, 5.0, 7.0])
直方图
直方图是一种条形图, 它将数据分组为bin. 创建Histogram对象后, 可以通过更改直方图的属性值修改它的各个方面. 这对快速修改bin属性或更改显示特别有用
HistogramPlot(X) 基于X创建直方图. histogram函数使用自动分bin算法, 然后返回均匀宽度的bin, 这些bin可涵盖 X 中的元素范围并显示分布的基本形状. histogram将bin显示为矩形条, 这样每个矩形的高度就表示bin中的元素数量.
HistogramPlot(X,nbins)指定bin的数量.
HistogramPlot(X,edges) 将X划分为在向量中指定bin边界的bin.
X — 要分布到各 bin 的数据, 向量, 矩阵
nbins — bin 数量, 正整数
edges — bin 边界, 向量
电商用户消费金额分析: 分析双11期间用户的单笔订单金额分布
HistogramPlot([89,120,150,199,210,299,320,399,450,499,520,599,620,799,850,999,1200,1500],[0,200,400,600,800,1000,2000])
工厂零件尺寸质检: 检测流水线生产的螺栓直径（标准值10.0±0.2mm）
HistogramPlot([9.7,9.8,9.85,9.9,9.95,10.0,10.0,10.05,10.1,10.1,10.15,10.2,10.2,10.25,10.3,10.4],8)
APP用户活跃时长分析: 统计每日用户使用时长（分钟）
HistogramPlot([2,5,8,10,12,15,15,18,20,25,30,35,40,45,50,60,75,90,120],[0,1,5,15,30,60,180])
医院患者年龄分布: 某科室疾病患者的年龄分布
HistogramPlot([32,35,38,40,42,45,48,50,52,55,58,60,62,65,68,70,72,75,80,85],[30,40,50,60,70,90])
网络请求响应时间监控: API接口的响应时间（毫秒）
HistogramPlot([12,15,18,20,22,25,28,30,35,40,45,50,60,75,100,150,200,500],[10,20,50,100,200,1000])
二元直方图
二元直方图是一种数值数据条形图, 它将数据分组到二维bin中. 创建Histogram2Plot对象后, 可以通过更改直方图的属性值修改它的各个方面. 这对快速修改bin属性或更改显示特别有用
histogram2Plot(X,Y) 创建 X 和 Y 的二元直方图. histogram2 函数使用自动分 bin 算法, 然后返回均匀面积的 bin, 这些 bin 可涵盖 X 和 Y 中的元素范围并显示分布的基本形状. histogram2 将 bin 显示为二维矩形条形，这样每个条形的高度就表示 bin 中的元素数量.
histogram2Plot(X,Y,nbins) 指定要在每个维度中使用的 bin 数量.
X,Y — 要分布到各 bin 的数据, 向量, 矩阵
nbins — bin 数量, 正整数
身高体重关系（健康数据分析）:
深色区域显示最常见的身高体重组合（如170cm/65kg），右上角稀疏区可能是异常值
Histogram2Plot([165,170,178,162,185], # 身高(cm)
[55,68,72,60,80], # 体重(kg)
12)
广告点击分析（数字营销）:
观看时长8-15秒且点击率1-2%的区域密度最高，短时长(<5s)点击率波动大
Histogram2Plot([12,8,25,3,17], # 广告观看时长(秒)
[0.2,1.5,0.8,3.2,1.6], # 点击率(%)
10)
气象模式分析（气候科学）:
高温(>30℃)与低降雨(<2mm)强相关，25℃中温区降雨分布最广
Histogram2Plot([28,32,25,30], # 日间温度(℃)
[0,5.2,12.8,3.4], # 降雨量(mm)
15)
股票收益率分析（金融量化）:
高波动率(>0.2)时收益率分布更分散，低波动率时收益率集中在±5%
Histogram2Plot([0.15,0.22,0.08], # 波动率
[-0.03,0.12,0.05], # 日收益率
20)
城市通勤研究（交通规划）:
10km/25min区域密度最高，20km以上出现明显的交通效率下降
Histogram2Plot([8,15,5,20], # 通勤距离(km)
[30,45,20,60], # 通勤时间(分钟)
12)
霍纳法则(秦九韶算法)
把表达式转成嵌套形式，提高计算效率.
horner(p)返回多项式p的霍纳形式,
horner(p，var)使用var中的变量.
p是符号表达式,符号函数,符号表达式数组,符号函数数组.
var是变量,符号变量,符号变量数组.
示例1：数值计算优化 - 减少乘法次数
print(horner(3*x^4 + 5*x^3 - 2*x^2 + 7*x + 1))
# 输出: x*(x*(x*(3*x + 5) - 2) + 7) + 1
# 原始需要10次乘法，霍纳形式只需4次乘法
示例2：多项式插值
print(horner(2*x^3 - 3*x^2 + 4*x - 5))
# 输出: x*(x*(2*x - 3) + 4) - 5
示例3：控制系统中的传递函数
print(horner(s^3 + 2*s^2 + 3*s + 4))
# 输出: s*(s*(s + 2) + 3) + 4
示例4：矩阵多项式 - 在数值线性代数中常见
matrix_poly = [[x^3 + 2*x^2 + x + 1, x^2 + 3*x + 2], [4*x^3 - x^2 + 2*x, x^4 + x + 1]]
result = horner(matrix_poly)
print(result)
# 输出: [[x*(x*(x + 2) + 1) + 1, x*(x + 3) + 2],
[x*(x*(4*x - 1) + 2), x*(x^3 + 1) + 1]]
示例5：多元多项式 - 指定变量
print(horner((x*y^3 + 2*x*y^2 + 3*x*y + 4*x, y)))
# 输出: 4*x + y*(3*x + y*(x*y + 2*x))
示例6：多元多项式（指定变量x）
print(horner((x^3*y + 2*x^2*y + 3*x*y + 4*y, x)))
# 输出: x*(x*(x*y + 2*y) + 3*y) + 4*y
示例7：高次多项式 - 计算效率提升明显
print(horner(x^5 + 4*x^4 + 6*x^3 + 4*x^2 + x + 1))
# 输出: x*(x*(x*(x*(x + 4) + 6) + 4) + 1) + 1
示例8：特征多项式
print(horner(x^3 - 6*x^2 + 11*x - 6))
# 输出: x*(x*(x - 6) + 11) - 6
示例9：数值稳定性测试（接近零的系数）
print(horner(x^4 + 0.001*x^3 + 0.0001*x^2 + x + 1))
# 输出: x*(x*(x*(1.0*x + 0.001) + 0.0001) + 1.0) + 1.0
示例10：复杂矩阵表达式
complex_matrix = [[x^2 + y*x + y^2, x^3 + 2*x*y], [x + y^2, x^2*y + x + 1]]
print(horner(complex_matrix, y))
#输出: [[x^2 + y*(x + y), x^3 + 2*x*y],
[x + y^2, x^2*y + x + 1]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.matrices import MatrixBase, MutableDenseMatrix
from sympy.polys.polyfuncs import horner
def horner_nested_polynomial(input_str):
"""
将输入的多项式或矩阵中的每个元素转换为霍纳嵌套形式。
参数:
input_str: 字符串，表示SymPy可解析的表达式。可以是多项式、矩阵或元组（矩阵或表达式，变量）。
返回:
SymPy表达式或矩阵，表示霍纳形式的多项式。若输入错误则返回错误信息。
示例:
>>> x, y = sp.symbols('x y')
>>> horner_nested_polynomial("x**3 + 2*x**2 + 3*x +4")
x*(x*(x + 2) + 3) + 4
>>> horner_nested_polynomial("[[x**3, x**2 + 2*x], [x+1, 5]]")
Matrix([
[ x**3, x*(x + 2)],
[ x + 1, 5]])
>>> horner_nested_polynomial("(x**2 + x*y + y**2, y)")
y*(x + y) + x**2
>>> horner_nested_polynomial("([[x**3 + y*x, x + y]], y)")
Matrix([[x*y + x**3, x + y]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def eval_horner(expression, wrt_var=None):
"""内部函数：应用霍纳法则到单个表达式"""
if wrt_var is not None:
poly = sp.poly(expression, wrt_var)
else:
poly = sp.poly(expression)
return horner(poly, wrt=wrt_var)
# 处理元组输入（表达式/矩阵，变量）
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.free_symbols for e in expr):
result = eval_horner(*expr)
else:
error = True
elif expr.free_symbols:
# 处理普通表达式
result = eval_horner(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
print(horner_nested_polynomial("9*x**4+8*x**3+7*x**2+6*x+5"))
# x*(x*(x*(9*x + 8) + 7) + 6) + 5
希尔伯特变换
H = htrans(f) 返回符号函数 f 的希尔伯特变换. 默认情况下, 独立变量为 t, 变换变量为 x.
f — 输入, 符号表达式, 符号函数
示例1: 正弦信号 - 基本测试
result = htrans(sin(t))
t_vals, H, y = result
plt.figure(figsize=(12, 8))
plt.subplot(2, 3, 1)
plt.plot(t_vals, y, 'b-', label='Original: sin(t)')
plt.plot(t_vals, H, 'r--', label='Hilbert: -cos(t)')
plt.plot(t_vals, -np.cos(t_vals), 'g:', label='Theoretical: -cos(t)')
plt.legend()
plt.title('正弦信号')
plt.grid(True)
示例2: 余弦信号
result = htrans(cos(t))
t_vals, H, y = result
plt.subplot(2, 3, 2)
plt.plot(t_vals, y, 'b-', label='Original: cos(t)')
plt.plot(t_vals, H, 'r--', label='Hilbert: sin(t)')
plt.plot(t_vals, np.sin(t_vals), 'g:', label='Theoretical: sin(t)')
plt.legend()
plt.title('余弦信号')
plt.grid(True)
示例3: 调幅信号 (AM信号)
result = htrans((1+0.5*cos(2*t))*cos(10*t))
t_vals, H, y = result
plt.subplot(2, 3, 3)
plt.plot(t_vals, y, 'b-', label='AM Signal')
plt.plot(t_vals, H, 'r--', label='Hilbert Transform')
plt.legend()
plt.title('调幅信号')
plt.grid(True)
# 计算包络
envelope = np.abs(signal.hilbert(y))
plt.plot(t_vals, envelope, 'g-', linewidth=2, label='Envelope')
plt.plot(t_vals, -envelope, 'g-', linewidth=2)
示例4: 调频信号 (FM信号)
result = htrans(cos(10*t + 2*sin(2*t)))
t_vals, H, y = result
plt.subplot(2, 3, 4)
plt.plot(t_vals, y, 'b-', label='FM Signal')
plt.plot(t_vals, H, 'r--', label='Hilbert Transform')
plt.legend()
plt.title('调频信号')
plt.grid(True)
示例5: 脉冲信号 (矩形脉冲)
创建一个矩形脉冲
t_vals = np.linspace(0, 2 * np.pi, 1000)
y_rect = np.where((t_vals > np.pi / 2) & (t_vals < 3 * np.pi / 2), 1, 0)
H_rect = np.imag(signal.hilbert(y_rect))
plt.subplot(2, 3, 5)
plt.plot(t_vals, y_rect, 'b-', label='Rectangular Pulse')
plt.plot(t_vals, H_rect, 'r--', label='Hilbert Transform')
plt.legend()
plt.title('矩形脉冲')
plt.grid(True)
示例6: 多频信号
result = htrans(sin(t) + 0.5*sin(3*t))
t_vals, H, y = result
plt.subplot(2, 3, 6)
plt.plot(t_vals, y, 'b-', label='Multi-frequency')
plt.plot(t_vals, H, 'r--', label='Hilbert Transform')
plt.legend()
plt.title('多频信号')
plt.grid(True)
plt.tight_layout()
plt.show()
示例7: 解析信号和瞬时频率计算
result = htrans(cos(5*t + sin(2*t))) # FM信号
t_vals, H, y = result
analytic_signal = signal.hilbert(y)
# 计算瞬时幅度（包络）
instantaneous_amplitude = np.abs(analytic_signal)
# 计算瞬时相位
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
# 计算瞬时频率
instantaneous_frequency = (np.diff(instantaneous_phase) /
(2.0 * np.pi * (t_vals[1] - t_vals[0])))
plt.figure(figsize=(15, 10))
plt.subplot(3, 1, 1)
plt.plot(t_vals, y, 'b-', label='Original Signal')
plt.plot(t_vals, instantaneous_amplitude, 'r--', label='Instantaneous Amplitude')
plt.legend()
plt.title('信号和瞬时包络')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_vals, instantaneous_phase, 'g-')
plt.title('瞬时相位')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t_vals[:-1], instantaneous_frequency, 'm-')
plt.title('瞬时频率')
plt.xlabel('Time')
plt.grid(True)
plt.tight_layout()
plt.show()
示例8: 单边带调制 (SSB)
基带信号
baseband_freq = 2
carrier_freq = 20
t_vals = np.linspace(0, 2 * np.pi, 1000)
baseband_signal = np.cos(baseband_freq * t_vals)
carrier = np.cos(carrier_freq * t_vals)
生成解析信号（单边带）
analytic_baseband = signal.hilbert(baseband_signal)
ssb_signal = np.real(analytic_baseband * np.exp(1j * carrier_freq * t_vals))
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(t_vals, baseband_signal)
plt.title('基带信号')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t_vals, ssb_signal)
plt.title('单边带调制信号')
plt.grid(True)
# 频谱分析
from scipy.fft import fft, fftfreq
plt.subplot(3, 1, 3)
fft_ssb = fft(ssb_signal)
freqs = fftfreq(len(ssb_signal), t_vals[1] - t_vals[0])
plt.plot(freqs[:len(freqs) // 2], np.abs(fft_ssb[:len(fft_ssb) // 2]))
plt.title('单边带信号频谱')
plt.xlabel('Frequency (Hz)')
plt.grid(True)
plt.tight_layout()
plt.show()
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
def hilbert_transform_func(input_str):
"""
对标MATLAB的htrans函数，计算输入表达式的数值希尔伯特变换。
参数:
input_str (str): 以't'为变量的数学表达式，例如'sin(t)'、'exp(t)*cos(t)'等。
返回:
tuple: (t_values, H_values) 包含时间序列和对应的希尔伯特变换数值结果。
或 str: 错误信息。
"""
try:
# 步骤1: 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
t = sp.symbols('t')
# 步骤2: 检查表达式是否含有符号't'，若无则报错
if not expr.has(t):
return f"错误: 表达式必须包含变量't'，但输入为'{expr}'"
# 步骤3: 替换所有其他符号为1（例如，处理如'a*sin(t)'中的'a'）
free_symbols = expr.free_symbols - {t}
subs_dict = {symbol: 1 for symbol in free_symbols}
expr_substituted = expr.subs(subs_dict)
# 步骤4: 将SymPy表达式转换为NumPy函数
expr_numpy = sp.lambdify(t, expr_substituted, modules='numpy')
# 步骤5: 生成时间序列（0到2π，1000个点）
t_vals = np.linspace(0, 2 * np.pi, 1000)
# 步骤6: 计算原始信号y(t)
y = expr_numpy(t_vals)
# 步骤7: 使用SciPy计算希尔伯特变换（注意：返回的是解析信号，虚部为Hilbert变换）
analytic_signal = signal.hilbert(y)
H = np.imag(analytic_signal)
return (t_vals, H, y) # 同时返回原始信号用于比较
except sp.SympifyError:
return f"错误: 无法解析输入表达式 '{input_str}'"
except Exception as e:
return f"错误: {str(e)}"
# 示例使用
if __name__ == "__main__":
# 示例1: 计算sin(t)的希尔伯特变换（理论值为-cos(t)）
input_expr = "sin(t)"
result = hilbert_transform_func(input_expr)
if isinstance(result, tuple):
import matplotlib.pyplot as plt
t_vals, H, y = result
plt.figure(figsize=(10, 6))
plt.plot(t_vals, H, '--', label='Computed Hilbert')
plt.xlabel('t')
plt.legend()
plt.title('Hilbert Transform Demo')
plt.show()
else:
print(result)
赫尔维茨黎曼函数
Z=Hurwitz-zeta(s,a)为数字或符号输入s和a计算Hurwitz zeta函数.只有当s不为1且a既不是0也不是负整数时,才定义Hurwitz zeta函数.
Z=hurwitzZeta(n,s,a)返回参照变量s的n阶hurwitz-zeta(s,a)的导数.
s是数字,数组,符号数,符号变量,符号函数,符号表达式,符号数组.
a是数字,数组,符号数,符号变量,符号函数,符号表达式,符号数组.
n是导数的阶数,非负整数.
示例1: 黎曼ζ函数（特殊情况）
print("ζ(2) =", hurwitzZeta(2, 1)
print("ζ(4) =", hurwitzZeta(4, 1)
print("ζ(0) =", hurwitzZeta(0, 1)
#ζ(2) = 1.6449340668482264
#ζ(4) = 1.0823232337111381
#ζ(0) = -0.5
示例2: 狄利克雷级数
print("ζ(2, 0.5) =", hurwitzZeta(2, 0.5)
print("ζ(3, 0.25) =", hurwitzZeta(3, 0.25)
#ζ(2, 0.5) = 4.934802200544679
#ζ(3, 0.25) = 64.66386996876847
示例3: 导数计算
print("ζ''(2, 1) =", hurwitzZeta(2, 1, 2)
#ζ''(2, 1) = oo
示例4: 数值数组计算
s_values = [2, 3, 4]
a_values = [0.5, 1.0, 1.5]
for s in s_values:
for a in a_values:
result = hurwitzZeta(s, a)
print(f"ζ({s}, {a}) = {result}")
#ζ(2, 0.5) = 4.934802200544679
#ζ(2, 1.0) = 1.6449340668482264
#ζ(2, 1.5) = 0.9348022005446793
#ζ(3, 0.5) = 8.41439832211716
#ζ(3, 1.0) = 1.2020569031595942
#ζ(3, 1.5) = 0.41439832211716
#ζ(4, 0.5) = 16.234848505667074
#ζ(4, 1.0) = 1.0823232337111381
#ζ(4, 1.5) = 0.23484850566707288
示例5: 在数论中的应用 - 狄利克特征
L(s, χ) = ∑ χ(n)/n^s 可以用Hurwitz ζ函数表示
def dirichlet_L(s, chi, modulus):
"""使用Hurwitz ζ函数计算狄利克雷L函数"""
result = 0
for a in range(1, modulus + 1):
if chi(a) != 0:
result += chi(a) * hurwitzZeta(s, a / modulus)
return result / (modulus ** s)
# 定义模4的特征
def chi_mod4(n):
if n % 4 == 1:
return 1
elif n % 4 == 3:
return -1
else:
return 0
L_value = dirichlet_L(1, chi_mod4, 4)
print(f"L(1, χ) = {L_value}")
#L(1, χ) = nan
示例6: 在物理学中的应用- 玻色-爱因斯坦统计:
在理想玻色气体中，粒子数密度与ζ(3/2)相关
zeta_3_2 = hurwitzZeta(1.5, 1)
print(f"ζ(3/2) = {zeta_3_2}")
#ζ(3/2) = 2.612375348685488
#这个值出现在玻色-爱因斯坦凝聚的临界温度公式中
示例7: 在概率论中的应用 - Zipf分布
Zipf分布的归一化常数与ζ函数相关
def zipf_constant(s, N):
"""计算Zipf分布的归一化常数"""
return 1.0 / abs(hurwitzZeta(s, 1)
s_zipf = 2.0
norm_const = zipf_constant(s_zipf, 100)
print(f"Zipf分布(s={s_zipf})的归一化常数: {norm_const}")
#Zipf分布(s=2.0)的归一化常数: 0.6079271018540267
示例8: 解析延拓
计算s<1时的ζ函数值（通过解析延拓）
negative_s = [-1, -2, -3]
for s in negative_s:
result = hurwitzZeta(s, 1)
print(f"ζ({s}) = {result}")
# 这些值与伯努利数相关
#ζ(-1) = -0.08333333333333333
#ζ(-2) = 0
#ζ(-3) = 0.008333333333333333
示例9: 特殊函数关系
ζ(s, a) 与多伽玛函数的关系
from mpmath import polygamma
s_poly = 3
a_poly = 0.5
zeta_val = hurwitzZeta(s_poly, a_poly)
polygamma_val = (-1) ** s_poly * polygamma(s_poly - 1, a_poly) / factorial(s_poly - 1)
print(f"ζ({s_poly}, {a_poly}) = {zeta_val}")
print(f"通过多伽玛函数: {polygamma_val}")
#ζ(3, 0.5) = 8.41439832211716
#通过多伽玛函数: 8.41439832211716
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from mpmath import zeta as mp_zeta
def hurwitz_zeta_function(input_expr):
"""
对标MATLAB的hurwitzZeta函数，计算Hurwitz zeta函数或其导数。
Hurwitz zeta函数定义: ζ(s, a) = ∑_{k=0}^∞ 1/(k + a)^s
导数定义: ζ^{(n)}(s, a) = d^n/ds^n ζ(s, a)
参数:
input_expr (str/tuple): 输入表达式，可以是以下形式：
- "(s, a)" 字符串，计算ζ(s, a)
- "(s, a, n)" 字符串，计算n阶导数ζ^{(n)}(s, a)
- 包含s和a的元组或矩阵
n (int, optional): 导数阶数，如果输入元组中已包含则优先使用元组中的n
返回:
SymPy表达式/矩阵或错误信息。
示例:
>>> hurwitz_zeta_function("(2, 0.5)")
-1.46035450880959
>>> hurwitz_zeta_function("(s, a, 1)", s=2, a=0.5)
-2.40411380631918
"""
try:
expr = sp.sympify(input_expr)
error = False
result = None
def hurwitz_zeta_sym(s_val, a_val, n_val=None):
zeta_expr = sp.zeta(s_val, a_val)
if n_val is not None and s_val.free_symbols:
vars = s_val.free_symbols
result = zeta_expr.rewrite(sp.lerchphi).diff(*vars, n_val)
else:
result = zeta_expr.rewrite(sp.lerchphi)
return result
def hurwitz_zeta_math(s_val, a_val, n_val=None):
if n_val:
n = float(n_val)
result = mp_zeta(s_val, a_val, n)
else:
result = mp_zeta(s_val, a_val)
return result
if isinstance(expr, tuple):
if len(expr) == 2:
derivative_order = None
elif len(expr) == 3 and expr[2].is_integer and expr[2] > 0:
derivative_order = expr[2]
else:
error = True
if all(e.is_number for e in expr[:2]):
params = tuple(complex(e) for e in expr[:2])
result = hurwitz_zeta_math(*params, derivative_order)
elif any(e.free_symbols for e in expr[:2]):
result = hurwitz_zeta_sym(*expr)
else:
error = True
return result if not error else f"输入错误: {input_expr}"
except Exception as e:
return f"错误:{e}"
# 示例使用
if __name__ == "__main__":
# 示例1: 标量计算
print("ζ(2, 0.5) =", hurwitz_zeta_function("(2,0.5)"))
# 4.93480220054468
# 示例2: 导数计算
s = sp.Symbol('s')
expr = hurwitz_zeta_function("s, 0.5, 1")
print("ζ'(s, 0.5) =", expr)
# Derivative(lerchphi(1, s, 0.5), s)
超几何函数的累积分布函数
hygecdf（x，M，K，N）使用相应的总体大小M,总体中具有所需特征的项目数量K和提取的样本数量N, 计算x中每个值的超几何cdf. x,M,K和N的向量或矩阵输入必须具有相同的大小. 标量输入被扩展为与其他输入具有相同维数的常数矩阵.
x -- 标量, 符号变量, 向量, 矩阵, 多维数组
M, K, N -- 数值标量, 向量, 矩阵, 多维数组
示例1: 基本的标量计算
print("hygecdf(2, 50, 10, 5) =", hygecdf(2, 50, 10, 5)
#hygecdf(2, 50, 10, 5) = 0.9517396968037909
#解释: 从50个物品中(其中10个是目标物品)，抽取5个，最多抽到2个目标物品的概率
示例2: 质量检测应用
场景: 一批100个产品中有5个次品，随机抽取10个检测
M = 100 # 总体大小
K = 5 # 次品数
n = 10 # 抽样数
质量检测场景: 100个产品中5个次品，抽检10个
for x in range(0, 6):
prob = hygecdf(x, M, K, n)
print(f" 最多发现 {x} 个次品的概率: {float(prob):.4f}")
# 最多发现 0 个次品的概率: 0.5838
# 最多发现 1 个次品的概率: 0.9231
# 最多发现 2 个次品的概率: 0.9934
# 最多发现 3 个次品的概率: 0.9997
# 最多发现 4 个次品的概率: 1.0000
# 最多发现 5 个次品的概率: 1.0000
示例3: 生物统计 - 基因研究
场景: 种群中有1000个体，其中50个有特定基因，随机捕捉100个研究
M_bio = 1000
K_bio = 50
n_bio = 100
基因研究: 1000个体中50个有特定基因，随机捕捉100个
thresholds = [0, 1, 2, 5, 10]
for x in thresholds:
prob = hygecdf(x, M_bio, K_bio, n_bio)
print(f" 最多发现 {x} 个有特定基因个体的概率: {float(prob):.6f}")
# 最多发现 0 个有特定基因个体的概率: 0.004476
# 最多发现 1 个有特定基因个体的概率: 0.030773
# 最多发现 2 个有特定基因个体的概率: 0.105637
# 最多发现 5 个有特定基因个体的概率: 0.616636
# 最多发现 10 个有特定基因个体的概率: 0.992314
示例4: 社会调查 - 投票预测
场景: 城市有10000选民，其中4000支持某候选人，随机调查500人
M_survey = 10000
K_survey = 4000
n_survey = 500
投票调查: 10000选民中4000支持候选人，调查500人
support_levels = [180, 200, 220, 240, 260]
for x in support_levels:
prob = hygecdf(x, M_survey, K_survey, n_survey)
print(f" 最多 {x} 人支持的概率: {float(prob):.4f}")
# 最多 180 人支持的概率: 0.0334
# 最多 200 人支持的概率: 0.5198
# 最多 220 人支持的概率: 0.9721
# 最多 240 人支持的概率: 0.9999
# 最多 260 人支持的概率: 1.0000
示例5: 数组输入计算
x_values = [0, 1, 2, 3, 4]
result_array = hygecdf(x_values, 20, 5, 10)
print("x值:", x_values)
print("CDF值:", [f"{p:.4f}" for p in result_array])
#x值: [0, 1, 2, 3, 4]
#CDF值: ['0.0163', '0.1517', '0.5000', '0.8483', '0.9837']
示例6: 可视化 - CDF曲线
plt.figure(figsize=(15, 10))
# 子图1: 不同总体大小下的CDF
plt.subplot(2, 2, 1)
M_values = [20, 50, 100]
K_ratio = 0.2 # 成功比例
n_ratio = 0.3 # 抽样比例
for M in M_values:
K = int(M * K_ratio)
n = int(M * n_ratio)
x_range = range(0, min(n, K) + 1)
cdf_values = []
for x in x_range:
prob = hygecdf(x, M, K, n)
cdf_values.append(float(prob))
plt.plot(x_range, cdf_values, 'o-', label=f'M={M}, K={K}, n={n}')
plt.xlabel('x (成功次数)')
plt.ylabel('CDF')
plt.title('不同总体大小下的超几何CDF')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图2: 不同抽样比例下的CDF
plt.subplot(2, 2, 2)
M_fixed = 100
K_fixed = 20
n_ratios = [0.1, 0.3, 0.5, 0.7]
for ratio in n_ratios:
n = int(M_fixed * ratio)
x_range = range(0, min(n, K_fixed) + 1)
cdf_values = []
for x in x_range:
prob = hygecdf(x, M_fixed, K_fixed, n)
cdf_values.append(float(prob))
plt.plot(x_range, cdf_values, 'o-', label=f'n={n} (ratio={ratio})')
plt.xlabel('x (成功次数)')
plt.ylabel('CDF')
plt.title('不同抽样比例下的超几何CDF')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图3: 质量控制决策图
plt.subplot(2, 2, 3)
# 场景: 决定可接受的质量水平(AQL)
M_quality = 1000
n_test = 50
defect_rates = [0.01, 0.02, 0.05, 0.1] # 次品率
acceptance_limits = range(0, 11)
for rate in defect_rates:
K_defects = int(M_quality * rate)
prob_accept = []
for limit in acceptance_limits:
prob = hygecdf(limit, M_quality, K_defects, n_test)
prob_accept.append(float(prob))
plt.plot(acceptance_limits, prob_accept, 'o-',
label=f'次品率={rate * 100}%')
plt.xlabel('接受上限 (最大允许次品数)')
plt.ylabel('接受概率')
plt.title('质量控制: OC曲线')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图4: 与二项分布的比较（当总体很大时）
plt.subplot(2, 2, 4)
from scipy.stats import binom
M_large = 10000
K_large = 2000
n_large = 100
p_binom = K_large / M_large # 二项分布的成功概率
x_range = range(0, 31)
hyper_cdf = []
binom_cdf = []
for x in x_range:
# 超几何分布
prob_hyper = hygecdf(x, M_large, K_large, n_large)
hyper_cdf.append(float(prob_hyper))
# 二项分布近似
prob_binom = binom.cdf(x, n_large, p_binom)
binom_cdf.append(prob_binom)
plt.plot(x_range, hyper_cdf, 'bo-', label='超几何分布', markersize=4)
plt.plot(x_range, binom_cdf, 'r--', label='二项分布近似', linewidth=2)
plt.xlabel('x')
plt.ylabel('CDF')
plt.title('超几何分布 vs 二项分布近似 (大总体)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
示例7: 假设检验应用
Fisher精确检验的场景
# 2x2列联表: 治疗组和对照组的成功/失败
# 治疗组: 成功8人，总20人
# 对照组: 成功2人，总20人
M_fisher = 40 # 总人数
K_fisher = 10 # 总成功人数 (8+2)
n_fisher = 20 # 治疗组人数
计算p值 (单边检验)
x_observed = 8
p_value = 1 - float(hygecdf(x_observed - 1, M_fisher, K_fisher, n_fisher)
print(f"Fisher精确检验p值: {p_value:.6f}")
#Fisher精确检验p值: 0.032417
示例8: 符号计算
x_sym, M_sym, K_sym, n_sym = sp.symbols('x M K n')
symbolic_result = hygecdf(x_sym, M_sym, K_sym, n_sym)
符号表达式:
print(symbolic_result)
#Lambda(_k, Sum(Piecewise((binomial(K, _ki)*binomial(-K + M, -_ki + n)/binomial(M, n), (_ki >= 0) & (_ki <= n)), (0, True)), (_ki, 0, _k)))
示例9: 风险评估
金融风险: 投资组合中有100个项目，其中10个高风险
M_portfolio = 100
K_risky = 10
n_selected = 15
投资组合风险评估:
risk_levels = [0, 1, 2, 3, 4, 5]
for x in risk_levels:
prob = hygecdf(x, M_portfolio, K_risky, n_selected)
risk = 1 - float(prob) # 超过x个高风险项目的概率
print(f" 选择超过 {x} 个高风险项目的概率: {risk:.4f}")
# 选择超过 0 个高风险项目的概率: 0.8192
# 选择超过 1 个高风险项目的概率: 0.4625
# 选择超过 2 个高风险项目的概率: 0.1705
# 选择超过 3 个高风险项目的概率: 0.0408
# 选择超过 4 个高风险项目的概率: 0.0063
# 选择超过 5 个高风险项目的概率: 0.0006
示例10: 抽样设计
确定合适的样本量来检测一定比例的问题
M_population = 500
K_problem = 25 # 5%的问题率
confidence = 0.95
抽样设计: 需要多大的样本才能以95%置信度检测到至少一个问题?
sample_sizes = range(10, 101, 10)
for n_sample in sample_sizes:
# 检测到至少一个问题的概率 = 1 - 检测到0个问题的概率
prob_zero = hygecdf(0, M_population, K_problem, n_sample)
prob_at_least_one = 1 - float(prob_zero)
print(f" 样本量 {n_sample}: 检测到至少一个问题的概率 = {prob_at_least_one:.4f}")
if prob_at_least_one >= confidence:
print(f" → 样本量 {n_sample} 满足95%置信度要求")
break
# 样本量 10: 检测到至少一个问题的概率 = 0.4041
# 样本量 20: 检测到至少一个问题的概率 = 0.6488
# 样本量 30: 检测到至少一个问题的概率 = 0.7954
# 样本量 40: 检测到至少一个问题的概率 = 0.8822
# 样本量 50: 检测到至少一个问题的概率 = 0.9330
# 样本量 60: 检测到至少一个问题的概率 = 0.9624
# → 样本量 60 满足95%置信度要求
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Hypergeometric, cdf
from scipy.stats import hypergeom
def hypergeometric_cumulative_distribution(input_str):
try:
# 将输入的字符串转换为SymPy表达式
expr = sp.sympify(input_str)
# 初始化错误标志为False
error = False
# 初始化结果为None
result = None
# 检查表达式是否为长度为4的元组
if isinstance(expr, tuple) and len(expr) == 4:
if all(e.is_integer for e in expr):
params = tuple(int(e) for e in expr)
"""
计算超几何分布的累积分布函数值
参数:
x_val (int): 累积概率的上限值
m_val (int): 总体大小
k_val (int): 总体中成功的次数
n_val (int): 样本大小
返回:
float: 累积分布函数值，即P(X ≤ x_val)
"""
# 直接调用scipy.stats.hypergeom.cdf函数计算累积分布函数值
result = hypergeom.cdf(*params)
elif any(e.free_symbols for e in expr):
"""
计算超几何分布的累积分布函数值。
参数:
x_val: 随机变量的值
m_val: 总体中的成功元素个数
k_val: 总体的元素个数
n_val: 抽样的元素个数
返回:
超几何分布的累积分布函数值或符号表达式
"""
# 创建一个超几何分布随机变量X
X = Hypergeometric(*expr)
# 返回累积分布函数的符号表达式
result = cdf(X)
else:
error = True
else:
# 如果输入不是长度为4的元组，设置错误标志为True
error = True
# 如果没有错误，返回计算结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 如果在计算过程中出现异常，返回错误信息
return f"错误: {e}"
# 标量输入示例
print(hypergeometric_cumulative_distribution("(2, 50, 10, 5)"))
# 0.9517396968037909
# 符号输入示例
print(hypergeometric_cumulative_distribution("(x, a, b, c)"))
# Lambda(_k, Sum(Piecewise((binomial(b, _ki)*binomial(a - b, -_ki + c)/binomial(a, c), (_ki >= 0) & (_ki <= c)), (0, True)), (_ki, 0, _k)))
超几何函数的逆累积分布函数
Y = hygeinv(P,M,K,N)返回最小的整数X,使得在X处计算的超几何cdf等于或超过P.您可以将P视为在N个图纸中观察到X个有缺陷的项目而不替换K有缺陷的M个项目组的概率.
import matplotlib.pyplot as plt
示例1: 基本的标量计算
print("hygeinv(0.5, 100, 20, 10) =", hygeinv(0.5, 100, 20, 10)
#hygeinv(0.5, 100, 20, 10) = 2.0
#解释: 从100个物品中(20个目标)，抽10个，50%分位数是2，即P(X≤2)≥0.5
示例2: 质量控制 - 确定接受标准
场景: 一批1000个产品，允许有2%的次品率(20个次品)，抽检50个
M = 1000
K = 20 # 允许的最大次品数
n = 50 # 抽检数量
confidence_levels = [0.90, 0.95, 0.99]
质量控制: 确定不同置信水平下的接受标准
for conf in confidence_levels:
threshold = hygeinv(conf, M, K, n)
print(f" {conf * 100}%置信水平: 最多接受 {int(threshold)} 个次品")
#90.0%置信水平: 最多接受 2 个次品
#95.0%置信水平: 最多接受 3 个次品
#99.0%置信水平: 最多接受 4 个次品
示例3: 抽样检验设计
确定在不同次品率下的接受阈值
M_design = 500
n_sample = 30
defect_rates = [0.01, 0.02, 0.05, 0.10] # 1%, 2%, 5%, 10%次品率
抽样设计: 不同次品率下的95%置信接受标准
for rate in defect_rates:
K_defects = int(M_design * rate)
threshold = hygeinv(0.95, M_design, K_defects, n_sample)
print(f" 次品率{rate * 100}%: 接受阈值 = {int(threshold)}")
# 次品率1.0%: 接受阈值 = 1
# 次品率2.0%: 接受阈值 = 2
# 次品率5.0%: 接受阈值 = 4
# 次品率10.0%: 接受阈值 = 6
示例4: 数组输入计算
p_values = [0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99]
results = hygeinv(p_values, 100, 25, 20)
print("概率分位数:", p_values)
print("对应阈值: ", results)
#概率分位数: [0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99]
#对应阈值: [3. 4. 5. 6. 7. 8. 9.]
示例5: 风险管理应用
金融投资: 100个项目中10个高风险，选择15个项目组合
M_risk = 100
K_risky = 10
n_portfolio = 15
risk_tolerances = [0.05, 0.10, 0.25, 0.50]
风险管理: 不同风险容忍度下的高风险项目上限
for tolerance in risk_tolerances:
threshold = hygeinv(1 - tolerance, M_risk, K_risky, n_portfolio)
print(f" {tolerance * 100}%风险容忍度: 最多 {int(threshold)} 个高风险项目")
# 5.0%风险容忍度: 最多 3 个高风险项目
# 10.0%风险容忍度: 最多 3 个高风险项目
# 25.0%风险容忍度: 最多 2 个高风险项目
# 50.0%风险容忍度: 最多 1 个高风险项目
示例6: 可视化 - 分位数函数曲线
plt.figure(figsize=(15, 10))
# 子图1: 不同置信水平下的接受阈值
plt.subplot(2, 2, 1)
M_viz = 1000
n_viz = 100
defect_rates_viz = [0.01, 0.02, 0.05, 0.10]
confidence_levels_viz = np.linspace(0.5, 0.99, 50)
for rate in defect_rates_viz:
K_viz = int(M_viz * rate)
thresholds = []
for conf in confidence_levels_viz:
threshold = hygeinv(conf, M_viz, K_viz, n_viz)
thresholds.append(threshold)
plt.plot(confidence_levels_viz, thresholds, label=f'次品率{rate * 100}%')
plt.xlabel('置信水平')
plt.ylabel('接受阈值')
plt.title('不同次品率下的接受阈值')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图2: 不同样本量对分位数的影响
plt.subplot(2, 2, 2)
M_fixed = 500
K_fixed = 25
p_fixed = 0.95
sample_sizes = range(10, 101, 5)
thresholds_sample = []
for n_size in sample_sizes:
threshold = hygeinv(p_fixed, M_fixed, K_fixed, n_size)
thresholds_sample.append(threshold)
plt.plot(sample_sizes, thresholds_sample, 'bo-', markersize=4)
plt.xlabel('样本量')
plt.ylabel('95%分位数')
plt.title('样本量对分位数的影响')
plt.grid(True, alpha=0.3)
# 子图3: 运营决策曲线
plt.subplot(2, 2, 3)
# 库存管理: 确定最大可接受的缺陷数量
M_inventory = 200
n_inspect = 20
confidence_ops = np.linspace(0.8, 0.99, 50)
defect_scenarios = [2, 5, 10] # 不同缺陷数量场景
for defects in defect_scenarios:
thresholds_ops = []
for conf in confidence_ops:
threshold = hygeinv(conf, M_inventory, defects, n_inspect)
thresholds_ops.append(threshold)
plt.plot(confidence_ops, thresholds_ops, label=f'{defects}个缺陷')
plt.xlabel('置信水平')
plt.ylabel('最大观察缺陷数')
plt.title('运营决策: 缺陷检测阈值')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图4: 比较不同总体大小的影响
plt.subplot(2, 2, 4)
p_compare = 0.9
K_ratio = 0.1 # 10%的成功比例
n_ratio = 0.2 # 20%的抽样比例
population_sizes = [50, 100, 200, 500, 1000]
thresholds_pop = []
for M_pop in population_sizes:
K_pop = int(M_pop * K_ratio)
n_pop = int(M_pop * n_ratio)
threshold = hygeinv(p_compare, M_pop, K_pop, n_pop)
thresholds_pop.append(threshold)
plt.plot(population_sizes, thresholds_pop, 'ro-', markersize=6)
plt.xlabel('总体大小')
plt.ylabel('90%分位数')
plt.title('总体大小对分位数的影响')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
示例7: 医疗检测应用
疾病筛查: 人群中5%患病，检测100人
M_medical = 10000
K_disease = 500 # 5%患病率
n_test = 100
detection_levels = [0.8, 0.9, 0.95, 0.99]
疾病筛查: 不同检测置信度下的预期最大病例数
for level in detection_levels:
threshold = hygeinv(level, M_medical, K_disease, n_test)
print(f" {level * 100}%置信度: 最多检测到 {int(threshold)} 个病例")
# 80.0%置信度: 最多检测到 7 个病例
# 90.0%置信度: 最多检测到 8 个病例
# 95.0%置信度: 最多检测到 9 个病例
# 99.0%置信度: 最多检测到 11 个病例
示例8: 教育评估
班级中有30名学生，预计60%能达到优秀(18人)，随机抽查10人评估
M_class = 30
K_excellent = 18
n_assess = 10
confidence_edu = [0.8, 0.9, 0.95]
教育评估: 不同置信水平下的优秀学生数量预期
for conf in confidence_edu:
threshold = hygeinv(conf, M_class, K_excellent, n_assess)
print(f" {conf * 100}%置信水平: 最多 {int(threshold)} 个优秀学生")
# 80.0%置信水平: 最多 7 个优秀学生
# 90.0%置信水平: 最多 8 个优秀学生
# 95.0%置信水平: 最多 8 个优秀学生
示例9: 环境监测
湖泊中有1000条鱼，污染导致10%死亡(100条)，抽样检测50条
M_lake = 1000
K_dead = 100
n_monitor = 50
alert_levels = [0.95, 0.99]
环境监测: 污染警报阈值
for level in alert_levels:
threshold = hygeinv(level, M_lake, K_dead, n_monitor)
print(f" {level * 100}%置信度警报: 发现超过 {int(threshold)} 条死鱼时触发")
# 95.0%置信度警报: 发现超过 9 条死鱼时触发
# 99.0%置信度警报: 发现超过 10 条死鱼时触发
示例10: 生产优化应用
优化生产线: 1000个产品中目标有50个瑕疵，抽检80个
M_production = 1000
K_defect_target = 50
n_quality = 80
improvement_levels = [0.1, 0.25, 0.5, 0.75, 0.9]
生产优化: 质量改进目标
for improvement in improvement_levels:
# 改进后的缺陷数
K_improved = int(K_defect_target * (1 - improvement))
threshold_improved = hygeinv(0.95, M_production, K_improved, n_quality)
threshold_current = hygeinv(0.95, M_production, K_defect_target, n_quality)
print(f" 改进{improvement * 100}%: 接受阈值从{int(threshold_current)}降至{int(threshold_improved)}")
# 改进10.0%: 接受阈值从7降至7
# 改进25.0%: 接受阈值从7降至6
# 改进50.0%: 接受阈值从7降至4
# 改进75.0%: 接受阈值从7降至3
# 改进90.0%: 接受阈值从7降至1
示例11: 边界情况测试
boundary_cases = [
(0.0, 100, 20, 10),
(0.0001, 100, 20, 10),
(0.9999, 100, 20, 10),
(1.0, 100, 20, 10),
(0.5, 100, 0, 10),
(0.5, 100, 100, 10),
]
for p_val, M_val, K_val, n_val in boundary_cases:
result = hygeinv(p_val, M_val, K_val, n_val)
print(f" p={p_val}, M={M_val}, K={K_val}, n={n_val} -> {result}")
# p=0.0, M=100, K=20, n=10 -> 0.0
# p=0.0001, M=100, K=20, n=10 -> 0.0
# p=0.9999, M=100, K=20, n=10 -> 7.0
# p=1.0, M=100, K=20, n=10 -> 10.0
# p=0.5, M=100, K=0, n=10 -> 0.0
# p=0.5, M=100, K=100, n=10 -> 10.0
示例12: 实际决策支持
客户投诉处理: 1000个客户中有50个可能投诉，随机联系100个了解情况
M_customers = 1000
K_complaints = 50
n_contact = 100
客户投诉风险评估:
risk_levels = [0.05, 0.1, 0.2]
for risk in risk_levels:
threshold = hygeinv(1 - risk, M_customers, K_complaints, n_contact)
print(f" 可接受{risk * 100}%风险: 最多遇到 {int(threshold)} 个投诉客户")
if threshold >= 10: # 假设10个投诉是处理上限
print(f" ⚠️ 警告: 可能超出处理能力!")
# 可接受5.0%风险: 最多遇到 9 个投诉客户
# 可接受10.0%风险: 最多遇到 8 个投诉客户
# 可接受20.0%风险: 最多遇到 7 个投诉客户
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.stats import hypergeom
# 定义计算单个元素的逆CDF函数
def compute_inv_cdf(p_val, M_val, K_val, n_val):
# 参数验证
if not (0 <= p_val <= 1):
raise ValueError(f"概率p必须在[0,1]范围内，当前为{p_val}")
# 转换为整数并验证
try:
M_int = int(M_val)
K_int = int(K_val)
n_int = int(n_val)
if M_int != M_val or K_int != K_val or n_int != n_val:
raise ValueError("参数M, K, n必须为整数")
except (TypeError, ValueError):
raise ValueError("参数M, K, n必须为整数")
if M_int <= 0 or K_int < 0 or n_int < 0:
raise ValueError("参数必须满足 M>0, K≥0, n≥0")
if K_int > M_int or n_int > M_int:
raise ValueError("K和n不能超过M")
# 处理边界情况
if p_val <= 0:
return 0
if p_val >= 1:
return min(K_int, n_int)
# 创建分布对象并查找分位数
hg = hypergeom(M_int, K_int, n_int)
max_k = min(K_int, n_int)
for k in range(0, max_k + 1):
if hg.cdf(k) >= p_val - 1e-10: # 浮点容差
return k
return max_k # 理论上不会执行到此
def hypergeometric_inverse_cumulative_distribution(input_str):
"""
计算超几何分布的逆累积分布函数（分位数值），对标Matlab的hygeinv函数。
参数:
input_str: 输入字符串，格式为四元组 (p, M, K, n)，各参数可为矩阵或标量。
返回:
数值、SymPy矩阵或错误信息字符串。
示例:
>>> hypergeometric_inverse_cumulative_distribution("(0.5, 10, 5, 3)")
1
>>> hypergeometric_inverse_cumulative_distribution("([[0.1, 0.5], [0.9, 1.0]], 10, 5, 3)")
Matrix([[1, 1], [2, 3]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 4:
if all(e.is_number for e in expr):
p, M, K, n = expr # 参数顺序: p, M总体大小, K成功数, n抽样数
# 标量计算
try:
p_val = float(p)
M_val = int(M)
K_val = int(K)
n_val = int(n)
result = compute_inv_cdf(p_val, M_val, K_val, n_val)
except Exception as e:
error = True
result = f"参数错误: {str(e)}"
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 标量计算示例
print(hypergeometric_inverse_cumulative_distribution("0.99,1000,10,50"))
# 3
print(hypergeometric_inverse_cumulative_distribution("0.5,1000,10,50"))
# 0
超几何函数的概率密度函数
Y = hygepdf（X，M，K，N）使用相应的总体大小M, 总体中具有所需特征的项目数量K和提取的样本数量N计算X中每个值的超几何pdf. X,M,K和N可以是向量,矩阵,它们都具有相同的大小.标量输入被扩展为与其他输入具有相同维度的常数数组.
M,K和N中的参数必须都是正整数,N≤M. X中的值必须小于或等于所有参数值.
X -- 标量, 符号变量, 向量, 矩阵, 多维数组
M, K, N -- 数值标量, 向量, 矩阵, 多维数组
import matplotlib.pyplot as plt
示例1: 基本的标量计算
result = hygepdf(2, 100, 20, 10)
print(f"hygepdf(2, 100, 20, 10) = {result}")
#hygepdf(2, 100, 20, 10) = 0.3181706296100318
#解释: 从100个物品中(20个目标)，抽10个，恰好抽到2个目标物品的概率
示例2: 质量检测应用
场景: 一批100个产品中有5个次品，随机抽取10个检测
M = 100
K = 5 # 次品数
n = 10 # 抽样数
质量检测: 100个产品中5个次品，抽检10个
for x in range(0, 6):
prob = hygepdf(x, M, K, n)
print(f" 恰好发现 {x} 个次品的概率: {float(prob):.6f}")
# 恰好发现 0 个次品的概率: 0.583752
# 恰好发现 1 个次品的概率: 0.339391
# 恰好发现 2 个次品的概率: 0.070219
# 恰好发现 3 个次品的概率: 0.006384
# 恰好发现 4 个次品的概率: 0.000251
# 恰好发现 5 个次品的概率: 0.000003
示例3: 生物统计 - 稀有物种研究
场景: 森林中有200只动物，其中10只是目标物种，随机观察30只
M_bio = 200
K_bio = 10
n_bio = 30
稀有物种研究: 200只动物中10只目标物种，观察30只
observation_counts = [0, 1, 2, 3, 4, 5]
for x in observation_counts:
prob = hygepdf(x, M_bio, K_bio, n_bio)
print(f" 恰好观察到 {x} 只目标物种的概率: {float(prob):.6f}")
# 恰好观察到 0 只目标物种的概率: 0.188941
# 恰好观察到 1 只目标物种的概率: 0.352065
# 恰好观察到 2 只目标物种的概率: 0.283608
# 恰好观察到 3 只目标物种的概率: 0.129914
# 恰好观察到 4 只目标物种的概率: 0.037430
# 恰好观察到 5 只目标物种的概率: 0.007078
示例4: 社会调查 - 选民倾向
场景: 城市有5000选民，其中2000支持某政策，随机调查100人
M_survey = 5000
K_survey = 2000
n_survey = 100
选民倾向调查: 5000选民中2000支持，调查100人
support_levels = [35, 40, 45, 50, 55, 60]
for x in support_levels:
prob = hygepdf(x, M_survey, K_survey, n_survey)
print(f" 恰好 {x} 人支持的概率: {float(prob):.8f}")
# 恰好 35 人支持的概率: 0.04908623
# 恰好 40 人支持的概率: 0.08204364
# 恰好 45 人支持的概率: 0.04780618
# 恰好 50 人支持的概率: 0.01001614
# 恰好 55 人支持的概率: 0.00076205
# 恰好 60 人支持的概率: 0.00002085
示例5: 数组输入计算
x_values = [0, 1, 2, 3, 4, 5]
result_array = hygepdf(x_values, 50, 10, 5)
print("x值:", x_values)
print("概率密度:", [f"{p:.6f}" for p in result_array])
#x值: [0, 1, 2, 3, 4, 5]
#概率密度: ['0.310563', '0.431337', '0.209840', '0.044177', '0.003965', '0.000119']
示例6: 可视化 - 概率质量函数
plt.figure(figsize=(15, 10))
# 子图1: 不同总体大小下的概率分布
plt.subplot(2, 2, 1)
M_values = [20, 50, 100]
K_ratio = 0.2 # 成功比例
n_ratio = 0.3 # 抽样比例
for M in M_values:
K = int(M * K_ratio)
n = int(M * n_ratio)
x_range = range(0, min(n, K) + 1)
pmf_values = []
for x in x_range:
prob = hygepdf(x, M, K, n)
pmf_values.append(float(prob))
plt.plot(x_range, pmf_values, 'o-', label=f'M={M}, K={K}, n={n}')
plt.xlabel('x (成功次数)')
plt.ylabel('概率密度 P(X=x)')
plt.title('不同总体大小下的超几何分布')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图2: 不同抽样比例的影响
plt.subplot(2, 2, 2)
M_fixed = 100
K_fixed = 20
n_ratios = [0.1, 0.3, 0.5, 0.7]
for ratio in n_ratios:
n = int(M_fixed * ratio)
x_range = range(0, min(n, K_fixed) + 1)
pmf_values = []
for x in x_range:
prob = hygepdf(x, M_fixed, K_fixed, n)
pmf_values.append(float(prob))
plt.plot(x_range, pmf_values, 'o-', label=f'n={n} (ratio={ratio})')
plt.xlabel('x (成功次数)')
plt.ylabel('概率密度 P(X=x)')
plt.title('不同抽样比例下的超几何分布')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图3: 质量控制中的概率分布
plt.subplot(2, 2, 3)
# 场景: 不同缺陷率下的质量分布
M_quality = 1000
n_test = 50
defect_rates = [0.01, 0.02, 0.05, 0.1] # 次品率
x_range_quality = range(0, 11)
for rate in defect_rates:
K_defects = int(M_quality * rate)
pmf_values = []
for x in x_range_quality:
prob = hygepdf(x, M_quality, K_defects, n_test)
pmf_values.append(float(prob))
plt.plot(x_range_quality, pmf_values, 'o-',
label=f'次品率={rate * 100}%', markersize=3)
plt.xlabel('检测到的次品数')
plt.ylabel('概率密度')
plt.title('质量控制: 不同缺陷率的概率分布')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图4: 与二项分布的比较
plt.subplot(2, 2, 4)
from scipy.stats import binom
M_large = 10000
K_large = 2000
n_large = 100
p_binom = K_large / M_large # 二项分布的成功概率
x_range = range(0, 31)
hyper_pmf = []
binom_pmf = []
for x in x_range:
# 超几何分布
prob_hyper = hygepdf(x, M_large, K_large, n_large)
hyper_pmf.append(float(prob_hyper))
# 二项分布近似
prob_binom = binom.pmf(x, n_large, p_binom)
binom_pmf.append(prob_binom)
plt.plot(x_range, hyper_pmf, 'bo-', label='超几何分布', markersize=4)
plt.plot(x_range, binom_pmf, 'r--', label='二项分布近似', linewidth=2)
plt.xlabel('x')
plt.ylabel('概率密度')
plt.title('超几何分布 vs 二项分布近似 (大总体)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
示例7: 医学研究 - 疾病检测
场景: 人群中有2%患病，检测50人
M_medical = 10000
K_disease = 200 # 2%患病率
n_test = 50
疾病检测: 10000人中200人患病，检测50人
case_counts = [0, 1, 2, 3, 4, 5]
for x in case_counts:
prob = hygepdf(x, M_medical, K_disease, n_test)
print(f" 恰好检测到 {x} 个病例的概率: {float(prob):.6f}")
#恰好检测到 0 个病例的概率: 0.363257
#恰好检测到 1 个病例的概率: 0.372533
#恰好检测到 2 个病例的概率: 0.186248
#恰好检测到 3 个病例的概率: 0.060498
#恰好检测到 4 个病例的概率: 0.014357
#恰好检测到 5 个病例的概率: 0.002654
示例8: 库存管理
场景: 仓库有500件商品，其中20件有瑕疵，随机检查25件
M_inventory = 500
K_defective = 20
n_inspect = 25
库存检查: 500件商品中20件瑕疵，检查25件
defect_found = [0, 1, 2, 3, 4, 5]
for x in defect_found:
prob = hygepdf(x, M_inventory, K_defective, n_inspect)
print(f" 恰好发现 {x} 件瑕疵品的概率: {float(prob):.6f}")
# 恰好发现 0 件瑕疵品的概率: 0.351194
# 恰好发现 1 件瑕疵品的概率: 0.385081
# 恰好发现 2 件瑕疵品的概率: 0.192119
# 恰好发现 3 件瑕疵品的概率: 0.057887
# 恰好发现 4 件瑕疵品的概率: 0.011792
# 恰好发现 5 件瑕疵品的概率: 0.001723
示例9: 教育评估
场景: 班级有30名学生，其中12名数学优秀，随机抽查8名
M_class = 30
K_excellent = 12
n_sample = 8
学生能力评估: 30人中12人优秀，抽查8人
excellent_counts = range(0, 9)
for x in excellent_counts:
prob = hygepdf(x, M_class, K_excellent, n_sample)
print(f" 恰好 {x} 人优秀的概率: {float(prob):.6f}")
# 恰好 0 人优秀的概率: 0.007476
# 恰好 1 人优秀的概率: 0.065247
# 恰好 2 人优秀的概率: 0.209335
# 恰好 3 人优秀的概率: 0.322054
# 恰好 4 人优秀的概率: 0.258794
# 恰好 5 人优秀的概率: 0.110419
# 恰好 6 人优秀的概率: 0.024154
# 恰好 7 人优秀的概率: 0.002436
# 恰好 8 人优秀的概率: 0.000085
示例10: 符号计算和理论分析
symbolic_result = hygepdf(x, M, K, n)
超几何分布的符号表达式:
print(symbolic_result)
#binomial(K, x)*binomial(-K + M, n - x)/binomial(M, n)
示例11: 选举预测
场景: 选区有1000选民，其中550支持候选人A，随机调查100人
M_election = 1000
K_support = 550
n_poll = 100
print("选举预测: 1000选民中550人支持A，调查100人")
support_levels = [45, 50, 55, 60, 65]
for x in support_levels:
prob = hygepdf(x, M_election, K_support, n_poll)
print(f" 恰好 {x} 人支持A的概率: {float(prob):.8f}")
# 恰好 45 人支持A的概率: 0.00907471
# 恰好 50 人支持A的概率: 0.04804026
# 恰好 55 人支持A的概率: 0.08431185
# 恰好 60 人支持A的概率: 0.04858466
# 恰好 65 人支持A的概率: 0.00891474
示例12: 环境监测
场景: 湖泊有800条鱼，污染导致40条死亡，抽样检测60条
M_lake = 800
K_dead = 40
n_sample = 60
环境污染监测: 800条鱼中40条死亡，检测60条
dead_found = [0, 1, 2, 3, 4, 5, 6, 7, 8]
for x in dead_found:
prob = hygepdf(f"({x}, {M_lake}, {K_dead}, {n_sample})")
print(f" 恰好发现 {x} 条死鱼的概率: {float(prob):.6f}")
# 恰好发现 0 条死鱼的概率: 0.040749
# 恰好发现 1 条死鱼的概率: 0.139510
# 恰好发现 2 条死鱼的概率: 0.228642
# 恰好发现 3 条死鱼的概率: 0.238941
# 恰好发现 4 条死鱼的概率: 0.178951
# 恰好发现 5 条死鱼的概率: 0.102345
# 恰好发现 6 条死鱼的概率: 0.046509
# 恰好发现 7 条死鱼的概率: 0.017254
# 恰好发现 8 条死鱼的概率: 0.005328
示例13: 最可能值分析
分析在给定参数下最可能出现的成功次数
M_analysis = 100
K_analysis = 30
n_analysis = 20
max_prob = 0
most_likely = 0
分析: M=100, K=30, n=20
for x in range(max(0, n_analysis + K_analysis - M_analysis), min(K_analysis, n_analysis) + 1):
prob = hygepdf(x, M_analysis, K_analysis, n_analysis)
prob_val = float(prob)
print(f" P(X={x}) = {prob_val:.6f}")
if prob_val > max_prob:
max_prob = prob_val
most_likely = x
# P(X=0) = 0.000302
# P(X=1) = 0.003553
# P(X=2) = 0.018826
# P(X=3) = 0.059674
# P(X=4) = 0.126808
# P(X=5) = 0.191826
# P(X=6) = 0.214091
# P(X=7) = 0.180287
# P(X=8) = 0.116176
# P(X=9) = 0.057760
# P(X=10) = 0.022238
# P(X=11) = 0.006628
# P(X=12) = 0.001523
# P(X=13) = 0.000268
# P(X=14) = 0.000036
# P(X=15) = 0.000004
# P(X=16) = 0.000000
# P(X=17) = 0.000000
# P(X=18) = 0.000000
# P(X=19) = 0.000000
# P(X=20) = 0.000000
print(f"最可能值: X = {most_likely} (概率 = {max_prob:.6f})")
expected = n_analysis * K_analysis / M_analysis
print(f"期望值: {expected:.2f}")
#最可能值: X = 6 (概率 = 0.214091)
#期望值: 6.00
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Hypergeometric, density
from scipy.stats import hypergeom
def hypergeometric_probability_density(input_str):
"""
计算超几何分布概率密度函数，对标Matlab的hygepdf函数。
参数:
input_str: 输入字符串，格式为四元组 (x, M, K, n)，各参数可为矩阵或标量。
参数说明:
- x: 成功次数
- M: 总体大小
- K: 总体中成功数量
- n: 抽样数量
返回:
数值、SymPy矩阵或错误信息字符串。
示例:
>>> hypergeometric_probability_density("(2, 100, 20, 5)")
0.2022
>>> hypergeometric_probability_density("([0,1,2], 10, 5, 3)")
Matrix([[0.0833, 0.4167, 0.4167]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def eval_hypergeometric_pdf(x_val, M_val, K_val, n_val):
"""计算单个元素的概率密度"""
# 参数验证
# 符号计算直接返回表达式
X = Hypergeometric('X', M_val, K_val, n_val)
return density(X)(x_val)
if isinstance(expr, tuple) and len(expr) == 4:
if all(e.is_integer for e in expr):
params = tuple(int(e) for e in expr)
"""
计算超几何分布的概率质量函数（PMF），对标Matlab的hygepdf函数
参数:
x_val (int): 成功次数
M_val (int): 总体大小
K_val (int): 总体中成功的元素个数
n_val (int): 抽取的样本大小
返回:
float: 超几何分布在x_val处的概率密度值
"""
# 直接调用scipy.stats.hypergeom.pmf计算概率密度
result = hypergeom.pmf(*params)
elif any(e.free_symbols for e in expr):
"""计算单个元素的概率密度"""
# 参数验证
X = Hypergeometric('X', *expr[1:])
result = density(X)(expr[0])
else:
error = True
return result if not error else f"输入错误:{input_str}"
except Exception as e:
return f"系统错误: {str(e)}"
# 测试示例
if __name__ == "__main__":
# 标量测试
print(hypergeometric_probability_density("(2, 100, 20, 5)"))
# 输出: 0.20734379349990548
# 符号计算测试
print(hypergeometric_probability_density("(k, 10, 5, 3)"))
# 输出: binomial(5, k)*binomial(5, 3 - k)/120
超几何函数的随机数
R = hygernd(M,K,N)从具有相应总体大小的超几何分布中生成随机数,M,总体中具有所需特征的项目数量K, 以及提取的样本数量N. M,K和N可以是向量,矩阵,它们都具有相同的大小, 这也是R的大小.
M,K或N的标量输入被扩展为与其他输入具有相同维度的常数数组.
R=hyernd(M,K,N,M,N)生成一个M-by-N阵列. M,K,N参数可以是标量或与R大小相同的数组.
M, K, N -- 数值标量, 向量, 矩阵, 多维数组
示例1. 产品质量抽检
场景：从100个产品中抽检20个，已知有15个次品
产品质量抽检:
print(hyernd(100, 15, 20)) # 单次抽检发现的次品数
#[2]
print(hyernd(100, 15, 20, 5)) # 5次抽检的结果
#[2 3 2 3 3]
print(hyernd(100, 15, 20, 3, 2)) # 3×2的抽检结果矩阵
#[[3 3]
[4 6]
[5 6]]
示例2. 生物多样性调查
场景：森林中有200只动物，标记了30只，重新捕获25只进行调查
print(hyernd(200, 30, 25)) # 重捕到的标记动物数量
#[0]
print(hyernd(200, 30, 25, 10)) # 10次调查的结果
#[2 1 2 7 5 2 6 5 1 4]
示例3. 选举民意调查
场景：城市有50000选民，其中18000支持A候选人，调查500人
print(hyernd(50000, 18000, 500)) # 调查中支持A的人数
#[166]
print(hyernd(50000, 18000, 500, 8)) # 8次独立调查结果
#[177 176 210 205 184 195 172 166]
示例4. 医学临床试验
场景：药物试验，100名患者中40人有副作用，随机观察30名
print(hyernd(100, 40, 30)) # 观察到有副作用的患者数
#[13]
print(hyernd(100, 40, 30, 6)) # 6个医疗中心的观察结果
#[12 10 8 12 11 12]
示例5. 库存管理
场景：仓库有不同批次的商品，检查各批次的缺陷品
M_matrix = [[80, 120], [95, 105]] # 各批次商品总数
K_matrix = [[8, 15], [12, 10]] # 各批次缺陷品数
N_matrix = [[10, 15], [12, 8]] # 各批次抽检数
库存管理:
result = hyernd(M_matrix, K_matrix, N_matrix)
各批次发现的缺陷品数:
print(result)
#[[0,2]
[1,2]]
示例6. 教育质量评估
场景：从不同班级随机抽取学生进行测试
全校1000名学生，优秀生200名，随机抽取50名
print(f"优秀生数量: {hyernd(1000, 200, 50)}")
#优秀生数量: [8]
5个班级的模拟结果
print(hyernd(1000, 200, 50, 5))
#[11 7 9 6 9]
示例7. 网络安全检测
场景：检测网络流量中的恶意数据包
10000个数据包中有150个恶意包，随机检查200个
print(f"检测到的恶意包: {hyernd(10000, 150, 200)}")
#[4]
8小时连续监测结果
print(hyernd(10000, 150, 200, 8))
#[3 4 2 1 6 6 5 2]
示例8. 多维度应用场景
复杂场景：不同部门不同时间段的抽检
各部门参数矩阵
departments_M = [[50, 75, 60], [80, 45, 90]] # 各部门物品总数
departments_K = [[5, 8, 6], [12, 4, 9]] # 各部门问题物品数
departments_N = [[10, 15, 12], [16, 8, 18]] # 各部门抽检数
result = hyernd(departments_M, departments_K, departments_N)
各部门抽检发现的问题数:
print(result)
#[[0,2,1]
[1,2,0]]
示例9. 科学研究中的随机抽样
场景：基因研究中稀有基因的发现
1000个基因中有50个目标基因，随机测序100个
for i in range(5):
result = hyernd(1000, 50, 100)
print(f"实验{i + 1}发现的目标基因数: {result}")
#实验1发现的目标基因数: [6]
#实验2发现的目标基因数: [4]
#实验3发现的目标基因数: [2]
#实验4发现的目标基因数: [6]
#实验5发现的目标基因数: [8]
示例10. 市场调研
场景：不同城市消费者偏好调查
各城市参数
cities_population = [[800, 1200], [600, 900]] # 各城市样本量
cities_preference = [[240, 360], [180, 270]] # 偏好某产品的人数
cities_survey = [[80, 120], [60, 90]] # 调查人数
result = hyernd(cities_population, cities_preference, cities_survey)
各城市调查中偏好该产品的人数:
print(result)
#[[18,38]
[17,30]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.stats import hypergeom
import numpy as np
def hypergeometric_random_number(input_str):
"""
生成超几何分布随机数，对标Matlab的hygernd函数。
参数格式:
(M, K, N) -> 单个随机数
(M, K, N, size) -> 指定形状的数组
(M_mat, K_mat, N_mat) -> 矩阵参数，逐元素生成随机数
参数说明:
M: 总体大小（必须为正整数）
K: 成功项数量（0 ≤ K ≤ M）
N: 抽样数量（0 ≤ N ≤ M）
size: 输出形状（整数或元组）
返回:
随机数、SymPy矩阵、numpy数组或错误信息
示例:
>>> hypergeometric_random_number("(100, 20, 5)")
3
>>> hypergeometric_random_number("(100, 20, 5, 3)")
Matrix([[2, 4, 3]])
>>> hypergeometric_random_number("([[90,100],[110,120]], 20, 5)")
Matrix([[2, 3], [1, 4]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def generate_rvs(M_val, K_val, N_val, size=None):
"""核心生成函数"""
# 参数验证
M = int(M_val)
K = int(K_val)
N = int(N_val)
if M <= 0 or K < 0 or N < 0:
raise ValueError("参数必须满足 M>0, K≥0, N≥0")
if K > M or N > M:
raise ValueError("K和N不能超过M")
# 生成随机数
if size:
arr = hypergeom.rvs(M, K, N, size=size)
return sp.Matrix(arr) if isinstance(arr, np.ndarray) else arr
else:
return hypergeom.rvs(M, K, N)
# 解析输入参数
if isinstance(expr, tuple):
# 处理标量参数+size参数 ------------------------------------------------
if len(expr) == 3: # 单个随机数
return generate_rvs(*expr)
elif len(expr) >= 4: # 带size参数
M, K, N = expr[:3]
size = tuple(int(x) for x in expr[3:]) if len(expr) > 4 else int(expr[3])
return generate_rvs(M, K, N, size=size)
raise ValueError("输入格式错误，应为 (M,K,N) 或 (M,K,N,size)")
except Exception as e:
return f"错误: {str(e)}"
# 测试示例
if __name__ == "__main__":
# 标量测试
print(hypergeometric_random_number("(100, 20, 5)"))
# 0
# 一维数组测试
print(hypergeometric_random_number("(100, 20, 5, 3)"))
# Matrix([[2],
# [1],
# [0]])
超几何函数的均值和方差
[MN，V] = hygestat（M, K, N）返回具有相应总体大小的超几何分布的均值和方差, M, 总体中具有所需特征的项目数量K, 以及提取的样本数量N. M,K和N的向量或矩阵输入必须具有相同的大小,这也是MN和V的大小.
M,K或N的标量输入被扩展为与其他输入具有相同维度的常数矩阵.
具有参数M,K和N的超几何分布的平均值为NK/M,方差为NK(M-K)(M-N)/[M^2(M-1)].
M, K, N -- 数值标量, 向量, 矩阵.
示例1：质量控制场景
在生产线质量检测中，从100个产品中随机抽取10个，已知有15个次品
mean, var = hygestat(100, 15, 10)
print(f"质量控制 - 平均抽到次品数: {mean}, 方差: {var}")
#质量控制 - 平均抽到次品数: [1.5], 方差: [1.15909091]
示例2：生物多样性研究
在生态调查中，从50只鸟中捕捉20只，已知其中有标记的鸟有8只
mean, var = hygestat(50, 8, 20)
print(f"生态调查 - 平均捕捉标记鸟数: {mean}, 方差: {var}")
#生态调查 - 平均捕捉标记鸟数: [3.2], 方差: [1.64571429]
示例3：选举预测
在选区中，从10000名选民中抽样500人，已知支持某候选人的选民有4500人
mean, var = hygestat(10000, 4500, 500)
print(f"选举预测 - 平均支持者数: {mean}, 标准差: {var ** 0.5}")
#选举预测 - 平均支持者数: [225.], 标准差: [10.84316639]
示例4：库存管理
仓库中有200件商品，其中30件有瑕疵，随机检查25件
mean, var = hygestat(200, 30, 25)
print(f"库存检查 - 平均瑕疵品数: {mean}, 方差: {var}")
#库存检查 - 平均瑕疵品数: [3.75], 方差: [2.80307789]
示例5：科学研究应用
基因研究中，从1000个基因中筛选50个，已知有致病风险的基因有80个
mean, var = hygestat(1000, 80, 50)
print(f"基因研究 - 平均风险基因数: {mean}, 95%置信区间: [{mean - 2 * var ** 0.5}, {mean + 2 * var ** 0.5}]")
#基因研究 - 平均风险基因数: [4.], 95%置信区间: [[0.25861015], [7.74138985]]
示例6：教育评估
从300名学生中随机选取30名评估，已知优秀学生有45名
mean, var = hygestat(300, 45, 30)
print(f"教育评估 - 平均优秀学生数: {mean}, 变异系数: {(var ** 0.5 / mean) * 100}%")
#教育评估 - 平均优秀学生数: [4.5], 变异系数: [41.29994696]%
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.stats import hypergeom
def hypergeometric_mean_variance(input_str):
"""
计算超几何分布的均值与方差，对标Matlab的hygestat函数。
参数格式:
(M, K, N) -> 返回标量均值与方差
(M_mat, K_mat, N_mat) -> 返回同形状的均值矩阵和方差矩阵
参数说明:
M: 总体大小（必须为正整数）
K: 成功项数量（0 ≤ K ≤ M）
N: 抽样数量（0 ≤ N ≤ M）
返回:
(均值, 方差) 元组，或错误信息字符串
示例:
>>> hypergeometric_mean_variance("(100, 20, 5)")
(4.0, 1.2626262626262625)
>>> hypergeometric_mean_variance("([100, 150], [20, 30], [5, 10])")
(Matrix([[4.0, 6.0]]), Matrix([[1.2626, 2.2222]]))
>>> hypergeometric_mean_variance("(50, 60, 30)")
'错误：K必须在0到M之间'
"""
try:
expr = sp.sympify(input_str)
mean_matrix, var_matrix = None, None
error = False
def compute_stats(M_val, K_val, N_val):
"""核心计算函数"""
# 参数验证
try:
M = int(M_val)
K = int(K_val)
N = int(N_val)
except (TypeError, ValueError):
raise ValueError("参数M、K、N必须为整数")
if M <= 0:
raise ValueError("M必须为正整数")
if not (0 <= K <= M):
raise ValueError("K必须在0到M之间")
if not (0 <= N <= M):
raise ValueError("N必须在0到M之间")
# 计算统计量
dist = hypergeom(M, K, N)
return dist.mean(), dist.var()
# 解析输入参数 -----------------------------------------------------------
if isinstance(expr, tuple) and len(expr) == 3:
M, K, N = expr # 参数顺序对应Matlab的hygestat(M, K, N)
# 标量计算模式
try:
mean, var = compute_stats(M, K, N)
return (mean, var)
except Exception as e:
return f"错误：{str(e)}"
else:
return "输入格式错误：应为三元素元组 (M, K, N)"
except Exception as e:
return f"系统错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 标量测试
print(hypergeometric_mean_variance("(100, 20, 5)"))
# (1.0, 0.7676767676767676)
广义超几何函数
hypergeom(a,b,z) 根据输入是浮点还是符号，超几何返回浮点或符号结果。
计算这些数字的超几何函数。因为这些数字是浮点数，所以hypergeom返回浮点结果。
a —— 超几何函数的上参数,数字,向量,符号数字,符号变量,符号表达式,符号函数,符号向量
b —— 超几何函数的下参数,数字,向量,符号数字,符号变量,符号表达式,符号函数,符号向量
z —— 标量,数值,符号变量
示例1：基本数学函数表示
指数函数 e^x 的超几何表示
result1 = hypergeom([], [], x)
print(f"e^x = {result1}")
#e^x = exp(x)
正弦函数的超几何表示
result2 = hypergeom([], [1/2], -x**2/4)
print(f"sin(x) = {result2}")
#sin(x) = cos(x)
贝塞尔函数的超几何表示
result3 = hypergeom([], [1], x)
print(f"J_0(x) = {result3}")
#J_0(x) = besseli(0, 2*sqrt(x))
示例2：概率统计应用
二项分布矩生成函数的超几何表示
result4 = hypergeom([-n], [], -x)
print(f"二项分布MGF: {result4}")
#二项分布MGF: (x + 1)**n
泊松分布相关函数
result5 = hypergeom([1], [2], x)
print(f"泊松相关: {result5}")
#泊松相关: exp(x)/x - 1/x
示例3：物理学应用
量子力学 - 谐振子波函数
result6 = hypergeom([-n], [1/2], x**2)
print(f"谐振子波函数: {result6}")
#谐振子波函数: hyper((-n,), (1/2,), x**2)
电磁学 - 勒让德多项式
result7 = hypergeom([-n, n+1], [1], (1-x)/2)
print(f"勒让德多项式: {result7}")
#勒让德多项式: hyper((-n, n + 1), (1,), 1/2 - x/2)
示例4：工程计算
热传导方程的解
result8 = hypergeom([1/2], [3/2], -x**2)
print(f"误差函数: {result8}")
#误差函数: sqrt(pi)*erf(x)/(2*x)
流体力学 - 边界层方程
result9 = hypergeom([1/3], [4/3], x)
print(f"Airy函数相关: {result9}")
#Airy函数相关: exp(-I*pi/3)*lowergamma(1/3, x*exp_polar(I*pi))/(3*x**(1/3))
示例5：数值计算验证
验证特殊值
import math
计算 log(1+x)/x 在 x=1 处的值
result10 = hypergeom([1, 1], [2], -1)
print(f"log(2) = {result10} = {float(result10.evalf())}")
print(f"验证: math.log(2) = {math.log(2)}")
#log(2) = log(2) = 0.6931471805599453
#验证: math.log(2) = 0.6931471805599453
计算 arctan(x)/x 在 x=1 处的值
result11 = hypergeom([1/2], [3/2], -1)
print(f"π/4 = {result11} = {float(result11.evalf())}")
print(f"验证: math.pi/4 = {math.pi / 4}")
#π/4 = sqrt(pi)*erf(1)/2 = 0.746824132812427
#验证: math.pi/4 = 0.7853981633974483
示例6：组合数学应用
二项式系数生成函数
result12 = hypergeom([-n], [], x)
print(f"二项式定理: {result12}")
#二项式定理: (1 - x)**n
超几何分布的概率生成函数
result13 = hypergeom([-a, -b], [c], 1)
print(f"超几何分布PGF: {result13}")
#超几何分布PGF: gamma(c)*gamma(a + b + c)/(gamma(a + c)*gamma(b + c))
示例7：复杂参数计算
复数参数的计算
result14 = hypergeom([1, 2], [3], 0.5)
print(f"数值计算: {result14.evalf()}")
#数值计算: 1.54517744447956
符号计算保持精确
result15 = hypergeom([1, 2], [3], x)
print(f"符号结果: {result15}")
#符号结果: -2/x - 2*log(1 - x)/x**2
分数参数
result16 = hypergeom([1/2, 1/3], [4/3], 1)
print(f"分数参数: {result16} = {result16.evalf()}")
#分数参数: sqrt(pi)*gamma(4/3)/gamma(5/6) = 1.40218210532545
示例8：级数展开验证
验证超几何级数展开
result17 = hypergeom([1, 2], [3], x)
级数展开
series_expansion = sp.series(result17, x, 0, 6)
print(f"级数展开: {series_expansion}")
#级数展开: 1 + 2*x/3 + x**2/2 + 2*x**3/5 + x**4/3 + 2*x**5/7 + O(x**6)
直接计算前几项验证
manual_sum = 1 + (1 * 2) / (1 * 3) * x + (1 * 2 * 3) / (1 * 2 * 3 * 4) * x ** 2 / 2
print(f"手动计算前3项: {manual_sum}")
#手动计算前3项: 0.125*x**2 + 0.666666666666667*x + 1
示例9：特殊函数关系
合流超几何函数（Kummer函数）
result18 = hypergeom([a], [b], x)
print(f"Kummer函数: {result18}")
#Kummer函数: hyper((a,), (b,), x)
高斯超几何函数
result19 = hypergeom([a, b], [c], x)
print(f"高斯超几何函数: {result19}")
#高斯超几何函数: hyper((a, b), (c,), x)
广义超几何函数 pFq
result20 = hypergeom([1, 1, 1], [2, 2], x)
print(f"3F2函数: {result20}")
#3F2函数: polylog(2, x)/x
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def hypergeometric_function_construct(input_str):
"""
构造并简化广义超几何函数，对标MATLAB的hypergeom函数。
参数:
input_str (str): 输入字符串，格式应为'([a1, a2, ...], [b1, b2, ...], z)'，其中：
- a1, a2, ... 是上参数列表
- b1, b2, ... 是下参数列表
- z 是变量或数值
返回:
Sympy表达式: 简化后的超几何函数表达式，或错误信息。
"""
try:
# 将输入字符串解析为Sympy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为三元组（上参数，下参数，z）
if isinstance(expr, tuple) and len(expr) == 3:
# 处理上参数列表
a = tuple(expr[0]) if isinstance(expr[0], list) else (expr[0],)
# 处理下参数列表
b = tuple(expr[1]) if isinstance(expr[1], list) else (expr[1],)
z = expr[2]
# 构造超几何函数
hyper_func = sp.hyper(a, b, z)
# 尝试展开为更简单的表达式
result = sp.hyperexpand(hyper_func)
else:
error = True
return result if not error else f"输入格式错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：符号参数，返回未简化形式
input1 = '(a,[],x)'
result1 = hypergeometric_function_construct(input1)
print(f"示例1结果: {result1}")
# (1 - x)**(-a)
# 示例2：可简化的超几何函数
input2 = '([1, 1], [2], -x)'
result2 = hypergeometric_function_construct(input2)
print(f"示例2结果: {result2}")
# log(x + 1)/x
# 示例3：数值参数，返回符号表达式（如需要数值结果可调用.evalf()）
input3 = '([1, 1], [2], -1)'
result3 = hypergeometric_function_construct(input3)
print(f"示例3符号结果: {result3}")
# log(2)
print(f"示例3数值结果: {result3.evalf()}")
# 0.693147180559945
平方和的平方根(斜边)
C = hypot(A,B)返回以下方程的计算结果(计算是为了避免下溢或上溢)： C = sqrt(abs(A)^2 + abs(B)^2)
A,B — 输入数组,标量,向量,矩阵,多维数组
示例1：几何计算
计算直角三角形的斜边
result1 = hypot(3, 4)
print(f"直角三角形斜边: {result1}")
#直角三角形斜边: 5.0
计算三维空间中的对角线长度
result2 = hypot(3, 4)
z = 12
distance_3d = np.sqrt(result1**2 + z**2)
print(f"三维空间对角线: {distance_3d}")
#三维空间对角线: 13.0
示例2：物理学应用
计算力的合成（两个垂直分力）
force_x = 8 # x方向分力
force_y = 6 # y方向分力
result3 = hypot(force_x, force_y)
print(f"合力大小: {result3} N")
#合力大小: 10.0 N
计算速度合成
v_x = 3 # x方向速度
v_y = 4 # y方向速度
result4 = hypot(v_x, v_y)
print(f"合速度大小: {result4} m/s")
#合速度大小: 5.0 m/s
示例3：工程测量
计算建筑物对角线测量
building_length = 50 # 长度
building_width = 30 # 宽度
result5 = hypot(building_length, building_width)
print(f"建筑物对角线: {result5} 米")
#建筑物对角线: 58.309518948453004 米
计算斜坡的实际长度
horizontal_distance = 100 # 水平距离
vertical_rise = 10 # 垂直升高
result6 = hypot(horizontal_distance, vertical_rise)
print(f"斜坡实际长度: {result6} 米")
#斜坡实际长度: 100.4987562112089 米
示例4：计算机图形学
计算像素点到原点的距离
pixel_x = 1920
pixel_y = 1080
result7 = hypot(pixel_x, pixel_y)
print(f"屏幕对角线像素数: {result7}")
#屏幕对角线像素数: 2202.9071700822983
计算向量长度（多个点）
vectors = ([[1, 2], [3, 4]], [[4, 3], [2, 1]])
result8 = hypot(vectors)
print("向量长度矩阵:")
print(result8)
#向量长度矩阵:
#[[4.12310563,3.60555128]
[3.60555128,4.12310563]]
示例5：统计学应用
计算欧几里得距离（二维数据点）
point1 = [2, 3]
point2 = [5, 7]
diff_x = point2[0] - point1[0]
diff_y = point2[1] - point1[1]
result9 = hypot(diff_x, diff_y)
print(f"两点间欧几里得距离: {result9}")
#两点间欧几里得距离: 5.0
计算误差向量的模
measurement_errors = ([0.1, -0.2, 0.05], [0.15, 0.1, -0.08])
result10 = hypot(measurement_errors)
print(f"误差向量模长: {result10}")
#误差向量模长: [0.18027756 0.2236068 0.09433981]
示例6：电子工程
计算交流电路的阻抗
resistance = 3 # 电阻 (Ω)
reactance = 4 # 电抗 (Ω)
result11 = hypot(resistance, reactance)
print(f"电路阻抗: {result11} Ω")
#电路阻抗: 5.0 Ω
计算信号幅度（I/Q分量）
I_component = 0.8 # 同相分量
Q_component = 0.6 # 正交分量
result12 = hypot(I_component, Q_component)
print(f"信号幅度: {result12}")
#信号幅度: 1.0
示例7：导航和定位
计算GPS坐标间的平面距离（近似）
lat_diff = 0.01 # 纬度差（度）
lon_diff = 0.02 # 经度差（度）
result13 = hypot(lat_diff, lon_diff)
print(f"坐标间近似距离: {result13} 度")
#坐标间近似距离: 0.022360679774997897 度
计算机器人移动路径
movement_x = [1, 2, 3, 4] # x方向移动
movement_y = [2, 3, 1, 5] # y方向移动
result14 = hypot(movement_x, movement_y)
print("各步移动距离:")
print(result14)
#各步移动距离:
#[2.23606798,3.60555128,3.16227766,6.40312424]
示例8：金融分析
计算投资组合的风险（简化模型）
stock_risk = 0.15 # 股票风险
bond_risk = 0.08 # 债券风险
假设股票和债券风险不相关
result15 = hypot(stock_risk, bond_risk)
print(f"投资组合总风险: {result15:.3f}")
#投资组合总风险: 0.170
计算价格变动的幅度
price_changes_x = [1.5, -2.0, 0.8] # 资产A价格变化
price_changes_y = [0.5, 1.2, -0.3] # 资产B价格变化
result16 = hypot(price_changes_x, price_changes_y)
print("价格变动幅度:")
print(result16)
#价格变动幅度:
#[1.58113883,2.33238076,0.85440037]
示例9：批量处理实际数据
处理多个三角形的斜边计算
triangles = ([[3, 6, 9], [5, 8, 12]], [[4, 8, 12], [12, 15, 16]])
result17 = hypot(triangles)
print("多个三角形的斜边:")
print(result17)
#多个三角形的斜边:
#[[5,10,15]
[13,17,20]]
计算复杂形状的对角线
complex_shape = ([10, 20, 30, 40], [15, 25, 35, 45])
result18 = hypot(complex_shape)
print("复杂形状对角线:")
print(result18)
#复杂形状对角线:
#[18.02775638,32.01562119,46.09772229,60.20797289]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def hypotenuse_calculate(input_str):
"""
对标MATLAB hypot函数计算直角三角形的斜边
参数:
input_str (str): SymPy可解析的字符串，格式为：
(n, x) 其中n和x可以是标量或矩阵
返回:
Matrix/float: 计算结果矩阵或标量
"""
try:
expr = sp.sympify(input_str)
error = False
# 检查输入是否为二元组
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入错误: 需要两个参数 {input_str}"
# 参数解析
a, b = expr[0], expr[1]
# 转换为numpy数组（支持标量和矩阵）
def to_ndarray(val):
if isinstance(val, list):
return np.array(val, dtype=float)
elif isinstance(val, sp.Number):
return float(val)
else:
return None
np_a = to_ndarray(a)
np_b = to_ndarray(b)
if np_a is None or np_b is None:
return f"参数类型错误: 需要数值或矩阵 {input_str}"
# 计算斜边（自动处理广播）
result = np.hypot(np_a, np_b)
# 转换回SymPy格式
if isinstance(result, np.ndarray):
return sp.Matrix(result.tolist())
else:
return float(result)
except Exception as e:
return f"错误: {str(e)}"
# 案例1：标量计算
print(hypotenuse_calculate("(3, 4)"))
# → 5.0
# 案例2：矩阵与标量
print(hypotenuse_calculate("([3,4], 4)"))
# Matrix([[5.00000000000000],
# [5.65685424949238]])
# 案例3：矩阵广播
print(hypotenuse_calculate("([[1,2],[3,4]], [5,6])"))
# Matrix([[5.09901951359278, 6.32455532033676],
# [5.83095189484530, 7.21110255092798]])