将角的单位从弧度转换为度
D = rad2deg(R) 将 R 中每个元素的角单位从弧度转换为度.
R — 以弧度为单位的角,标量,向量,矩阵
示例1：常见特殊角度的转换
special_angles = [
0, # 0弧度 = 0度
@pi/6, # π/6弧度 = 30度
@pi/4, # π/4弧度 = 45度
@pi/3, # π/3弧度 = 60度
@pi/2, # π/2弧度 = 90度
@pi, # π弧度 = 180度
3*@pi/2, # 3π/2弧度 = 270度
2*@pi # 2π弧度 = 360度
]
for angle in special_angles:
result = rad2deg(angle)
print(f"{angle:8} 弧度 = {result} 度")
#0 弧度 = 0.0 度
#pi/6 弧度 = 30.000000000000004 度
#pi/4 弧度 = 45.0 度
#pi/3 弧度 = 60.00000000000001 度
#pi/2 弧度 = 90.0 度
#pi 弧度 = 180.0 度
#3*pi/2 弧度 = 270.0 度
#2*pi 弧度 = 360.0 度
示例2：三角函数反函数的输出转换
计算常见三角函数值的反函数并转换为角度
trig_values = [
asin(0.5), # arcsin(0.5) ≈ 0.5236弧度 = 30度
acos(0.5), # arccos(0.5) ≈ 1.0472弧度 = 60度
atan(1), # arctan(1) = π/4弧度 = 45度
asin(sqrt(2)/2) # arcsin(√2/2) = π/4弧度 = 45度
]
for trig_expr in trig_values:
# 先计算弧度值
rad_value = sympify(trig_expr)
# 转换为角度
deg_result = rad2deg(rad_value)
print(f"{trig_expr:15} = {deg_result:.2f} 度")
#asin(0.5) = 30.00 度
#acos(0.5) = 60.00 度
#atan(1) = 45.00 度
#asin(sqrt(2)/2) = 45.00 度
示例3：机器人关节角度转换
机器人逆运动学计算得到的关节弧度值转换为角度
joint_angles_rad = [
[0.7854, 1.0472, 0.5236], # 三个关节的弧度值
[1.5708, 0.0000, 2.0944],
[0.3491, 2.6180, 0.8727]
]
result = rad2deg(joint_angles_rad)
print("机器人关节角度（度）:")
print(result)
#机器人关节角度（度）:
#[[45.0001052295749, 60.0001403060998, 30.0000701530499],
[90.0002104591497, 0.0, 120.0002806122],
[20.001956628017, 150.00035076525, 50.0020267810669]]
示例4：计算机图形学中的旋转矩阵角度
3D变换中的欧拉角（弧度）转换为角度
euler_angles_rad = [
[1.0472, 0.0000, 0.0000], # X轴旋转
[0.0000, 0.7854, 0.0000], # Y轴旋转
[0.0000, 0.0000, 0.5236] # Z轴旋转
]
result = rad2deg(euler_angles_rad)
print("欧拉旋转角度（度）:")
print(result)
#欧拉旋转角度（度）:
#[[60.0001403060998, 0.0, 0.0],
[0.0, 45.0001052295749, 0.0],
[0.0, 0.0, 30.0000701530499]]
示例5：地理坐标系统转换
地理坐标系统中常用弧度表示，但日常使用度
geo_coordinates_rad = [
[0.017453, 0.034907], # 约1°E, 2°N
[0.087266, 0.139626], # 约5°E, 8°N
[0.174533, 0.261799] # 约10°E, 15°N
]
result = rad2deg(geo_coordinates_rad)
print("地理坐标（度）:")
print(result)
#地理坐标（度）:
#[[0.999983239841826, 2.00002377546316],
[4.99997349498864, 7.99998051029363],
[10.0000042857568, 14.9999777807454]]
示例6：物理学的角度量转换
角速度、相位角等物理量从弧度转换为度
physics_angles = [
@pi/12, # 15度 - 小角度近似
@pi/8, # 22.5度
@pi/5, # 36度
2*pi/7" # 约102.857度
]
物理学角度转换:
for angle in physics_angles:
result = rad2deg(angle)
print(f"{angle:8} = {result:.2f} 度")
#pi/12 = 15.00 度
#pi/8 = 22.50 度
#pi/5 = 36.00 度
#2*pi/7 = 51.43 度
示例7：天文学中的角度测量
天体中常用的弧度测量转换为度分秒格式
astronomy_angles = [
0.000290888, # 约1角分（arcminute）
0.000004848, # 约1角秒（arcsecond）
0.008726646, # 约30角分
0.000145444 # 约30角秒
]
天文学角度:
for i, angle in enumerate(astronomy_angles):
deg = rad2deg(angle)
# 转换为度分秒格式
degrees = int(deg)
minutes = (deg - degrees) * 60
seconds = (minutes - int(minutes)) * 60
print(f"{angle:.9f} 弧度 ≈ {degrees}°{int(minutes)}'{seconds:.2f}\"")
#0.000290888 弧度 ≈ 0°0'60.00"
#0.000004848 弧度 ≈ 0°0'1.00"
#0.008726646 弧度 ≈ 0°29'60.00"
#0.000145444 弧度 ≈ 0°0'30.00"
示例8：工程测量中的角度转换
建筑工程、机械工程中的角度测量
engineering_angles = [
[0.05236, 0.10472], # 约3°, 6°
[0.15708, 0.20944], # 约9°, 12°
[0.26180, 0.31416] # 约15°, 18°
]
result = rad2deg(engineering_angles)
print("工程角度测量（度）:")
print(result)
#工程角度测量（度）:
#[[3.00000701530499, 6.00001403060998],
[9.00002104591497, 12.00002806122],
[15.000035076525, 18.0000420918299]]
示例9：复数相位角转换
将复数的相位角（弧度）转换为度
complex_phases = [
atan2(1, 1), # 45度
atan2(1, 0), # 90度
atan2(0, -1), # 180度
atan2(-1, -1) # -135度（或225度）
]
复数相位角:
for phase_expr in complex_phases:
rad_value = sympify(phase_expr)
deg_result = rad2deg(rad_value)
print(f"{phase_expr:12} = {deg_result:.2f} 度")
#atan2(1, 1) = 45.00 度
#atan2(1, 0) = 90.00 度
#atan2(0, -1) = 180.00 度
#atan2(-1, -1) = -135.00 度
示例10：导航系统的航向角转换
将弧度表示的航向角转换为度（0-360度范围）
heading_angles_rad = [
0.0, # 北
@pi/2, # 东
@pi, # 南
3*@pi/2, # 西
0.5236, # 约30度 - 东北偏北
1.0472 # 约60度 - 东北偏东
]
导航航向角:
for angle in heading_angles_rad:
deg_result = rad2deg(str(angle))
print(f"{angle:8} 弧度 = {deg_result:.2f} 度")
#0.0 弧度 = 0.00 度
#pi/2 弧度 = 90.00 度
#pi 弧度 = 180.00 度
#3*pi/2 弧度 = 270.00 度
#0.5236 弧度 = 30.00 度
#1.0472 弧度 = 60.00 度
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def radian_to_degree(input_str):
"""
将输入字符串中的弧度值转换为角度值，支持标量和矩阵。
参数:
input_str (str): 输入的弧度表达式，可以是数字、矩阵或其他形式
返回:
sp.Matrix/float/str: 转换后的结果，若输入不合法则返回错误字符串
示例:
>>> radian_to_degree("1.0471975511965976")
60.0
>>> radian_to_degree("[[0.5235987755982988, 0.7853981633974483], [1.5707963267948966, 3.141592653589793]]")
Matrix([[30.0, 45.0], [90.0, 180.0]])
"""
try:
# 将输入字符串解析为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
def evaluate_rad2deg_function(z):
"""内部函数：处理单个元素的转换"""
if z.is_number:
return np.rad2deg(float(z)) # 使用NumPy的rad2deg函数进行转换
else:
return None # 非数值元素返回None
# 处理元组类型（不支持）
if isinstance(expr, tuple):
error = True
# 处理标量数值
if expr.is_number:
result = evaluate_rad2deg_function(expr)
# 其他类型视为错误
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
if __name__ == "__main__":
# 示例测试用例
test_cases = [
"1.0471975511965976",
# 59.99999999999999
]
for case in test_cases:
print(f"输入：'{case}'")
result = radian_to_degree(case)
print("输出：")
if isinstance(result, sp.Matrix):
# 格式化矩阵输出
print(result)
else:
print(result)
print("-" * 40)
随机均匀矩阵
X = rand(n) 返回一个由均匀分布的随机数组成的 n×n 矩阵
X = rand(m,n) 返回一个由均匀分布的随机数组成的 m×n 矩阵.
X = rand(m,n,min) 返回一个由均匀分布的随机数组成的 m×n 矩阵, 数值在min与inf之间的矩阵.
X = rand(m,n,min,max) 返回一个由均匀分布的随机数组成的 m×n 矩阵, 数值在min与max之间的矩阵.
m, n — 方阵的大小,整数值
min,max,整数值
示例1：蒙特卡洛模拟 - 风险评估
生成1000个随机场景的输入参数矩阵，用于金融风险评估
result1 = rand(1000, 5, 0, 1) # 1000个场景，5个风险因子，每个因子在0-1之间
print(f"蒙特卡洛模拟矩阵: {result1.rows}×{result1.cols}")
#蒙特卡洛模拟矩阵: 1000×5
#应用: 模拟不同风险因子的组合，评估投资组合的潜在损失
示例2：神经网络权重初始化
为神经网络隐藏层初始化权重矩阵
result2 = rand(10, 5, -0.1, 0.1) # 10个神经元到5个神经元的连接，小随机权重
print(f"神经网络权重矩阵: {result2.rows}×{result2.cols}")
#神经网络权重矩阵: 10×5
#应用: 打破对称性，确保神经元学习不同的特征
示例3：物理系统模拟 - 粒子初始位置
在3D空间中随机初始化100个粒子的位置
result3 = rand(100, 3, 0, 10) # 100个粒子，3D坐标，空间范围0-10单位
print(f"粒子位置矩阵: {result3.rows}×{result3.cols}")
#粒子位置矩阵: 100×3
#应用: 模拟气体扩散、流体动力学或天体运动
示例4：图像处理 - 随机噪声添加
生成与图像尺寸相同的随机噪声矩阵
result4 = rand(480, 640, 0, 0.1) # 480×640图像，添加轻微噪声(0-0.1强度)
print(f"图像噪声矩阵: {result4.rows}×{result4.cols}")
#图像噪声矩阵: 480×640
#应用: 数据增强，提高机器学习模型的鲁棒性
示例5：优化算法 - 初始种群生成
为遗传算法生成初始随机种群
result5 = rand(50, 10, -5, 5) # 50个个体，每个有10个参数，搜索范围-5到5
print(f"遗传算法种群矩阵: {result5.rows}×{result5.cols}")
#遗传算法种群矩阵: 50×10
#应用: 在多维参数空间中随机初始化搜索起点
示例6：传感器校准 - 随机测试信号
生成用于传感器校准的随机输入信号
result6 = rand(20, 1, 0, 5) # 20个测试点，电压范围0-5V
print(f"传感器测试信号: {result6.rows}×{result6.cols}")
#传感器测试信号: 20×1
#应用: 测试传感器在不同输入下的响应特性
示例7：游戏开发 - 随机地形生成
生成随机地形高度图
result7 = rand(32, 32, 0, 100) # 32×32地形网格，高度范围0-100单位
print(f"地形高度矩阵: {result7.rows}×{result7.cols}")
#地形高度矩阵: 32×32
#应用: 创建随机但自然的游戏环境
示例8：化学模拟 - 随机浓度分布
模拟化学反应器中不同物质的初始随机浓度
result8 = rand(8, 5, 0, 1) # 8种物质在5个区域的初始浓度(0-1 mol/L)
print(f"化学浓度矩阵: {result8.rows}×{result8.cols}")
#化学浓度矩阵: 8×5
#应用: 研究非均匀初始条件对反应动力学的影响
示例9：经济模型 - 随机需求模拟
模拟不同产品在不同地区的随机需求
result9 = rand(5, 10, 100, 1000) # 5种产品在10个地区的需求(100-1000单位)
print(f"市场需求矩阵: {result9.rows}×{result9.cols}")
#市场需求矩阵: 5×10
#应用: 供应链优化和库存管理决策
示例10：质量控制 - 制造公差模拟
模拟批量生产中产品的随机尺寸变化
result11 = rand(100, 3, -0.05, 0.05) # 100个产品，3个关键尺寸，公差±0.05mm
print(f"产品公差矩阵: {result11.rows}×{result11.cols}")
#产品公差矩阵: 100×3
#应用: 分析制造过程的能力和稳定性
示例11：计算机图形学 - 随机纹理生成
生成随机纹理的基础噪声
result12 = rand(256, 256) # 256×256纹理，灰度值0-1
print(f"纹理噪声矩阵: {result12.rows}×{result12.cols}")
#纹理噪声矩阵: 256×256
#应用: 创建自然材质如木纹、石纹的基础图案
示例12：机器人学 - 随机目标位置
为机器人路径规划算法生成随机目标点
result13 = rand(20, 2, 0, 10) # 20个目标点，2D坐标，工作空间0-10单位
print(f"目标位置矩阵: {result13.rows}×{result13.cols}")
#目标位置矩阵: 20×2
#应用: 测试路径规划算法在不同场景下的性能
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def matrix_rand(input_str):
"""
根据输入字符串生成均匀分布的随机矩阵，对标MATLAB的rand函数，支持扩展参数。
参数：
input_str (str): 表示矩阵维度和范围的字符串。支持以下格式：
- "n"：生成n×n的0~1随机矩阵。
- "(rows, cols)"：生成rows×cols的0~1随机矩阵。
- "(rows, cols, a)"：生成元素在[a, a+1)区间的随机矩阵。
- "(rows, cols, a, b)"：生成元素在[a, b)区间的随机矩阵。
返回：
sympy.Matrix or str: 生成的随机矩阵，若输入无效则返回错误信息。
示例：
>>> matrix_rand("3") # 3x3矩阵，元素0~1
>>> matrix_rand("(2,3)") # 2x3矩阵，元素0~1
>>> matrix_rand("(2,2,5)") # 2x2矩阵，元素5~6
>>> matrix_rand("(3,2,5,10)") # 3x2矩阵，元素5~10
"""
try:
expr = sp.sympify(input_str) # 解析输入字符串为SymPy表达式
error = False
result = None
if isinstance(expr, tuple):
# 处理元组参数
if len(expr) == 4:
# 参数格式：(rows, cols, a, b)，生成[a, b)区间的随机矩阵
rows, cols, a, b = expr
random_matrix = np.random.rand(rows, cols) * (b - a) + a
result = sp.Matrix(random_matrix)
elif len(expr) == 3:
# 参数格式：(rows, cols, a)，生成[a, a+1)区间的随机矩阵
rows, cols, a = expr
random_matrix = np.random.rand(rows, cols) + a
result = sp.Matrix(random_matrix)
elif len(expr) == 2:
# 参数格式：(rows, cols)，生成0~1的随机矩阵
rows, cols = expr
random_matrix = np.random.rand(rows, cols)
result = sp.Matrix(random_matrix)
else:
error = True # 参数长度不符合要求
elif expr.is_integer:
# 参数格式：整数n，生成n×n的0~1随机矩阵
n = int(expr)
random_matrix = np.random.rand(n, n)
result = sp.Matrix(random_matrix)
else:
error = True # 无效的输入类型
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：3x3矩阵，元素0~1
print("示例1:", matrix_rand("3"))
# Matrix([[0.183112236962607, 0.212394217857657, 0.664021485720149],
# [0.669792434360979, 0.0240691129401981, 0.463604975941356],
# [0.288865834582140, 0.330471827282497, 0.923577809087900]])
# 示例2：2x3矩阵，元素0~1
print("示例2:", matrix_rand("(2,3)"))
# Matrix([[0.595556557127196, 0.867659443372348, 0.0504270317200585],
# [0.776355720672051, 0.656865338639271, 0.102445830808784]])
# 示例3：2x2矩阵，元素5~6
print("示例3:", matrix_rand("(2,2,5)"))
# Matrix([[5.61650241924118, 5.87854171655421],
# [5.05252724807980, 5.82949043459748]])
# 示例4：3x2矩阵，元素5~10
print("示例4:", matrix_rand("(3,2,5,10)"))
# Matrix([[9.78184517353882, 8.41039417069258],
# [7.18826504860124, 5.54331991813526],
# [5.50597244331105, 8.61157846090020]])
均匀分布的伪随机整数
X = randi(imax) 返回一个介于 1 和 imax 之间的伪随机整数标量
X = randi(imax,n) 返回 n×n 矩阵，其中包含从区间 [1,imax] 的均匀离散分布中得到的伪随机整数.
X = randi(imax,sz1,...,szN) 返回 sz1×...×szN 数组，其中 sz1,sz2指示每个维度的大小。例如，randi(10,3,4) 返回一个由介于 1 和 10 之间的伪随机整数组成的 3×4 数组.
imax — 示例区间中的最大整数,整数值
imin — 示例区间中的最小整数(默认)
n — 方阵的大小,整数值
示例1：骰子游戏模拟
模拟投掷多个骰子
result1 = randi(6, 5) # 投掷5个骰子，每个骰子1-6点
print(f"骰子投掷结果: {result1}")
#骰子投掷结果:
#[[3, 4, 1, 4, 1],
[6, 4, 4, 3, 1],
[6, 4, 3, 1, 3],
[4, 3, 5, 2, 6],
[2, 5, 3, 1, 6]]
#应用: 桌面游戏、赌博游戏的概率模拟
示例2：彩票号码生成
生成大乐透前区号码（1-35选5个）
result2 = randi(35, 5) # 从1-35中选5个不重复号码（注：这里可能会有重复，实际彩票应确保不重复）
print(f"彩票前区号码: {result2}")
#彩票前区号码:
#[[20, 20, 2, 19, 25],
[27, 7, 4, 18, 22],
[9, 17, 32, 24, 29],
[25, 2, 15, 1, 2],
[28, 7, 15, 20, 17]]
#应用: 随机选号，避免人为偏见
示例3：密码生成 - PIN码
生成4位数字PIN码
result3 = randi(9, 4) # 4位数字，每位0-9（注意：这里生成的是1-9，实际PIN码包含0）
print(f"随机PIN码: {result3}")
#随机PIN码:
#[[9, 7, 9, 7],
[9, 2, 5, 6],
[6, 9, 4, 5],
[3, 9, 5, 6]]
#应用: 账户安全、门禁系统
示例4：学生随机分组
将30名学生随机分成6组，每组5人
result4 = randi(30, 6, 5) # 6×5矩阵，每个元素代表学生编号
print(f"学生分组矩阵: {result4}")
#学生分组矩阵:
#[[2, 29, 26, 24, 28],
[7, 14, 16, 5, 30],
[15, 6, 26, 28, 9],
[10, 9, 14, 30, 7],
[5, 22, 15, 16, 4],
[15, 12, 4, 8, 19]]
#应用: 确保分组的随机性和公平性
示例5：产品质量抽检
从1000个产品中随机选择20个进行质量检查
result5 = randi(1000, 20) # 20个产品编号
print(f"抽检产品编号: {result5}")
#抽检产品编号:
#[[244, 659, 689, 544, 644, 467, 213, 918, 476, 31, 327, 81, 993, 538, 753, 232, 273, 832, 577, 61],
[608, 490, 831, 588, 461, 825, 393, 53, 830, 333, 496, 652, 328, 678, 487, 915, 68, 67, 152, 763],
[631, 144, 374, 90, 90, 277, 575, 875, 648, 545, 909, 611, 259, 686, 651, 445, 688, 583, 802, 113],
[114, 230, 578, 372, 125, 809, 550, 492, 273, 781, 536, 890, 173, 741, 177, 246, 234, 718, 4, 57],
[415, 780, 181, 143, 223, 863, 148, 434, 313, 755, 536, 726, 505, 243, 228, 218, 211, 445, 644, 717],
[151, 761, 840, 664, 421, 271, 737, 661, 574, 408, 78, 186, 236, 34, 557, 646, 152, 672, 699, 847],
[984, 204, 520, 58, 130, 336, 998, 556, 237, 39, 847, 570, 873, 17, 638, 253, 514, 242, 550, 479],
[83, 14, 226, 186, 70, 980, 609, 132, 122, 881, 905, 801, 772, 133, 689, 266, 471, 716, 478, 905],
[582, 742, 773, 118, 583, 368, 155, 727, 217, 729, 743, 243, 897, 931, 93, 280, 888, 805, 869, 472],
[582, 527, 746, 899, 519, 369, 631, 637, 693, 306, 964, 277, 817, 793, 752, 878, 563, 701, 138, 63],
[334, 826, 246, 220, 865, 963, 589, 677, 553, 741, 910, 479, 419, 920, 245, 284, 220, 649, 875, 243],
[23, 5, 907, 847, 882, 745, 490, 502, 230, 942, 842, 56, 514, 491, 773, 762, 94, 440, 98, 888],
[974, 399, 875, 499, 898, 120, 900, 151, 289, 856, 622, 413, 803, 295, 329, 307, 801, 529, 739, 289],
[971, 36, 155, 328, 625, 423, 692, 742, 721, 136, 734, 352, 355, 860, 976, 500, 916, 635, 800, 358],
[755, 683, 159, 550, 196, 450, 576, 19, 551, 827, 361, 566, 551, 351, 573, 994, 66, 863, 89, 301],
[474, 980, 330, 62, 687, 816, 194, 663, 984, 550, 696, 254, 884, 74, 396, 116, 123, 45, 311, 519],
[231, 77, 848, 661, 442, 373, 889, 107, 792, 212, 207, 673, 694, 797, 806, 191, 567, 540, 978, 666],
[287, 121, 995, 787, 874, 637, 715, 660, 605, 281, 323, 242, 940, 990, 137, 167, 785, 146, 969, 350],
[535, 587, 12, 697, 280, 204, 960, 15, 145, 533, 526, 198, 161, 707, 945, 897, 918, 369, 409, 867],
[441, 655, 99, 378, 519, 198, 756, 581, 418, 374, 311, 564, 222, 966, 241, 313, 485, 854, 83, 800]]
#应用: 质量控制，确保抽样代表性
示例6：随机迷宫生成
生成10×10迷宫的初始随机墙位置
result6 = randi(2, 10, 10) # 10×10矩阵，1表示墙，2表示通路
print(f"迷宫初始矩阵 (1=墙, 2=通路):\n{result6}")
#迷宫初始矩阵 (1=墙, 2=通路):
#[[1, 1, 1, 2, 1, 1, 1, 2, 1, 1],
[2, 1, 1, 1, 1, 1, 1, 1, 1, 2],
[2, 1, 2, 1, 2, 1, 1, 1, 2, 1],
[2, 2, 2, 1, 1, 1, 1, 1, 1, 1],
[2, 1, 2, 1, 2, 1, 1, 1, 1, 1],
[2, 2, 1, 1, 2, 2, 1, 1, 2, 1],
[1, 1, 2, 1, 2, 2, 2, 2, 2, 1],
[1, 1, 2, 2, 1, 1, 1, 2, 2, 1],
[2, 1, 2, 1, 2, 2, 2, 2, 1, 2],
[1, 2, 1, 2, 1, 2, 1, 2, 1, 1]]
#应用: 游戏开发、算法测试
示例7：实验条件随机分配
将40名参与者随机分配到4种实验条件
result7 = randi(4, 10) # 10×1矩阵，每个元素1-4代表实验条件
print(f"实验条件分配: {result7}")
#实验条件分配:
#[[4, 3, 2, 1, 2, 3, 4, 3, 2, 1],
[4, 1, 3, 1, 2, 3, 1, 3, 3, 1],
[3, 1, 2, 4, 2, 1, 1, 4, 4, 3],
[4, 1, 1, 4, 3, 2, 2, 3, 2, 2],
[3, 3, 2, 4, 3, 3, 2, 3, 2, 1],
[3, 3, 4, 3, 4, 4, 1, 4, 3, 3],
[3, 1, 1, 2, 4, 3, 4, 4, 1, 1],
[3, 3, 4, 4, 3, 3, 4, 1, 4, 2],
[2, 3, 3, 1, 2, 2, 1, 1, 4, 4],
[4, 2, 2, 1, 2, 1, 4, 2, 2, 2]]
#应用: 消除顺序效应，确保实验有效性
示例8：随机测试题目顺序
为50道题的试卷生成随机题目顺序
result8 = randi(50, 50) # 1-50的随机排列
print(f"题目随机顺序 (前10个): {result8[:10, :] if hasattr(result8, 'shape') else result8[:10]}")
#题目随机顺序 (前10个):
#[[5, 25, 49, 36, 50, 39, 41, 32, 5, 11, 41, 25, 17, 45, 28, 33, 15, 3, 15, 7, 40, 4, 25, 12, 28, 28, 10, 47, 36, 25, 24, 2, 39, 38, 44, 9, 45, 3, 22, 39, 26, 38, 27, 33, 33, 27, 6, 13, 35, 9],
[26, 12, 38, 13, 29, 6, 7, 45, 48, 29, 30, 45, 27, 13, 6, 41, 43, 14, 27, 41, 1, 36, 3, 41, 45, 38, 7, 20, 17, 3, 26, 21, 5, 35, 39, 26, 7, 37, 33, 50, 34, 11, 48, 32, 33, 40, 41, 26, 4, 5],
[47, 44, 20, 29, 32, 30, 41, 20, 15, 46, 21, 5, 38, 28, 29, 11, 14, 40, 49, 17, 38, 44, 2, 2, 49, 15, 30, 9, 1, 32, 10, 15, 29, 2, 30, 27, 2, 22, 21, 25, 5, 3, 15, 23, 35, 48, 39, 43, 4, 50],
[6, 15, 24, 33, 22, 35, 8, 50, 50, 2, 17, 36, 26, 30, 16, 16, 2, 18, 17, 1, 13, 49, 45, 16, 47, 3, 35, 17, 20, 21, 15, 3, 13, 7, 42, 49, 43, 9, 20, 49, 17, 21, 13, 34, 10, 22, 3, 11, 42, 21],
[5, 14, 17, 14, 34, 12, 5, 10, 37, 21, 18, 22, 46, 16, 19, 37, 21, 6, 26, 8, 16, 1, 26, 42, 4, 14, 37, 4, 49, 2, 16, 18, 14, 28, 37, 9, 29, 25, 17, 23, 15, 50, 13, 34, 18, 8, 45, 14, 44, 23],
[32, 50, 13, 48, 25, 50, 4, 9, 49, 38, 18, 11, 18, 16, 8, 15, 50, 46, 45, 2, 8, 22, 33, 41, 26, 11, 6, 22, 16, 3, 18, 50, 4, 27, 40, 22, 40, 16, 18, 29, 14, 23, 10, 35, 31, 46, 17, 40, 41, 9],
[25, 7, 17, 23, 1, 46, 26, 41, 25, 32, 24, 19, 42, 40, 19, 34, 27, 35, 32, 2, 24, 49, 42, 31, 12, 7, 47, 38, 31, 13, 22, 43, 11, 27, 5, 39, 34, 11, 42, 41, 10, 23, 19, 2, 29, 33, 39, 12, 47, 10],
[25, 27, 23, 8, 47, 39, 35, 23, 29, 12, 7, 35, 47, 40, 48, 23, 26, 15, 7, 48, 35, 11, 28, 49, 30, 40, 24, 48, 28, 29, 10, 27, 6, 47, 19, 39, 2, 21, 31, 33, 25, 43, 47, 13, 28, 2, 4, 10, 17, 14],
[23, 13, 36, 46, 47, 31, 7, 37, 19, 15, 47, 12, 48, 1, 37, 19, 30, 19, 22, 20, 11, 31, 23, 46, 6, 18, 48, 27, 48, 4, 18, 35, 29, 16, 49, 41, 32, 17, 17, 28, 17, 14, 42, 18, 30, 30, 5, 37, 36, 7],
[9, 47, 19, 29, 6, 12, 43, 33, 33, 3, 25, 49, 33, 13, 22, 45, 18, 25, 50, 8, 40, 13, 8, 44, 27, 12, 14, 21, 8, 34, 32, 33, 18, 30, 1, 3, 30, 37, 45, 20, 33, 21, 20, 38, 4, 4, 42, 49, 34, 30]]
#应用: 防止作弊，公平考试
示例9：游戏敌人属性生成
为游戏生成10个敌人的随机等级（1-10级）
result9 = randi(10, 10) # 10个敌人的等级
print(f"敌人等级: {result9}")
#敌人等级:
#[[9, 8, 9, 6, 4, 9, 2, 1, 5, 4],
[10, 2, 9, 10, 10, 6, 10, 7, 6, 10],
[2, 4, 7, 8, 4, 10, 7, 10, 10, 6],
[2, 4, 8, 5, 8, 6, 5, 2, 6, 5],
[6, 3, 4, 4, 8, 6, 3, 10, 9, 8],
[8, 7, 5, 6, 8, 4, 10, 6, 1, 9],
[9, 7, 5, 1, 2, 2, 6, 2, 8, 8],
[2, 1, 7, 9, 6, 1, 9, 4, 7, 1],
[8, 3, 5, 1, 10, 7, 6, 4, 1, 1],
[7, 3, 9, 3, 2, 4, 2, 10, 2, 5]]
#应用: 游戏平衡，增加游戏多样性
示例10：会议座位随机安排
为20人的会议生成5×4座位安排
result10 = randi(20, 5, 4) # 5排4列，每个座位分配1-20的编号
print(f"座位安排矩阵:\n{result10}")
#座位安排矩阵:
#[[9, 1, 13, 2],
[20, 14, 19, 11],
[14, 12, 15, 16],
[5, 17, 12, 2],
[17, 3, 14, 19]]
#应用: 活动策划，促进随机交流
示例11：模拟随机投票
模拟100人对3个候选人的随机投票
result11 = randi(3, 100) # 100个投票，每个投票选择候选人1-3
print(f"投票结果样本: {result11[:10] if hasattr(result11, 'shape') else result11[:10]}")
#投票结果样本: [2, 1, 2, 3, 3, 2, 2, 2, 1, 3]
#应用: 选举预测、政治学研究
示例12：随机颜色生成（RGB分量）
生成5个随机颜色的RGB分量（0-255）
result12 = randi(255, 5, 3) # 5个颜色，每个颜色有3个RGB分量
print(f"随机RGB颜色矩阵:\n{result12}")
#随机RGB颜色矩阵:
#[[135, 153, 154],
[66, 220, 21],
[221, 53, 252],
[211, 86, 53],
[211, 210, 117]]
#应用: 图形设计、数据可视化
示例13：随机音乐播放列表
从50首歌曲中随机选择10首播放
result13 = randi(50, 10) # 10首随机选择的歌曲编号
print(f"随机播放列表: {result13}")
#随机播放列表:
#[[15, 50, 37, 15, 42, 27, 13, 50, 22, 24],
[15, 36, 45, 5, 34, 37, 23, 5, 24, 22],
[45, 23, 6, 32, 21, 46, 28, 17, 25, 48],
[19, 21, 36, 6, 33, 41, 10, 24, 16, 16],
[39, 38, 43, 39, 45, 27, 30, 48, 9, 27],
[4, 14, 31, 11, 28, 49, 40, 8, 7, 1],
[14, 32, 2, 41, 28, 16, 23, 38, 35, 7],
[41, 5, 42, 45, 24, 48, 39, 48, 22, 2],
[21, 16, 33, 42, 13, 18, 5, 6, 1, 30],
[1, 3, 35, 38, 10, 19, 11, 17, 27, 44]]
#应用: 音乐播放器，避免重复模式
示例14：模拟随机网络延迟
模拟100个数据包的随机延迟（1-100ms）
result14 = randi(100, 100) # 100个延迟时间(ms)
print(f"网络延迟样本 (前10个): {result14[:10] if hasattr(result14, 'shape') else result14[:10]}")
#网络延迟样本 (前10个): [14, 4, 79, 46, 9, 8, 26, 34, 20, 33]
#应用: 网络性能测试、QoS评估
示例15：随机优惠券代码生成
生成10个6位数字的优惠券代码
result15 = randi(9, 10, 6) # 10个代码，每个6位数字（1-9）
print(f"优惠券代码矩阵:\n{result15}")
print("应用: 营销活动、客户激励")
#优惠券代码矩阵:
#[[5, 7, 2, 1, 2, 8],
[8, 7, 6, 2, 9, 4],
[4, 6, 6, 3, 8, 3],
[5, 6, 7, 2, 2, 1],
[7, 6, 9, 2, 1, 3],
[5, 1, 3, 3, 2, 8],
[1, 9, 1, 4, 6, 9],
[5, 8, 1, 1, 1, 4],
[1, 3, 7, 9, 1, 5],
[6, 5, 2, 1, 7, 7]]
#应用: 营销活动、客户激励
示例16：随机迷宫难度设置
为游戏生成5个关卡的随机难度（1-5级）
result17 = randi(5, 5) # 5个关卡的难度级别
print(f"关卡难度设置: {result17}")
#关卡难度设置:
#[[5, 4, 5, 1, 1],
[2, 4, 3, 3, 1],
[5, 2, 1, 5, 4],
[3, 3, 3, 1, 3],
[4, 3, 3, 3, 1]]
#应用: 游戏设计，增加重玩价值
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def random_imax_generation(input_str):
"""
根据输入字符串生成均匀分布的伪随机整数，对标MATLAB的randi函数
参数：
input_str (str): 表示生成范围和维度的字符串，支持格式：
- "max_val" : 生成单个1~max_val的整数
- "(max_val, size)" : 生成size×size的矩阵
- "(max_val, rows, cols)" : 生成rows×cols的矩阵
返回：
sympy.Matrix/int/str: 生成的整数或矩阵，错误时返回字符串
示例：
>>> random_imax_generation("5") # 生成1~5的随机整数
>>> random_imax_generation("(5, 2)") # 2×2矩阵，元素1~5
>>> random_imax_generation("(5, 2, 3)") # 2×3矩阵，元素1~5
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
# 处理元组参数
if len(expr) == 3:
# 格式：(max_val, rows, cols)
max_val = int(expr[0])
rows = int(expr[1])
cols = int(expr[2])
result = sp.Matrix(np.random.randint(1, max_val + 1, (rows, cols)))
elif len(expr) == 2:
# 格式：(max_val, size)
max_val = int(expr[0])
size = int(expr[1])
result = sp.Matrix(np.random.randint(1, max_val + 1, (size, size)))
else:
error = True # 无效元组长度
elif expr.is_integer:
# 格式：单个整数max_val
max_val = int(expr)
result = np.random.randint(1, max_val + 1)
else:
error = True # 无效输入类型
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：生成单个随机整数
print("示例1:", random_imax_generation("5"))
# 1
# 示例2：生成3×3矩阵
print("\n示例2:", random_imax_generation("(5, 3)"))
# Matrix([[5, 5, 3],
# [3, 5, 1],
# [2, 2, 2]])
# 示例3：生成2×4矩阵
print("\n示例3:", random_imax_generation("(10, 2, 4)"))
# Matrix([[3, 1, 3, 6],
# [5, 10, 8, 3]])
随机正态矩阵
X = randn(n) 返回一个由正态分布的随机数组成的 1×n 矩阵
X = randn(m,n) 返回一个由正态分布的随机数组成的 m×n 矩阵.
m, n — 方阵的大小,整数值
示例1：金融资产收益率模拟
模拟10只股票100天的日收益率，假设平均日收益率为0.1%，标准差为2%
result1 = randn(10, 100, 0.001, 0.02) # 10只股票，100天，均值0.001，标准差0.02
print(f"股票收益率矩阵: {result1.rows}×{result1.cols}")
#股票收益率矩阵: 10×100
#应用: 投资组合风险分析，蒙特卡洛模拟
示例2：传感器测量误差模拟
模拟5个传感器对同一物理量的100次测量，假设真值为10，测量误差标准差为0.1
result2 = randn(5, 100, 10, 0.1) # 5个传感器，100次测量，均值10，标准差0.1
print(f"传感器测量矩阵: {result2.rows}×{result2.cols}")
#传感器测量矩阵: 5×100
#应用: 评估传感器精度，误差分析
示例3：神经网络噪声注入
为图像数据添加高斯噪声，增强模型鲁棒性
result3 = randn(64, 64, 0, 0.1) # 64×64图像，添加均值为0，标准差0.1的噪声
print(f"图像噪声矩阵: {result3.rows}×{result3.cols}")
#图像噪声矩阵: 64×64
#应用: 防止过拟合，提高模型泛化能力
示例4：气候温度变化模拟
模拟10个城市365天的温度变化，基准温度20°C，日变化标准差3°C
result4 = randn(10, 365, 20, 3) # 10个城市，365天，均值20，标准差3
print(f"气温变化矩阵: {result4.rows}×{result4.cols}")
#气温变化矩阵: 10×365
#应用: 气候变化研究，能源需求预测
示例5：产品质量控制
模拟生产线上100个产品的5个关键尺寸，目标尺寸50mm，公差±0.5mm(3σ)
result5 = randn(100, 5, 50, 0.5/3) # 100个产品，5个尺寸，均值50，标准差0.1667
print(f"产品尺寸矩阵: {result5.rows}×{result5.cols}")
#产品尺寸矩阵: 100×5
#应用: 过程能力分析，质量控制
示例6：信号处理 - 噪声添加
模拟通信信道中的加性高斯白噪声(AWGN)
result6 = randn(1000, 1, 0, 0.5) # 1000个信号样本，均值为0，标准差0.5的噪声
print(f"信道噪声向量: {result6.rows}×{result6.cols}")
#信道噪声向量: 1000×1
#应用: 通信系统性能评估，误码率分析
示例7：经济学 - 宏观经济指标模拟
模拟5个宏观经济指标在50个时间点的随机波动
result7 = randn(5, 50, 0, 0.02) # 5个指标，50个时间点，均值0，标准差0.02
print(f"经济指标波动矩阵: {result7.rows}×{result7.cols}")
#经济指标波动矩阵: 5×50
#应用: 经济预测，政策影响分析
示例8：医学研究 - 生物标志物变异
模拟100名患者的3种生物标志物测量值，考虑个体差异和测量误差
result8 = randn(100, 3, 0, 1) # 100名患者，3种标志物，标准正态分布
print(f"生物标志物矩阵: {result8.rows}×{result8.cols}")
#生物标志物矩阵: 100×3
#应用: 临床试验数据分析，疾病诊断模型
示例9：计算机图形学 - 自然纹理生成
使用高斯随机场生成自然纹理的基础
result9 = randn(256, 256) # 256×256标准正态随机场
print(f"随机场矩阵: {result9.rows}×{result9.cols}")
#随机场矩阵: 256×256
#应用: 云层、地形、大理石等自然纹理生成
示例10：心理学 - 认知测试得分模拟
模拟1000名受试者在5项认知测试中的得分，平均分100，标准差15
result10 = randn(1000, 5, 100, 15) # 1000名受试者，5项测试，均值100，标准差15
print(f"认知测试得分矩阵: {result10.rows}×{result10.cols}")
#认知测试得分矩阵: 1000×5
#应用: 智力测验标准化，常模建立
示例11：物理学 - 粒子位置不确定性
模拟量子粒子位置的统计波动
result11 = randn(100, 3, 0, 0.1) # 100个粒子，3D位置，均值0，标准差0.1
print(f"粒子位置波动矩阵: {result11.rows}×{result11.cols}")
#粒子位置波动矩阵: 100×3
#应用: 量子力学模拟，不确定性原理演示
示例12：风险管理 - 极端事件模拟
使用厚尾分布模拟极端风险事件（通过正态分布近似）
result12 = randn(10000, 1, 0, 1) # 10000个场景，标准正态分布
extreme_events = np.sum(np.abs(result12) > 2) # 统计超过2σ的事件
print(f"极端事件数量(|值|>2σ): {extreme_events}/10000")
#极端事件数量(|值|>2σ): 454/10000
#应用: 风险价值(VaR)计算，压力测试
示例13：声学 - 环境噪声建模
模拟不同频率的环境噪声强度
result13 = randn(20, 100, 60, 10) # 20个频率点，100个时间样本，均值60dB，标准差10dB
print(f"噪声强度矩阵: {result13.rows}×{result13.cols}")
#噪声强度矩阵: 20×100
#应用: 声学设计，噪声控制
示例14：机器学习 - 潜在空间采样
在潜在空间中采样生成新数据点
result14 = randn(10, 32) # 10个样本，32维潜在空间，标准正态分布
print(f"潜在空间样本: {result14.rows}×{result14.cols}")
#潜在空间样本: 10×32
#应用: 生成模型，数据合成
示例15：工程学 - 结构载荷变化
模拟建筑物承受的随机风载荷
result15 = randn(100, 1, 0, 50) # 100个时间点，均值为0，标准差50kN的载荷变化
print(f"风载荷变化向量: {result15.rows}×{result15.cols}")
#风载荷变化向量: 100×1
#应用: 结构安全分析，疲劳寿命预测
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def matrix_randn(input_str):
"""
根据输入字符串生成正态分布的随机矩阵，对标MATLAB的randn函数，支持扩展参数。
参数：
input_str (str): 表示矩阵维度和分布参数的字符串。支持以下格式：
- "n"：生成n×n的标准正态分布矩阵（均值0，标准差1）
- "(rows, cols)"：生成rows×cols的标准正态分布矩阵
- "(rows, cols, mu)"：生成均值为mu，标准差为1的正态分布矩阵
- "(rows, cols, mu, sigma)"：生成均值为mu，标准差为sigma的正态分布矩阵
返回：
sympy.Matrix or str: 生成的随机矩阵，若输入无效则返回错误信息。
示例：
>>> matrix_randn("3") # 3x3标准正态分布矩阵
>>> matrix_randn("(2,3)") # 2x3标准正态分布矩阵
>>> matrix_randn("(2,2,5)") # 2x2矩阵，均值5，标准差1
>>> matrix_randn("(3,2,5,2)") # 3x2矩阵，均值5，标准差2
"""
try:
expr = sp.sympify(input_str) # 解析输入字符串为SymPy表达式
error = False
result = None
if isinstance(expr, tuple):
# 处理元组参数
if len(expr) == 2:
# 格式：(rows, cols)，生成标准正态分布
rows = int(expr[0])
cols = int(expr[1])
random_matrix = np.random.randn(rows, cols)
result = sp.Matrix(random_matrix)
elif len(expr) == 3:
# 格式：(rows, cols, mu)，生成N(mu,1)
rows = int(expr[0])
cols = int(expr[1])
mu = float(expr[2])
random_matrix = np.random.randn(rows, cols) + mu
result = sp.Matrix(random_matrix)
elif len(expr) == 4:
# 格式：(rows, cols, mu, sigma)，生成N(mu, sigma^2)
rows = int(expr[0])
cols = int(expr[1])
mu = float(expr[2])
sigma = float(expr[3])
random_matrix = sigma * np.random.randn(rows, cols) + mu
result = sp.Matrix(random_matrix)
else:
error = True # 不支持的参数长度
elif expr.is_integer:
# 格式：整数n，生成n×n标准正态分布矩阵
n = int(expr)
random_matrix = np.random.randn(n, n)
result = sp.Matrix(random_matrix)
else:
error = True # 无效的输入类型
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：3x3标准正态分布矩阵
print("示例1:", matrix_randn("3"))
# Matrix([[1.39013708918388, -0.324752480460865, -0.190469074314680],
# [0.924750724511891, 1.15842217229967, -0.111156310307857],
# [-0.142219586288319, 0.955105925465763, -0.331834544891866]])
# 示例2：2x3标准正态分布矩阵
print("示例2:", matrix_randn("(2,3)"))
# Matrix([[-1.03503301334510, 0.187320578973013, 0.00353779280429042],
# [1.73695585011864, 0.966429881059110, -1.43696858006492]])
# 示例3：2x2矩阵，均值5，标准差1
print("示例3:", matrix_randn("(2,2,5)"))
# Matrix([[6.15546425528767, 6.43937825874499],
# [5.58857185229244, 6.65723732023311]])
# 示例4：3x2矩阵，均值5，标准差2
print("示例4:", matrix_randn("(3,2,5,2)"))
# Matrix([[4.48575907428833, 5.43818803588287],
# [6.50850096480664, 6.91253704333354],
# [2.45699214062053, 3.87080672903311]])
整数的随机排列
p = randperm(n) 返回行向量,其中包含从1到n没有重复元素的整数随机排列.
p = randperm(n,k) 返回行向量,其中包含在1到n之间随机选择的k个唯一整数.
n — 采样间隔中的整数数量,正整数值
k — 所选整数的数量,正整数值
示例1：抽奖系统
从100名参与者中随机抽取3名获奖者
result1 = randperm(100, 3)
print(f"中奖者编号: {result1}")
#中奖者编号: [64 38 39]
#应用: 公平抽奖，确保每个参与者有同等机会
示例2：考试座位安排
为30名考生随机分配座位号
result2 = randperm(30)
print(f"考生座位安排: {result2}")
#考生座位安排: [14 24 13 17 9 15 12 18 30 3 8 6 1 21 27 26 22 29 11 5 16 23 20 28 10 19 7 25 2 4]
#应用: 防止作弊，确保座位分配的随机性
示例3：产品质检抽样
从生产批次的500个产品中随机抽取20个进行质检
result3 = randperm(500, 20)
print(f"抽检产品编号: {result3}")
#抽检产品编号: [263 209 41 362 136 50 225 302 124 181 221 175 327 500 85 42 246 247 58 489]
#应用: 随机抽样确保质量评估的代表性
示例4：音乐播放列表随机排序
对包含15首歌曲的播放列表进行随机排序
result4 = randperm(15)
print(f"随机播放顺序: {result4}")
#随机播放顺序: [ 3 4 7 8 2 12 15 14 5 9 13 1 6 11 10]
#应用: 避免重复听歌顺序，增加新鲜感
示例5：实验分组
将24只实验动物随机分为3组，每组8只
group_size = 8
total_animals = 24
result5 = randperm(total_animals, group_size)
print(f"第一组动物编号: {result5}")
#第一组动物编号: [23 6 1 20 11 14 4 5]
#应用: 随机分组消除选择偏差，确保实验有效性
示例6：密码生成
从26个字母中随机选择8个字符作为密码基础
result6 = randperm(26, 8)
转换为字母
letters = [chr(96 + num) for num in result6] # a=97 in ASCII
print(f"随机字符序列: {''.join(letters)}")
#随机字符序列: gsrtlkao
#应用: 增强密码的随机性和安全性
示例7：旅游路线规划
从10个必游景点中随机选择访问顺序
result7 = randperm(10)
print(f"景点访问顺序: {result7}")
#景点访问顺序: [ 9 5 1 7 6 2 8 3 4 10]
#应用: 为每次旅行创造独特的体验
示例8：研究调查对象选择
从1000人的总体中随机选择50人进行调查
result8 = randperm(1000, 50)
print(f"调查对象编号: {result8}")
#调查对象编号: [378 286 263 322 535 205 320 874 405 226 986 690 187 570 329 978 654 382
767 245 214 524 12 56 728 738 668 984 365 518 419 699 100 589 138 955
321 10 197 202 818 429 277 925 963 280 935 795 787 736]
#应用: 确保样本代表性，减少抽样偏差
示例9：比赛对阵抽签
为8支队伍生成随机比赛顺序
result9 = randperm(8)
print(f"队伍比赛顺序: {result9}")
#队伍比赛顺序: [4 5 2 3 7 6 8 1]
#应用: 公平抽签，避免人为操纵
示例10：记忆测试材料
为记忆测试准备20个单词的随机顺序
result10 = randperm(20)
print(f"测试单词顺序: {result10}")
#测试单词顺序: [ 7 16 15 1 18 20 14 3 6 19 13 4 8 10 17 2 9 12 11 5]
#应用: 消除顺序效应，提高实验准确性
示例11：机器学习训练集划分
从1000个样本中随机选择200个作为验证集
result11 = randperm(1000, 200)
print(f"验证集样本索引: {result11}")
#验证集样本索引:
#[568 738 65 902 429 70 332 460 661 219 442 854 148 799 860 246 12 116
885 638 189 871 654 683 379 101 549 934 693 92 163 210 89 787 602 714
759 80 361 615 212 330 952 596 426 688 646 584 679 417 309 913 311 181
238 630 585 182 479 57 552 886 443 345 667 505 524 851 411 336 310 501
663 273 300 565 130 908 59 282 516 225 405 242 390 898 472 998 757 297
91 226 784 286 342 500 690 526 551 503 402 792 367 344 40 8 509 468
875 407 532 974 255 836 237 227 906 589 844 469 821 925 222 365 119 74
697 969 781 984 732 785 901 748 113 359 478 211 169 724 822 884 957 139
999 677 520 448 729 985 437 401 343 964 369 543 514 938 709 768 241 700
877 730 295 250 917 645 608 35 706 275 287 916 44 81 750 894 518 93
879 881 995 340 597 22 258 328 284 682 556 648 127 816 456 659 366 16
180 575]
#应用: 确保训练集和验证集的随机划分
示例12：任务分配系统
将15项任务随机分配给团队成员
result12 = randperm(15)
print(f"任务执行顺序: {result12}")
#任务执行顺序: [10 15 2 6 4 5 9 12 3 8 14 1 7 11 13]
#应用: 公平分配工作量，避免偏见
示例13：游戏卡牌洗牌
模拟52张扑克牌的随机洗牌
result13 = randperm(52)
print(f"洗牌后牌序 (前10张): {result13[:10]}")
#洗牌后牌序 (前10张): [24 5 42 35 20 46 38 17 22 40]
#应用: 确保游戏公平性和不可预测性
示例14：A/B测试用户分组
将1000名用户随机分为A组（500人）进行新功能测试
result14 = randperm(1000, 500)
print(f"A测试组用户ID: {result14}")
#A测试组用户ID:
#[695 555 791 861 572 485 449 279 951 208 582 352 498 998 440 901 875 652
543 465 851 740 347 825 513 143 148 378 266 874 815 291 470 86 823 395
429 396 697 238 299 356 436 950 432 698 183 386 586 992 991 876 515 729
833 737 576 150 925 460 629 994 845 348 309 308 501 593 368 707 569 581
810 603 398 542 400 607 177 987 600 73 651 101 772 968 888 562 648 381
516 557 50 140 260 172 366 202 312 646 958 133 863 798 196 529 623 360
464 40 158 128 903 334 448 12 620 975 511 187 688 397 656 848 65 47
799 153 441 855 566 413 814 696 664 80 81 550 539 217 173 545 923 836
190 55 232 989 61 154 211 207 455 212 467 130 142 828 487 567 857 4
461 340 16 911 38 405 175 619 114 282 547 929 274 988 533 659 763 882
816 248 370 293 638 359 246 748 840 837 297 957 632 505 231 976 273 784
588 298 23 56 43 790 188 530 392 647 986 343 480 416 918 354 666 762
362 650 92 587 741 503 426 166 711 996 372 618 935 410 723 843 920 694
165 391 727 390 344 993 703 526 801 136 327 306 644 453 277 639 49 916
491 300 625 561 117 709 337 259 219 269 393 885 346 303 123 53 271 752
564 883 517 917 51 509 18 728 606 84 904 847 371 261 751 952 302 301
700 605 240 583 228 195 456 649 197 800 258 222 924 686 224 964 827 72
263 912 502 341 598 113 658 544 919 54 267 201 349 29 357 523 459 813
473 193 144 683 559 388 895 590 418 317 764 864 858 251 147 233 849 859
242 914 974 549 350 787 189 163 732 926 103 726 884 850 877 909 286 963
199 444 624 24 614 96 670 107 693 77 831 463 965 687 985 74 156 522
424 902 110 809 285 227 706 335 927 541 466 773 494 510 584 160 565 364
206 891 46 468 514 492 948 268 425 169 941 921 486 476 776 783 42 100
93 374 22 445 677 712 770 45 182 88 898 955 626 296 387 868 170 621
307 609 129 811 146 718 860 792 330 571 508 48 375 36 613 839 954 64
234 602 668 345 982 181 90 490 252 599 739 255 447 797 119 108 553 682
380 669 20 70 966 830 769 451 262 94 691 141 305 675 940 627 780 310
475 527 176 151 915 676 725 63 126 672 570 785 892 534]
#应用: 随机分组确保测试结果的可靠性
示例15：教学材料随机展示
从题库的50道题中随机选择10道展示
result15 = randperm(50, 10)
print(f"随机题目顺序: {result15}")
#随机题目顺序: [ 9 11 43 24 34 42 20 50 46 5]
#应用: 防止学生记忆答案顺序，提高学习效果
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def random_permutation_generation(input_str):
"""
生成整数的随机排列，对标MATLAB的randperm函数
参数：
input_str (str): 表示排列参数的字符串，支持格式：
- "n" : 生成1~n的随机排列
- "(n, k)" : 从1~n中随机选择k个不重复整数
返回：
numpy.ndarray/sympy.Matrix: 生成的随机排列数组或矩阵
str: 输入错误时的提示信息
示例：
>>> random_permutation_generation("5") # 生成1~5的随机排列
>>> random_permutation_generation("(10,3)") # 从1~10选3个随机数
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理元组参数 (n, k)
if isinstance(expr, tuple):
if len(expr) == 2:
# 验证参数是否为有效数字
if expr[0].is_integer and expr[1].is_integer:
n = int(expr[0])
k = int(expr[1])
# 参数有效性检查
if n <= 0 or k <= 0:
raise ValueError("参数必须为正整数")
if k > n:
raise ValueError("k不能大于n")
# 生成随机排列（MATLAB的randperm(n,k)行为）
result = np.random.choice(np.arange(1, n + 1), k, replace=False)
else:
error = True
else:
error = True
# 处理单个整数参数 n
elif expr.is_integer:
n = int(expr)
if n <= 0:
raise ValueError("参数必须为正整数")
# 生成完整随机排列（MATLAB的randperm(n)行为）
result = np.random.permutation(np.arange(1, n + 1))
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 正常用例
print("=== 正常用例测试 ===")
print("1. 基础排列:", random_permutation_generation("5"))
# [3 5 2 4 1]
print("\n2. 子集选择:", random_permutation_generation("(10,3)"))
# [ 9 10 2]
# 边界用例
print("\n=== 边界用例测试 ===")
print("3. 最小参数:", random_permutation_generation("(1,1)"))
# [1]
print("4. 最大选择:", random_permutation_generation("(5,5)"))
# [4 5 3 1 2]
随机列表或矩阵
Random(a,b,size) -- 生成指定范围内的均匀分布随机数，支持标量/向量/矩阵/高维数组
a,b,size -- 标量
示例1. 物理模拟应用
粒子物理模拟 - 生成随机初始位置和速度
生成100个粒子的初始位置 (0-10米范围内)
particle_positions = Random(0,10,100,3)
print(f"粒子位置矩阵形状: {particle_positions.shape}")
#粒子位置矩阵形状: (100, 3)
生成粒子初始速度 (-5到5 m/s)
particle_velocities = Random(-5,5,100,3)
print(f"粒子速度矩阵形状: {particle_velocities.shape}")
#粒子速度矩阵形状: (100, 3)
温度场模拟 - 2D温度分布
temperature_field = Random(20,35,50,50)
print(f"温度场矩阵形状: {temperature_field.shape}")
#温度场矩阵形状: (50, 50)
示例2. 金融风险评估
生成随机股票收益率 (-0.1到0.1，即±10%)
daily_returns = Random(-0.1,0.1,252)
print(f"252个交易日的收益率序列: {len(daily_returns)}个数据点")
#252个交易日的收益率序列: 252个数据点
投资组合相关性矩阵模拟
correlation_matrix = Random(0.2,0.9,10,10)
print(f"10种资产的相关性矩阵: {correlation_matrix.shape}")
#10种资产的相关性矩阵: (10, 10)
蒙特卡洛模拟 - 1000次情景测试
monte_carlo_scenarios = Random(0.8,1.2,1000,5)
print(f"蒙特卡洛模拟: {monte_carlo_scenarios.shape}")
#蒙特卡洛模拟: (1000, 5)
示例3. 机器学习数据集生成
生成分类问题的训练数据
X_train = Random(0,10,1000,5) # 1000个样本，5个特征
y_train = Random(0,2,1000) # 二分类标签
print(f"训练数据: X={X_train.shape}, y={len(y_train)}")
#训练数据: X=(1000, 5), y=1000
生成测试数据
X_test = Random(0,10,200,5)
print(f"测试数据: {X_test.shape}")
#测试数据: (200, 5)
神经网络权重初始化
neural_weights = Random(-0.5,0.5,128,64)
print(f"神经网络权重矩阵: {neural_weights.shape}")
#神经网络权重矩阵: (128, 64)
示例4. 工程优化设计
结构参数随机采样
beam_lengths = Random(1,10,50) # 梁长度
material_strengths = Random(200,500,50) # 材料强度
print(f"结构参数采样: {len(beam_lengths)}个设计方案")
#结构参数采样: 50个设计方案
3D打印参数优化
print_parameters = Random(0.1,0.3,20,4) # 20组打印参数
print(f"3D打印参数矩阵: {print_parameters.shape}")
#3D打印参数矩阵: (20, 4)
示例5. 游戏开发应用
生成随机地形高度图
terrain_heightmap = Random(0,100,256,256)
print(f"地形高度图: {terrain_heightmap.shape}")
#地形高度图: (256, 256)
NPC随机属性生成
npc_attributes = Random(50,100,100,6) # 100个NPC的6种属性
print(f"NPC属性矩阵: {npc_attributes.shape}")
#NPC属性矩阵: (100, 6)
道具掉落概率
drop_rates = Random(0.01,0.1,20)
print(f"道具掉落概率: {len(drop_rates)}种道具")
#道具掉落概率: 20种道具
示例6. 科学研究应用
实验数据模拟 - 3重复3因素
experimental_data = Random(10,50,3,3,3)
print(f"实验数据立方: {experimental_data.shape}")
#实验数据立方: (3, 3, 3)
传感器网络读数
sensor_readings = Random(20,30,24,10) # 24小时，10个传感器
print(f"传感器数据: {sensor_readings.shape}")
#传感器数据: (24, 10)
分子动力学模拟初始状态
molecular_positions = Random(0,1,1000,3,10) # 1000个分子，10个时间步
print(f"分子位置张量: {molecular_positions.shape}")
#分子位置张量: (1000, 3, 10)
示例7. 统计分析应用
生成多元正态分布的随机样本（通过均匀分布近似）
multivariate_data = Random(-2,2,500,4)
print(f"多元数据: {multivariate_data.shape}")
#多元数据: (500, 4)
bootstrap重采样索引
bootstrap_samples = Random(0,1000,100,1000)
print(f"Bootstrap样本矩阵: {bootstrap_samples.shape}")
#Bootstrap样本矩阵: (100, 1000)
随机对照试验分组
treatment_effects = Random(0.5,2.5,30)
control_effects = Random(0.1,0.8,30)
print(f"实验组效果: {len(treatment_effects)}个样本")
#实验组效果: 30个样本
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def random_generation(input_str):
"""
生成指定范围内的均匀分布随机数，支持标量/向量/矩阵/高维数组
参数：
input_str (str): 表示范围和维度的字符串，格式为：
- "(low, high)" : 生成单个 [low, high) 之间的随机数
- "(low, high, size)" : 生成一维数组
- "(low, high, rows, cols)" : 生成二维矩阵(sympy.Matrix)
- "(low, high, d1, d2,...,dn)" : 生成n维numpy数组
返回：
[float/sp.Matrix/np.ndarray]: 生成的随机数结构
示例：
>>> random_generation("(0,1)") # 单个0~1的随机数
>>> random_generation("(2,5,3)") # 3个元素的1D数组
>>> random_generation("(0.5,1.5,2,2)") # 2x2的矩阵
>>> random_generation("(10,20,3,2,4)") # 3x2x4的3D数组
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
# 参数必须包含至少low和high两个参数
if len(expr) >= 2:
low = float(expr[0]) # 转换为浮点数确保精度
high = float(expr[1])
# 处理不同维度情况
if len(expr) == 2:
# 格式1：生成单个随机数
result = np.random.uniform(low, high)
else:
# 转换维度参数为整数元组
shape = tuple(int(d) for d in expr[2:])
# 生成numpy数组
arr = np.random.uniform(low, high, shape)
# 二维矩阵转sympy类型以保持一致性
if len(shape) == 2:
result = sp.Matrix(arr)
else:
result = arr
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__":
# 标量测试
print("标量测试:", random_generation("(0,1)"))
# 0.7152685221222377
# 一维数组测试
print("\n一维数组:", random_generation("(2,5,3)"))
# [2.41453067 3.70404521 4.03557003]
# 二维矩阵测试
print("\n二维矩阵:", random_generation("(0.5,1.5,2,2)"))
# Matrix([[0.986714064531144, 1.14702670344779],
# [0.860102519975358, 0.934994265190028]])
# 三维数组测试
print("\n三维数组:", random_generation("(10,20,3,2,4)"))
# [[[18.00333876 10.34994515 16.50753459 17.67988035]
# [16.09993668 18.59615135 16.9802518 19.65361739]]
#
# [[11.24132715 13.53054531 16.41626786 18.35072047]
# [19.62384891 13.28545208 18.91747723 17.68413307]]
#
# [[15.77014999 19.95114792 18.64022946 13.05333766]
# [16.3620224 10.64245557 11.13328264 17.22049697]]]
矩阵的秩
k = rank(A) 返回矩阵 A 的秩.
A — 输入矩阵,矩阵
示例1. 线性方程组解的存在性分析
有唯一解的方程组
A1 = [[1,2,3],[4,5,6],[7,8,10]]
b1 = [[1],[2],[3]]
print(f"系数矩阵秩: {rank(A1)}")
print(f"增广矩阵秩: {rank([[1,2,3,1],[4,5,6,2],[7,8,10,3]])}")
#系数矩阵秩: 3
#增广矩阵秩: 3
#结论：有唯一解 ✓
无解的方程组
A2 = [[1,2],[2,4]]
b2 = [[1],[3]]
print(f"\n系数矩阵秩: {rank(A2)}")
print(f"增广矩阵秩: {rank([[1,2,1],[2,4,3]])}")
#系数矩阵秩: 1
#增广矩阵秩: 2
#结论：无解 ×
无穷多解的方程组
A3 = [[1,2,3],[2,4,6],[3,6,9]]
print(f"\n系数矩阵秩: {rank(A3)}")
print(f"增广矩阵秩: {rank([[1,2,3,1],[2,4,6,2],[3,6,9,3]])}")
#系数矩阵秩: 1
#增广矩阵秩: 1
#结论：无穷多解 ∞
示例2. 控制系统分析
状态空间模型的可控性矩阵秩检验
A = [[0,1],[-2,-3]] # 系统矩阵
B = [[0],[1]] # 输入矩阵
构建可控性矩阵 [B, AB, A²B, ...]
controllability_matrix = [[0,1],[1,-3]]
controllability_rank = rank(controllability_matrix)
print(f"可控性矩阵秩: {controllability_rank}")
#可控性矩阵秩: 2
#系统完全能控 ✓
可观性分析
C = [[1,0]] # 输出矩阵
observability_matrix = [[1,0],[0,1]]
print(f"可观性矩阵秩: {rank(observability_matrix)}")
#可观性矩阵秩: 2
#系统完全能观 ✓
示例3. 机器学习特征分析
特征矩阵的秩分析
feature_matrix =
[[1, 2, 3, 4],
[2, 4, 6, 8],
[3, 6, 9, 12],
[4, 8, 12, 16]]
rank = rank(feature_matrix)
print(f"特征矩阵秩: {rank}")
#特征矩阵秩: 1
#警告：存在线性相关特征，需要特征选择
实际数据集示例（模拟）
real_features =
[[1.0, 2.1, 3.2],
[2.1, 4.0, 6.1],
[3.0, 6.1, 9.0],
[4.2, 8.1, 12.3]]
print(f"真实特征矩阵秩: {rank(real_features)}")
#真实特征矩阵秩: 3
#特征线性独立，适合建模
主成分分析前的秩检验
high_dim_data = [[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7]]
print(f"高维数据秩: {rank(high_dim_data)}")
#高维数据秩: 2
#数据实际存在于低维子空间
示例4. 网络和图论应用
邻接矩阵的秩与连通分量
connected_graph =
[[0,1,0,1],
[1,0,1,0],
[0,1,0,1],
[1,0,1,0]]
print(f"连通图邻接矩阵秩: {rank(connected_graph)}")
#连通图邻接矩阵秩: 2
拉普拉斯矩阵分析
laplacian_matrix =
[[2,-1,0,-1],
[-1,2,-1,0],
[0,-1,2,-1],
[-1,0,-1,2]]
laplacian_rank = rank(laplacian_matrix)
print(f"拉普拉斯矩阵秩: {laplacian_rank}")
print(f"连通分量数: {4 - laplacian_rank}")
#拉普拉斯矩阵秩: 3
#连通分量数: 1
示例5. 结构工程分析
刚度矩阵的奇异性检验（判断结构是否稳定）
stiffness_matrix_stable =
[[2,-1,0],
[-1,2,-1],
[0,-1,2]]
print(f"稳定结构刚度矩阵秩: {rank(stiffness_matrix_stable)}")
#稳定结构刚度矩阵秩: 3
#结构稳定，无刚体位移
奇异刚度矩阵（机构）
stiffness_matrix_unstable =
[[1,-1],
[-1,1]]
print(f"不稳定结构刚度矩阵秩: {rank(stiffness_matrix_unstable)}")
#不稳定结构刚度矩阵秩: 1
#结构不稳定，存在刚体位移模式
示例6. 经济学投入产出分析
Leontief投入产出矩阵
input_output_matrix =
[[0.2, 0.3, 0.1],
[0.1, 0.4, 0.2],
[0.3, 0.1, 0.3]]
io_rank = rank(input_output_matrix)
print(f"投入产出矩阵秩: {io_rank}")
#投入产出矩阵秩: 3
#经济系统可解，存在均衡解
技术系数矩阵分析
technology_matrix =
[[0.8,-0.3,-0.1],
[-0.1,0.6,-0.2],
[-0.3,-0.1,0.7]]
print(f"技术矩阵秩: {rank(technology_matrix)}")
#技术矩阵秩: 3
示例7. 计算机视觉应用
基础矩阵秩检验（双目视觉）
fundamental_matrix =
[[0,0,0.01],
[0,0,-0.02],
[-0.01,0.02,0]]
print(f"基础矩阵秩: {rank(fundamental_matrix)}")
#基础矩阵秩: 2
#有效的基础矩阵（秩为2）
单应性矩阵
homography_matrix =
[[1.0, 0.1, 10],
[0.05, 0.9, 20],
[0.001, 0.002, 1]]
print(f"单应性矩阵秩: {rank(homography_matrix)}")
#单应性矩阵秩: 3
示例8. 符号矩阵分析
参数化矩阵的秩分析
symbolic_matrix1 = [[a,b],[c,d]]
print(f"2x2符号矩阵秩: {rank(symbolic_matrix1)}")
#2x2符号矩阵秩: 2
特殊情况的符号矩阵
symbolic_matrix2 = [[x,2*x],[y,2*y]]
print(f"线性相关符号矩阵秩: {rank(symbolic_matrix2)}")
#线性相关符号矩阵秩: 1
条件满秩的符号矩阵
symbolic_matrix3 = [[1,x],[x,1]]
print(f"条件满秩矩阵秩: {rank(symbolic_matrix3)}")
#条件满秩矩阵秩: 2
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def rank_of_matrix(input_str):
"""
计算矩阵秩，支持符号和数值矩阵，对标MATLAB的rank函数
参数:
input_str (str): 矩阵输入，支持格式：
- "[[1,2],[3,4]]" # 数值矩阵
- "Matrix([[x,1],[2,y]])" # 符号矩阵
tol (float): 可选参数，指定秩计算的容差
返回:
int: 矩阵的秩
str: 错误信息
示例:
>>> rank_of_matrix("[[1,2],[3,4]]") # 数值矩阵
>>> rank_of_matrix("Matrix([[x,1],[2,y]]") # 符号矩阵
>>> rank_of_matrix("[[1,0.999],[2,2]]") # 接近奇异的数值矩阵
"""
try:
# 转换为SymPy矩阵
expr = sp.sympify(input_str)
mat = sp.Matrix(expr) if isinstance(expr, list) else None
if mat is None:
return "错误: 输入无法转换为矩阵"
# 数值矩阵使用SVD计算（类似MATLAB）
if mat.is_symbolic() is False:
# 转换为numpy数组
np_mat = np.array(mat.tolist(), dtype=float)
# 计算SVD
s = np.linalg.svd(np_mat, compute_uv=False)
# 自动计算容差（MATLAB默认算法）
max_singular = np.max(s) if s.size > 0 else 0
tol = max(np_mat.shape) * np.spacing(max_singular)
# 统计大于容差的奇异值数量
rank = np.sum(s > tol)
# 处理数值精度问题
if rank == 0 and not np.allclose(np_mat, 0):
return 0 if np.allclose(np_mat, 0) else 1
return int(rank)
# 符号矩阵使用SymPy默认方法
else:
return mat.rank()
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 正常用例测试
print("=== 正常用例 ===")
print("1. 满秩矩阵:", rank_of_matrix("[[1,2],[3,4]]"))
# 2
print("2. 秩亏矩阵:", rank_of_matrix("[[1,2],[2,4]]"))
# 1
print("3. 符号矩阵:", rank_of_matrix("[[x,1],[2,y]]"))
# 2
# 数值边界测试
print("\n=== 数值边界 ===")
print("4. 零矩阵:", rank_of_matrix("[[0,0],[0,0]]"))
# 0
print("5. 接近奇异的矩阵:",
rank_of_matrix("[[1, 0.9999999999], [2, 2]]"))
# 2
有理分式近似值
R = rat(X) 返回 X 在默认容差 1e-6*norm(X(:),1) 内的有理分式近似值. 近似值是一个包含截断的连分式展开项的字符向量.
R = rat(X,tol) 在容差 tol 内求 X 的近似值
[N,D] = rat(___) 使用以上任何语法返回两个数组 N 和 D, 以使 N./D 求 X 的近似值
X — 输入数组,数值数组
tol — 容差,标量
R — 连分式,字符数组
N — 分子,数值数组
D — 分母,数值数组
示例1. 数学常数近似
print("π 近似:", rat(3.1415926535))
print("e 近似:", rat(2.718281828))
print("黄金比例:", rat(1.618033988))
#π 近似: 355/113
#e 近似: 1457/536
#黄金比例: 1597/987
示例2. 工程应用 - 电阻分压比
print("1:3分压:", rat(0.33333))
print("2:5分压:", rat(0.4))
print("3:8分压:", rat(0.375))
#1:3分压: 1/3
#2:5分压: 2/5
#3:8分压: 3/8
示例3. 金融计算 - 利率和比率
print("8.75%利率:", rat(0.0875))
print("16.67%比率:", rat(0.166667))
print("33.33%折扣:", rat(0.3333))
#8.75%利率: 7/80
#16.67%比率: 1/6
#33.33%折扣: 1/3
示例4. 机械设计 - 齿轮传动比
print("2.4传动比:", rat(2.4))
print("3.75减速比:", rat(3.75))
print("1.6增速比:", rat(1.6))
#2.4传动比: 12/5
#3.75减速比: 15/4
#1.6增速比: 8/5
示例5. 化学计量 - 摩尔比
print("水H2O氢氧比:", rat(0.1119))
print("CO2碳氧比:", rat(0.2727))
#水H2O氢氧比: 63/563
#CO2碳氧比: 3/11
示例6. 带容差的精确控制
print("高精度π:", rat(3.1415926535, 0.0001))
print("宽松近似:", rat(0.333, 0.1))
print("严格近似:", rat(0.333333, 0.000001))
#高精度π: 355/113
#宽松近似: 1/3
#严格近似: 333333/1000000
示例7. 矩阵应用 - 变换矩阵有理化
transform_matrix = [[0.866, -0.5, 0.5], [0.5, 0.866, -0.5], [0, 0, 1]]
N, D = rat(transform_matrix)
print("旋转变换矩阵分子:\n", N)
print("旋转变换矩阵分母:\n", D)
#旋转变换矩阵分子:
#[[433, -1, 1],
[1, 433, -1],
[0, 0, 1]]
#旋转变换矩阵分母:
#[[500, 2, 2],
[2, 500, 2],
[1, 1, 1]]
示例8. 电路分析 - 阻抗矩阵
impedance_matrix = [[1.414, 3.142], [2.718, 4.669]]
N, D = rat(impedance_matrix)
print("阻抗矩阵分子:\n", N)
print("阻抗矩阵分母:\n", D)
#阻抗矩阵分子:
#[[707, 1571],
[1359, 4669]]
#阻抗矩阵分母:
#[[500, 500],
[500, 1000]]
示例9. 控制系统 - 状态空间矩阵
state_matrix = [[0.7071, 0.7071], [-0.7071, 0.7071]] # 45度旋转
N, D = rat(state_matrix)
print("状态矩阵分子:\n", N)
print("状态矩阵分母:\n", D)
#状态矩阵分子:
#[[239, 239],
[-239, 239]]
#状态矩阵分母:
#[[338, 338],
[338, 338]]
示例10. 概率统计 - 转移概率矩阵
prob_matrix = [[0.333, 0.667], [0.25, 0.75]]
N, D = rat(prob_matrix)
print("概率矩阵分子:\n", N)
print("概率矩阵分母:\n", D)
#概率矩阵分子:
#[[333, 667],
[1, 3]]
#概率矩阵分母:
#[[1000, 1000],
[4, 4]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def rational_approximation(input_str):
"""
数值有理分式近似，对标MATLAB的rat函数
参数:
input_str (str): 输入表达式，支持格式：
- "x" : 单个数值近似
- "(x, tol)" : 带容差的单个数值近似
- "[[1.25,3.14],[2.71,4.44]]" : 矩阵近似
- "(matrix, tol)" : 带容差的矩阵近似
返回:
tuple/float: 矩阵返回(N,D)元组，标量返回分数形式
str: 错误信息
示例:
>>> rational_approximation("0.333") # 返回 1/3
>>> rational_approximation("(0.333, 0.01)") # 返回 1/3
>>> rational_approximation("[[0.333,1.414],[3.141,2.718]]")
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def rat_scalar(x, tol=0.001):
"""处理标量有理近似"""
return sp.nsimplify(x, tolerance=tol, rational=True)
# 处理元组输入 (参数, 容差)
if isinstance(expr, tuple):
if len(expr) != 2:
raise ValueError("元组参数必须包含两个元素")
# 标量+容差模式
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = rat_scalar(*params)
else:
error = True
# 处理标量输入
elif expr.is_number:
result = rat_scalar(float(expr))
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
print("=== 标量测试 ===")
print("1. 1/3近似:", rational_approximation("0.333333"))
# 1/3
print("2. 带容差:", rational_approximation("(0.333, 0.1)"))
# 1/3
有理输出
S = rats(X) 返回一个包含 X 的元素的有理近似值的字符向量. 星号指示无法在分配的空间中输出, 但与 X 中的其他元素相比并非可忽略不计的元素.
S = rats(X,tol) 在容差 tol 内求 X 的近似值
X — 输入数组,数值数组
tol — 容差,标量
示例1. 基础数值转换
print(f"0.333333 → {rats(0.333333)}")
print(f"0.666667 → {rats(0.666667)}")
print(f"0.142857 → {rats(0.142857)}")
print(f"0.875 → {rats(0.875)}")
#0.333333 → 1/3
#0.666667 → 2/3
#0.142857 → 1/7
#0.875 → 7/8
示例2. 数学常数近似
print(f"π ≈ 3.14159 → {rats(3.14159)}")
print(f"e ≈ 2.71828 → {rats(2.71828)}")
print(f"√2 ≈ 1.41421 → {rats(1.41421)}")
print(f"φ ≈ 1.61803 → {rats(1.61803)}")
#π ≈ 3.14159 → 355/113
#e ≈ 2.71828 → 1264/465
#√2 ≈ 1.41421 → 1055/746
#φ ≈ 1.61803 → 1364/843
示例3. 不同容差下的转换
print(f"π (tol=0.1) → {rats(3.14159, 0.1)}")
print(f"π (tol=0.01) → {rats(3.14159, 0.01)}")
print(f"π (tol=0.001) → {rats(3.14159, 0.001)}")
print(f"π (tol=1e-5) → {rats(3.14159, 0.00001)}")
#π (tol=0.1) → 22/7
#π (tol=0.01) → 311/99
#π (tol=0.001) → 355/113
#π (tol=1e-5) → 314159/100000
示例4. 矩阵转换 - 电路分析中的阻抗矩阵
impedance_matrix = [[1.5, 0.333, 2.25], [0.6667, 3.142, 1.414], [4.444, 2.718, 0.875]]
result = rats(impedance_matrix)
阻抗矩阵:
for row in result.tolist():
print(f" {row}")
#[3/2, 333/1000, 9/4]
#[2/3, 1571/500, 707/500]
#[1111/250, 1359/500, 7/8]
示例5. 控制系统中的传递函数系数
control_matrix = [[0.1667, 0.8333], [0.4286, 0.5714]]
result = rats(control_matrix)
状态空间矩阵:
for row in result.tolist():
print(f" {row}")
#[1/6, 5/6]
#[3/7, 4/7]
示例6. 概率统计中的概率值
probabilities = [[0.25, 0.125, 0.0625], [0.333, 0.1667, 0.1111]]
result = rats(probabilities)
概率矩阵:
for row in result.tolist():
print(f" {row}")
#[1/4, 1/8, 1/16]
#[333/1000, 1/6, 1/9]
示例7. 机械工程中的齿轮比
gear_ratios = [[2.4, 3.6], [1.8, 4.2]]
result = rats(gear_ratios)
齿轮比矩阵:
for row in result.tolist():
print(f" {row}")
#[12/5, 18/5]
#[9/5, 21/5]
示例8. 化学计量系数
stoichiometry = [[0.5, 1.333, 2.667], [0.75, 1.5, 3.0]]
result = rats(stoichiometry)
化学计量矩阵:
for row in result.tolist():
print(f" {row}")
#[1/2, 1333/1000, 2667/1000]
#[3/4, 3/2, 3]
示例9. 金融计算中的比率
financial_ratios = [[1.333, 2.5], [0.75, 1.6667]]
result = rats(financial_ratios)
财务比率矩阵:
for row in result.tolist():
print(f" {row}")
#[1333/1000, 5/2]
#[3/4, 5/3]
示例10. 带容差的矩阵转换
filter_coeffs = ([[0.1234, 0.5678, 0.9012], [0.3456, 0.7890, 0.2345]], 0.01)
result = rats(filter_coeffs)
滤波器系数矩阵:
for row in result.tolist():
print(f" {row}")
#[10/81, 46/81, 73/81]
#[28/81, 71/90, 19/81]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def rational_approximation_strlen(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
def rat_to_str(x, tol=0.001):
"""将数值转换为分数字符串"""
rational = sp.nsimplify(x, tolerance=tol, rational=True)
return rational
# 处理元组输入（带容差）
if isinstance(expr, tuple):
if len(expr) != 2:
raise ValueError("元组参数必须为 (输入值, 容差)")
input_val, tol_val = expr[0], float(expr[1])
# 处理标量输入
if input_val.is_number:
result = rat_to_str(float(input_val), tol_val)
else:
error = True
# 处理标量输入
elif expr.is_number:
result = rat_to_str(float(expr))
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
print("=== 标量测试 ===")
print("1. 基础测试:", rational_approximation_strlen("0.333333"))
# 1/3
print("2. 整数测试:", rational_approximation_strlen("2"))
# 2
print("3. 带容差测试:", rational_approximation_strlen("(3.1415926, 0.001)"))
# 355/113
有理模型拟合
有理模型的分子和分母均为多项式,其中分母的首项系数设置为1.模型名称为ratij,其中i是分子的次数,j是分母的次数.分子和分母的次数最高都为5.
rat02: Y = (p1)/(x^2+q1*x+q2)
rat21: Y = (p1*x^2+p2*x+p3)/(x+q1)
rat55: Y = (p1*x^5+...+p6)/(x^5+...+q5)
返回有理模型拟合方程式
示例1. 电路分析 - RC电路频率响应
模拟RC低通滤波器的频率响应: H(f) = 1/(1 + j2πfRC)
frequencies = np.logspace(-1, 2, 20) # 0.1Hz to 100Hz
RC = 0.01 # 时间常数
magnitude = 1 / np.sqrt(1 + (2 * np.pi * frequencies * RC)**2)
添加一些噪声
magnitude_noisy = magnitude * (1 + 0.05 * np.random.normal(size=len(magnitude)))
x_data = frequencies.tolist()
y_data = magnitude_noisy.tolist()
expr = ([x_data], [y_data], 2)
result = ratfit(expr)
print(f"RC电路频率响应拟合: {result}")
print(f"理论模型: 1/sqrt(1 + (2πf×0.01)²)")
#RC电路频率响应拟合: 1213.060195/(x**2 + 12.391964*x + 1190.676477)
#理论模型: 1/sqrt(1 + (2πf×0.01)²)
示例2. 机械系统 - 弹簧质量阻尼系统
模拟质量-弹簧-阻尼系统的传递函数
omega_n = 10 # 自然频率 (rad/s)
zeta = 0.1 # 阻尼比
frequencies_mech = np.linspace(0.1, 20, 25)
二阶系统的幅频特性: |H(ω)| = 1/√((1-(ω/ω_n)²)² + (2ζω/ω_n)²)
magnitude_mech = 1 / np.sqrt((1 - (frequencies_mech/omega_n)**2)**2 + (2*zeta*frequencies_mech/omega_n)**2)
magnitude_mech_noisy = magnitude_mech * (1 + 0.03 * np.random.normal(size=len(magnitude_mech)))
expr = ([frequencies_mech], [magnitude_mech_noisy], 2)
result = ratfit(expr)
print(f"机械系统拟合: {result}")
print(f"理论模型: 1/√((1-(ω/10)²)² + (0.2ω/10)²)")
#机械系统拟合: 29.272331/(x**2 - 19.217912*x + 99.233675)
#理论模型: 1/√((1-(ω/10)²)² + (0.2ω/10)²)
示例3. 化学动力学 - 酶催化反应Michaelis-Menten 动力学
Michaelis-Menten 方程: v = Vmax * [S] / (Km + [S])
substrate_concentration = np.linspace(0.1, 10, 20)
Vmax = 2.0
Km = 1.5
reaction_rate = Vmax * substrate_concentration / (Km + substrate_concentration)
reaction_rate_noisy = reaction_rate * (1 + 0.02 * np.random.normal(size=len(reaction_rate)))
使用 21 模型 (二次/一次) 来拟合 Michaelis-Menten 方程
expr = ([substrate_concentration], [reaction_rate_noisy], 21)
result = ratfit(expr)
print(f"酶动力学拟合: {result}")
print(f"理论模型: 2.0*[S]/(1.5 + [S])")
#动力学拟合: (0.007351*x**2 + 1.942883*x - 0.00011)/(x + 1.445078)
#理论模型: 2.0*[S]/(1.5 + [S])
示例4. 光学系统 - 透镜成像公式
薄透镜公式: 1/f = 1/u + 1/v
重新整理为: v = f*u/(u-f)
object_distance = np.linspace(1.5, 5, 15) # 物距
focal_length = 1.0
image_distance = focal_length * object_distance / (object_distance - focal_length)
image_distance_noisy = image_distance + 0.1 * np.random.normal(size=len(image_distance))
expr = ([object_distance], [image_distance_noisy], 21)
result = ratfit(expr)
print(f"透镜成像拟合: {result}")
print(f"理论模型: 1.0*u/(u-1.0)")
#透镜成像拟合: (-0.018153*x**2 + 1.2425*x - 0.788964)/(x - 1.137934)
#理论模型: 1.0*u/(u-1.0)
示例5. 经济学 - 需求弹性曲线
模拟价格弹性需求函数
prices = np.linspace(1, 20, 18)
需求函数: Q = a/(p + b) + c
a, b, c = 100, 2, 5
demand = a/(prices + b) + c
demand_noisy = demand + 2 * np.random.normal(size=len(demand))
expr = ([prices], [demand_noisy], 21)
result = ratfit(expr)
print(f"需求曲线拟合: {result}")
print(f"理论模型: 100/(p+2) + 5")
#需求曲线拟合: (-0.34233*x**2 + 15.277235*x + 22.630358)/(x - 0.043586)
#理论模型: 100/(p+2) + 5
示例6. 材料科学 - 粘弹性材料频率响应
模拟标准线性固体模型的频率响应
freq_material = np.logspace(-2, 3, 25)
复数模量的实部: E'(ω) = E_inf + (E_0 - E_inf)*ω²τ²/(1+ω²τ²)
E_inf, E_0, tau = 10, 1, 0.1
storage_modulus = E_inf + (E_0 - E_inf) * (freq_material**2 * tau**2) / (1 + freq_material**2 * tau**2)
storage_modulus_noisy = storage_modulus * (1 + 0.04 * np.random.normal(size=len(storage_modulus)))
expr = ([freq_material], [storage_modulus_noisy], 21)
result = ratfit(expr)
print(f"粘弹性材料拟合: {result}")
print(f"理论模型: 10 + (1-10)*ω²×0.01/(1+ω²×0.01)")
#粘弹性材料拟合: (0.001719*x**2 - 0.405037*x + 128.651422)/(x + 12.353449)
#理论模型: 10 + (1-10)*ω²×0.01/(1+ω²×0.01)
示例7. 高阶系统 - 五阶巴特沃斯滤波器
模拟五阶低通滤波器的频率响应
freq_filter = np.logspace(-1, 2, 30)
五阶巴特沃斯滤波器的近似幅频响应
w = 2 * np.pi * freq_filter
w_c = 10 # 截止频率
magnitude_filter = 1 / np.sqrt(1 + (w/w_c)**10)
magnitude_filter_noisy = magnitude_filter * (1 + 0.02 * np.random.normal(size=len(magnitude_filter)))
使用55模型进行高阶拟合
expr = ([freq_filter], [magnitude_filter_noisy], 55)
result = ratfit(expr)
print(f"五阶滤波器拟合: {result}")
print("理论模型: 五阶巴特沃斯滤波器响应")
#五阶滤波器拟合: (-0.002231*x**5 + 0.105323*x**4 - 1.034238*x**3 + 3.548568*x**2 - 5.448504*x + 0.773382)/(x*(1.952358*x**4 - 4.034387*x**3 + 4.833161*x**2 - 5.614258*x + 0.777087))
#理论模型: 五阶巴特沃斯滤波器响应
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.optimize import curve_fit
def rat_fit_nonlinear(input_str):
"""
有理函数非线性拟合，支持多种分子分母阶数组合
参数:
input_str (str): 输入表达式字符串，格式为：
"(x_data, y_data, model_code)"
其中：
- x_data: 自变量数据矩阵/列表
- y_data: 因变量数据矩阵/列表
- model_code: 模型代码，支持：
2 : p1/(x² + q1x + q2)
21 : (p1x² + p2x + p3)/(x + q1)
55 : 五阶分子/分母多项式
返回:
sympy.Expr: 拟合得到的有理函数表达式
str: 错误信息
示例:
>>> rat_fit_nonlinear("([[1,2,3]], [[0.5,0.3,0.2]], 2)")
>>> rat_fit_nonlinear("([[0.1,0.5,1]], [[1.2,0.8,0.5]], 21)")
"""
try:
expr = sp.sympify(input_str)
error = False
expression = None
maxfev = 100000 # 最大函数评估次数
# 定义有理函数模型 --------------------------------------------------------
def rat_model2(x, p1, q1, q2):
"""(2阶分母模型) p1 / (x² + q1*x + q2)"""
return p1 / (x ** 2 + q1 * x + q2)
def rat_model21(x, p1, p2, p3, q1):
"""(2阶分子/1阶分母模型) (p1x² + p2x + p3)/(x + q1)"""
return (p1 * x ** 2 + p2 * x + p3) / (x + q1)
def rat_model55(x, p1, p2, p3, p4, p5, p6, q1, q2, q3, q4, q5):
"""(5阶分子/5阶分母模型) (p1x⁵+...+p6)/(x⁵+...+q5)"""
numerator = p1 * x ** 5 + p2 * x ** 4 + p3 * x ** 3 + p4 * x ** 2 + p5 * x + p6
denominator = x ** 5 + q1 * x ** 4 + q2 * x ** 3 + q3 * x ** 2 + q4 * x + q5
return numerator / denominator
# 参数解析 -------------------------------------------------------------
if isinstance(expr, tuple):
# 解析输入数据
if len(expr) >= 3:
x, y, n = expr[0], expr[1], int(expr[2])
else: # 默认使用模型2
x, y = expr[0], expr[1]
n = 2
# 转换输入矩阵为numpy数组
x_data = np.array(x).ravel() if isinstance(x, list) else None
y_data = np.array(y).ravel() if isinstance(y, list) else None
if x_data is None or y_data is None:
raise ValueError("输入数据格式错误")
# 选择拟合模型 -----------------------------------------------------
if n == 2:
# 初始参数猜测 [p1, q1, q2]
initial_guess = [1.0, 0.5, 0.5]
params, _ = curve_fit(rat_model2, x_data, y_data,
p0=initial_guess, maxfev=maxfev)
p1, q1, q2 = params
expression = f"{p1:.6f}/(x**2 + {q1:.6f}*x + {q2:.6f})"
elif n == 21:
# 初始参数猜测 [p1, p2, p3, q1]
initial_guess = [1.0, -0.5, 0.3, 0.5]
params, _ = curve_fit(rat_model21, x_data, y_data,
p0=initial_guess, maxfev=maxfev)
p1, p2, p3, q1 = params
expression = f"({p1:.6f}*x**2 + {p2:.6f}*x + {p3:.6f})/(x + {q1:.6f})"
elif n == 55:
# 初始参数猜测 [p1-6, q1-5]
initial_guess = [1.0] + [0.5] * 5 + [0.3] * 5
params, _ = curve_fit(rat_model55, x_data, y_data,
p0=initial_guess, maxfev=maxfev)
p_params = params[:6]
q_params = params[6:]
num_str = " + ".join([f"{p:.6f}*x**{5 - i}" for i, p in enumerate(p_params)])
den_str = " + ".join([f"{q:.6f}*x**{5 - i}" for i, q in enumerate(q_params)])
expression = f"({num_str})/(x**5 + {den_str})"
else:
raise ValueError(f"不支持的模型代码: {n}")
return sp.simplify(expression) if not error else f"输入错误: {input_str}"
else:
return f"输入格式错误，应为元组"
except Exception as e:
return f"拟合失败: {str(e)}"
# 示例测试 =======================================================================
if __name__ == "__main__":
# 测试用例1: 二阶分母模型
print("=== 测试1: 二阶分母模型 ===")
x_data = [[1, 2, 3, 4, 5]]
y_data = [[1 / (i ** 2 + 2 * i + 3) for i in x_data[0]]]
expr = f"({x_data}, {y_data}, 2)"
print("拟合结果:", rat_fit_nonlinear(expr))
# 1.0/(x**2 + 2.0*x + 3.0)
# 测试用例2: 二次/一次模型
print("\n=== 测试2: 二次/一次模型 ===")
x_data = [[0.5, 1.0, 1.5, 2.0]]
y_data = [[(2 * i ** 2 - 3 * i + 1) / (i + 2) for i in x_data[0]]]
expr = f"({x_data}, {y_data}, 21)"
print("拟合结果:", rat_fit_nonlinear(expr))
# (2.0*x**2 - 3.0*x + 1.0)/(x + 2.0)
条件数倒数
C = rcond(A) 返回 A 的 1-范数条件数倒数估计值. 如果A的条件设置良好, rcond(A)接近1.0. 如果A是病态的, 则rcond(A)接近 0
A — 输入矩阵, 数值方阵
C — 条件数倒数, 标量
示例1. 理想情况 - 单位矩阵和正交矩阵
print(f"单位矩阵 2x2: {rcond([[1, 0], [0, 1]])}")
#单位矩阵 2x2: 1.0
print(f"单位矩阵 3x3: {rcond([[1, 0, 0], [0, 1, 0], [0, 0, 1]])}")
#单位矩阵 3x3: 1.0
print(f"正交矩阵: {rcond([[0.6, -0.8], [0.8, 0.6]])}")
#正交矩阵: 0.5102040816326532
print(f"对角矩阵: {rcond([[2, 0, 0], [0, 3, 0], [0, 0, 4]])}")
#对角矩阵: 0.5
示例2. 病态矩阵 - Hilbert矩阵 (经典的病态矩阵)
hilbert_3x3 = [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]]
print(f"Hilbert 3x3: {rcond(hilbert_3x3)}")
#Hilbert 3x3: 0.0013368983957219203
hilbert_4x4 = [[1, 1/2, 1/3, 1/4], [1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7]]
print(f"Hilbert 4x4: {rcond(hilbert_4x4)}")
#Hilbert 4x4: 3.5242290748899124e-05
hilbert_5x5 = [[1, 1/2, 1/3, 1/4, 1/5], [1/2, 1/3, 1/4, 1/5, 1/6], [1/3, 1/4, 1/5, 1/6, 1/7], [1/4, 1/5, 1/6, 1/7, 1/8], [1/5, 1/6, 1/7, 1/8, 1/9]]
print(f"Hilbert 5x5: {rcond(hilbert_5x5)}")
#Hilbert 5x5: 1.0597081987501896e-06
示例3. 数值分析中的典型矩阵
Vandermonde矩阵 (多项式插值中常见)
vandermonde = [[1, 1, 1], [1, 2, 4], [1, 3, 9]]
print(f"Vandermonde矩阵: {rcond(vandermonde)}")
#Vandermonde矩阵: 0.008928571428571428
对称正定矩阵
spd_matrix = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]] # 离散拉普拉斯算子
print(f"对称正定矩阵: {rcond(spd_matrix)}")
#对称正定矩阵: 0.125
Toeplitz矩阵
toeplitz = [[1, 0.5, 0.25], [0.5, 1, 0.5], [0.25, 0.5, 1]]
print(f"Toeplitz矩阵: {rcond(toeplitz)}")
#Toeplitz矩阵: 0.16666666666666669
示例4. 工程应用中的矩阵
结构力学中的刚度矩阵
stiffness = [[3, -1, -1], [-1, 3, -1], [-1, -1, 3]]
print(f"刚度矩阵: {rcond(stiffness)}")
#刚度矩阵: 0.2
电路分析中的导纳矩阵
admittance = [[1.5, -0.5, -0.5], [-0.5, 1.5, -0.5], [-0.5, -0.5, 1.5]]
print(f"导纳矩阵: {rcond(admittance)}")
#导纳矩阵: 0.2
控制系统中的状态矩阵
state_matrix = [[0.9, 0.1], [-0.2, 0.8]]
print(f"状态矩阵: {rcond(state_matrix)}")
#状态矩阵: 0.6727272727272726
示例5. 接近奇异的矩阵
nearly_singular1 = [[1, 1], [1, 1.000001]]
print(f"接近奇异1: {rcond(nearly_singular1)}")
#接近奇异1: 2.4999974997962085e-07
nearly_singular2 = [[1, 2, 3], [4, 5, 6], [7, 8, 9.000001]]
print(f"接近奇异2: {rcond(nearly_singular2)}")
#接近奇异2: 6.944444921500388e-09
nearly_singular3 = [[0.001, 0], [0, 0.001]]
print(f"接近奇异3: {rcond(nearly_singular3)}")
#接近奇异3: 1.0
示例6. 随机矩阵的条件数分析
np.random.seed(42) # 固定随机种子以便重现
良态随机矩阵
well_conditioned_random = np.random.rand(3, 3) + 10 * np.eye(3)
print(f"良态随机矩阵: {rcond(well_conditioned_random)}")
#良态随机矩阵: 0.7221874757898497
病态随机矩阵 (列几乎线性相关)
ill_conditioned_random = np.array([[1, 1.001, 0.999],
[2, 2.002, 1.998],
[3, 3.003, 2.997]])
print(f"病态随机矩阵: {rcond(ill_conditioned_random)}")
#病态随机矩阵: 0.0
示例7. 线性方程组求解的稳定性分析
A1 = [[10, -1, 2, 0], [-1, 11, -1, 3], [2, -1, 10, -1], [0, 3, -1, 8]]
rcond1 = rcond(A1)
print(f"良态方程组矩阵: rcond = {rcond1}")
print(f" 可安全求解: {rcond1 > 1e-10}")
#良态方程组矩阵: rcond = 0.31875
#可安全求解: True
A2 = [[1, 2, 3], [4, 5, 6], [7, 8, 9]] # 秩亏缺
rcond2 = rcond(A2)
print(f"病态方程组矩阵: rcond = {rcond2}")
print(f" 需要正则化: {rcond2 < 1e-12}")
#病态方程组矩阵: rcond = 1.5419764230904951e-18
#需要正则化: True
示例8. 矩阵求逆的数值稳定性
matrices = [
("良态可逆矩阵", [[4, 1], [1, 4]]),
("边界情况", [[1, 2], [2, 4.0001]]),
("不可逆矩阵", [[1, 2], [2, 4]])
]
for name, matrix_str in matrices:
rcond = rcond(matrix_str)
status = "稳定" if rcond > 1e-8 else "不稳定" if rcond > 1e-12 else "不可逆"
print(f" {name}: rcond = {rcond:.2e} ({status})")
#良态可逆矩阵: rcond = 6.00e-01 (稳定)
#边界情况: rcond = 2.78e-06 (稳定)
#不可逆矩阵: rcond = 0.00e+00 (不可逆)
示例9. 特征值问题的条件数
对称矩阵通常有较好的特征值条件数
symmetric = [[2, 1, 0], [1, 3, 1], [0, 1, 2]]
print(f"对称矩阵: {rcond(symmetric)}")
#对称矩阵: 0.2
非对称矩阵可能对特征值敏感
nonsymmetric = [[0.9, 0.1], [0.1, 0.8]]
print(f"非对称矩阵: {rcond(nonsymmetric)}")
#非对称矩阵: 0.7100000000000001
示例10. 有限元分析中的刚度矩阵
二维三角形单元的刚度矩阵
stiffness_2d = [[2, -1, -1, 0], [-1, 2, 0, -1], [-1, 0, 2, -1], [0, -1, -1, 2]]
rcond_stiffness = rcond(stiffness_2d)
print(f"刚度矩阵条件数倒数: {rcond_stiffness}")
print(f" 网格质量: {'优秀' if rcond_stiffness > 0.1 else '良好' if rcond_stiffness > 0.01 else '需要改进'}")
#刚度矩阵条件数倒数: 0.0
# 网格质量: 需要改进
示例11. 电路网络分析
RLC网络的导纳矩阵
rlc_admittance = [[2+1@i, -1-0.5@i, 0], [-1-0.5@i, 3+2@i, -1-0.5@i], [0, -1-0.5@i, 2+1@i]]
注意：当前函数不支持复数，这里使用实数近似
rlc_real = [[2, -1, 0], [-1, 3, -1], [0, -1, 2]]
rcond_circuit = rcond(rlc_real)
print(f"电路导纳矩阵: {rcond_circuit}")
print(f" 数值稳定性: {'高' if rcond_circuit > 0.1 else '中等' if rcond_circuit > 0.001 else '低'}")
#电路导纳矩阵: 0.2
# 数值稳定性: 高
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def reverse_condition_norm(input_str):
"""
计算矩阵条件数的倒数（对标MATLAB的rcond函数）
注意:
- 此处使用精确计算而非估计，与MATLAB的rcond实现方式不同
- MATLAB的rcond使用1-范数估计，而这里计算精确的1-范数条件数倒数
参数:
input_str: 矩阵的字符串表示（如"[[1, 0], [0, 1]]"）或Matrix对象
返回:
float: 条件数的倒数
str: 错误时返回错误信息
"""
try:
# 将输入转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 尝试转换为矩阵
if isinstance(expr, list):
# 转换为NumPy数组进行数值计算
A = np.array(expr, dtype=float)
# 计算1-范数条件数
cond_A = np.linalg.cond(A, p=1)
# 返回条件数的倒数
result = 1 / cond_A
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例代码
if __name__ == "__main__":
# 示例1：单位矩阵（条件数=1，倒数=1）
print("示例1:", reverse_condition_norm("[[1, 0], [0, 1]]"))
# 1.0
# 示例2：接近奇异的矩阵
ill_conditioned_matrix = "[[1, 1], [1, 1.0001]]"
print("示例2:", reverse_condition_norm(ill_conditioned_matrix))
# 2.4997500187484752e-05
# 示例3：直接使用Python列表输入
print("示例3:", reverse_condition_norm([[1, 2], [3, 4]]))
# 0.04761904761904764
# 示例4：使用SymPy矩阵语法
print("示例5:", reverse_condition_norm("[[2, 0], [0, 2]]"))
# 1.0
复数的实部
X = real(Z) 返回数组 Z 中每个元素的实部
Z - 输入数组,标量,向量,矩阵
示例1. 基础复数操作
print(f"3+4i → {real(3 + 4@i)}")
print(f"2-3i → {real(2 - 3@i)}")
print(f"-1+2i → {real(-1 + 2@i)}")
print(f"纯实数 5 → {real(5)}")
print(f"纯虚数 3i → {real(3@i)}")
#3+4i → 3
#2-3i → 2
#-1+2i → -1
#纯实数 5 → 5
#纯虚数 3i → 0
示例2. 复数矩阵操作
complex_matrix1 = [[1+2@i, 3-4@i], [5+6@i, 7-8@i]]
result1 = real(complex_matrix1)
print(f"2x2复数矩阵: {result1}")
#2x2复数矩阵: [[1, 3], [5, 7]]
complex_matrix2 = [[1, 2+3@i, 4], [5@i, 6, 7-8@i]]
result2 = real(complex_matrix2)
print(f"混合矩阵: {result2}")
#混合矩阵: [[1, 2, 4], [0, 6, 7]]
real_matrix = [[1, 2, 3], [4, 5, 6]]
result3 = real(real_matrix)
print(f"实数矩阵: {result3}")
#实数矩阵: [[1, 2, 3], [4, 5, 6]]
示例3. 工程应用中的复数
电路分析 - 阻抗
impedance = 10 + 5@i # 10Ω电阻 + 5Ω感抗
print(f"电路阻抗 {impedance} → {real(impedance)}")
#电路阻抗 10 + 5*I → 10
信号处理 - 复数信号
signal = 3*exp(@pi/4*@i) # 幅度3, 相位π/4
print(f"复数信号 {signal} → {real(signal)}")
#复数信号 3*exp(I*pi/4) → 3*sqrt(2)/2
控制系统 - 传递函数在特定频率的值
tf_value = 1/(1 + 0.1@i) # 在频率点处的传递函数值
print(f"传递函数值 {tf_value} → {real(tf_value)}")
#传递函数值 1/(1 + 0.1*I) → 0.990099009900990
示例4. 交流电路分析
RLC串联电路阻抗: Z = R + j(ωL - 1/ωC)
R, L, C, omega = 10, 0.1, 0.001, 100 # 参数
Z_RLC = R + omega * L - 1 / (omega * C)*@i
print(f"RLC电路阻抗: {Z_RLC}")
print(f"电阻分量: {real(Z_RLC)} Ω")
#RLC电路阻抗: 10 + 0.0*I
#电阻分量: 10 Ω
示例5. 交流功率计算
复功率 S = P + jQ, 其中P为有功功率(实部), Q为无功功率(虚部)
voltage = 220*exp(@pi/6*@i) # 电压相量
current = 10*exp(-@pi/12*@i) # 电流相量
复功率 S = V × conj(I)
S = (voltage) * conjugate(current)
S_real = real(S)
print(f"复功率 S = {S}")
print(f"有功功率 P = {S_real} W")
#复功率 S = (220*exp(I*pi/6)) * conjugate(10*exp(-I*pi/12))
#有功功率 P = 1100*sqrt(2) W
示例6. 滤波器频率响应
一阶RC低通滤波器: H(ω) = 1/(1 + jωRC)
RC = 0.01
frequencies = [10, 100, 1000] # Hz
for f in frequencies:
omega = 2 * sp.pi * f
H = 1/(1 + omega * RC*@i)
H_real = real(H)
print(f"f={f}Hz: H(ω)实部 = {float(H_real):.4f}")
#f=10Hz: H(ω)实部 = 0.7170
#f=100Hz: H(ω)实部 = 0.0247
#f=1000Hz: H(ω)实部 = 0.0003
示例7. 傅里叶变换系数
复数傅里叶系数表示信号的幅度和相位
fourier_coeffs = [[2+3@i, 1-2@i, 0.5+0.5@i], [-1+1@i, 3, 2-1@i]]
result = real(fourier_coeffs)
傅里叶系数矩阵实部:
print(result)
#[[2, 1, 0.5],
[-1, 3, 2]]
示例8. 复数滤波器系数
FIR滤波器的复数系数
fir_coeffs = [0.1+0.2@i, 0.3-0.1@i, 0.2+0.05@i, 0.1-0.15@i]
real_coeffs = real(fir_coeffs)
print(f"复数滤波器系数实部: {real_coeffs}")
#复数滤波器系数实部: [0.1, 0.3, 0.2, 0.1]
示例9. 频率响应分析
传递函数在特定频率的取值
G = 1 / (s ** 2 + 2 * s + 2) # 二阶系统
在s = jω处的取值
omega = 1
G_jw = G.subs(s, sp.I * omega)
print(f"传递函数 G(s) = {G}")
print(f"在ω={omega}处: G(jω) = {G_jw}")
print(f"实部: {real(G_jw)}")
#传递函数 G(s) = 1/(s**2 + 2*s + 2)
#在ω=1处: G(jω) = (1 - 2*I)/5
#实部: 1/5
示例10. 奈奎斯特图数据
开环传递函数在不同频率的取值
frequencies = [0.1, 1, 10]
for w in frequencies:
L_jw = 1 / ((sp.I * w) ** 2 + 2 * (sp.I * w) + 1)
L_real = real(L_jw)
print(f"ω={w}: 实部 = {float(L_real):.4f}")
#ω=0.1: 实部 = 0.9705
#ω=1: 实部 = 0.0000
#ω=10: 实部 = -0.0097
示例11. 实验数据中的复数处理
假设从测量设备获得的复数数据
experimental_data = [
[1.2 + 0.5@i, 2.1 - 0.3@i, 0.8 + 1.2@i],
[3.2 - 0.8@i, 1.5 + 0.6@i, 2.3 - 1.1@i],
[0.9 + 0.2@i, 1.8 - 0.4@i, 2.7 + 0.9@i]
]
real_data = real(experimental_data)
原始复数数据实部:
print(real_data)
#[[1.2, 2.1, 0.8],
[3.2, 1.5, 2.3],
[0.9, 1.8, 2.7]]
示例12. 从复数数据中提取物理量
例如：复数电压测量值 → 实际电压值
complex_measurements = [10.2+3.5@i, 8.7-2.1@i, 12.3+5.6@i, 9.1-1.8@i]
physical_values = real(complex_measurements)
print(f"测量值实部(物理量): {physical_values}")
#测量值实部(物理量): [10.2, 8.7, 12.3, 9.1]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def re_part_of_complex(input_str):
"""
提取复数或复数矩阵的实部（对标MATLAB的real函数）
参数:
input_str: 输入可以是以下类型：
- 复数/实数的字符串表示（如"3+4*I"）
- 矩阵的字符串表示（如"[[1, 2+I], [3, 4]]"）
- 直接传入Python列表或SymPy矩阵
返回:
sp.Expr/sp.Matrix: 实部计算结果
str: 错误时返回错误信息
注意:
- 支持符号运算，但需要明确声明复数属性（如使用Symbol('x', real=False)）
- 输入矩阵需符合SymPy矩阵语法
"""
try:
# 将输入转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
def eval_re_part(x):
"""递归处理元素的实部提取"""
return sp.re(x) if x.is_complex else x
if expr.is_number or expr.is_symbol:
# 处理标量或符号
result = eval_re_part(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例代码
if __name__ == "__main__":
# 示例1：复数标量
print("示例1:", re_part_of_complex("3 + 4*I"))
# 3
# 示例2：实数标量
print("示例2:", re_part_of_complex("5"))
# 5
非负实数数组的自然对数
Y = reallog(X)返回数组X中每个元素的自然对数.数组X只能包含非负实数.Y的大小与X的大小相同.
示例1：金融计算 - 复利计算
计算不同利率下的终值对数（用于对数收益率计算）
principal = 1000
future_values = [[1100, 1210, 1331], [1200, 1440, 1728]] # 不同投资方案的终值
print("对数收益率矩阵:", reallog(future_values))
#对数收益率矩阵: [[7.00306545878646, 7.09837563859079, 7.19368581839511],
[7.09007683577609, 7.27239839257005, 7.454719949364]]
示例2：物理学 - 声压级计算
声压级计算：Lp = 20*log10(P/P0)，这里计算对数部分
reference_pressure = 2e-5 # 参考声压20μPa
measured_pressures = [[0.002, 0.02, 0.2], [0.0002, 0.002, 0.02]] # 测量的声压(Pa)
print("声压对数:", reallog(measured_pressures))
#声压对数: [[-6.21460809842219, -3.91202300542815, -1.6094379124341],
[-8.51719319141624, -6.21460809842219, -3.91202300542815]]
示例3：化学 - pH值计算
pH = -log10([H+])，这里计算氢离子浓度的对数
hydrogen_ion_concentrations = [1e-7, 1e-6, 1e-5, 1e-4] # 氢离子浓度(mol/L)
print("浓度对数:", reallog(hydrogen_ion_concentrations))
#浓度对数: [-16.1180956509583, -13.8155105579643, -11.5129254649702, -9.21034037197618]
示例4：生物学 - 细胞生长模型
细胞数量增长的指数模型：N(t) = N0 * e^(kt)，取对数得线性关系
cell_counts = [[100, 271, 738], [50, 136, 369]] # 不同时间点的细胞数量
print("细胞数量对数:", reallog(cell_counts))
#细胞数量对数: [[4.60517018598809, 5.6021188208797, 6.60394382460047],
[3.91202300542815, 4.91265488573605, 5.91079664404053]]
示例5：经济学 - 价格指数计算
计算价格指数的对数变化率
price_indices = [[100, 110, 121], [100, 95, 90.25]] # 不同商品的价格指数序列
print("价格指数对数:", reallog(price_indices))
#价格指数对数: [[4.60517018598809, 4.70048036579242, 4.79579054559674],
[4.60517018598809, 4.55387689160054, 4.50258359721299]]
示例6：工程学 - 信号处理
计算信号功率的对数（分贝计算的基础）
signal_powers = [[1, 10, 100], [0.1, 1, 10]] # 信号功率
print("功率对数:", reallog(signal_powers))
#功率对数: [[0, 2.30258509299405, 4.60517018598809],
[-2.30258509299405, 0, 2.30258509299405]]
示例7：统计学 - 数据变换
对偏态分布数据进行对数变换，使其更接近正态分布
skewed_data = [[1, 4, 9, 16, 25], [2, 8, 18, 32, 50]] # 右偏分布数据
print("对数变换后:", reallog(skewed_data))
#对数变换后: [[0, 1.38629436111989, 2.19722457733622, 2.77258872223978, 3.2188758248682],
[0.693147180559945, 2.07944154167984, 2.89037175789616, 3.46573590279973, 3.91202300542815]]
示例8：机器学习 - 特征工程
对数值型特征进行对数变换，减少数值范围差异
features = [[1000, 2.5, 0.001], [500, 1.8, 0.005]] # 不同尺度的特征
print("对数变换特征:", reallog(features))
#对数变换特征: [[6.90775527898214, 0.916290731874155, -6.90775527898214],
[6.21460809842219, 0.587786664902119, -5.29831736654804]]
示例9：符号计算示例
symbolic_matrix = [[x, 2 * y], [x ** 2, y + 1]]
print("符号矩阵对数:", reallog(symbolic_matrix))
#符号矩阵对数: [[log(x), log(2*y)],
[log(x**2), log(y + 1)]]
示例10：标量计算
单个数值的对数计算
print("e的自然对数:", reallog(2.718281828459045))
print("10的自然对数:", reallog(10))
print("1的自然对数:", reallog(1))
#e的自然对数: 1.0
#10的自然对数: 2.30258509299405
#1的自然对数: 0
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def real_log_array(input_str):
"""
计算非负实数数组的自然对数（对标MATLAB的reallog函数）
参数:
input_str: 输入可以是以下类型：
- 非负实数/数值的字符串表示（如"5"）
- 矩阵的字符串表示（如"[[1, 2], [3, 4]]"）
- 直接传入Python列表或SymPy矩阵
返回:
sp.Matrix/sp.Expr: 自然对数计算结果
str: 错误时返回错误信息
注意:
- 输入必须全为非负实数，负数或复数将引发错误
- 符号变量需明确声明为非负实数（如Symbol('x', nonnegative=True)）
- 支持符号运算但需满足非负实数条件
"""
try:
expr = sp.sympify(input_str)
result = None
def check_nonnegative(e):
"""递归检查元素的非负实数属性"""
return bool(e.is_real and e.is_nonnegative) # 确保返回布尔值
if expr.is_number or expr.is_symbol:
# 检查标量是否为非负实数
if not check_nonnegative(expr):
return f"错误：输入 {expr} 包含非负实数以外的元素"
result = sp.log(expr).evalf()
else:
return f"输入类型错误: {type(expr)}"
return result
except Exception as e:
return f"错误:{e}"
# 示例代码
if __name__ == "__main__":
# 示例1：标量输入
print("示例1:", real_log_array("2.718"))
# 0.999896315728952
仅实数输出的数组幂
Z = realpow(X,Y) 将数组X的每个元素提升到数组Y中的对应元素的幂.数组X和Y的大小必须相同.realpow的范围是所有实数集合,即输出数组Z的所有元素必须为实数.
示例1：几何学 - 面积和体积计算
计算不同形状的面积和体积
side_lengths = [2, 3, 4, 5]
print("正方形面积 (指数2):", realpow(side_lengths, 2))
print("立方体体积 (指数3):", realpow(side_lengths, 3))
#正方形面积 (指数2): [ 4. 9. 16. 25.]
#立方体体积 (指数3): [ 8. 27. 64. 125.]
示例2：物理学 - 运动学计算
自由落体距离：h = 0.5 * g * t^2
times = [1, 2, 3, 4] # 时间(s)
g = 9.8 # 重力加速度
print("下落距离 (t^2):", realpow(times, 2)
print("实际下落距离:", 0.5 * g * np.array(realpow(times, 2))
#下落距离 (t^2): [ 1. 4. 9. 16.]
#实际下落距离: [ 4.9 19.6 44.1 78.4]
示例3：金融学 - 复利计算
复利公式：A = P * (1 + r)^t
principal = 1000
rates = [0.05, 0.08, 0.1] # 年利率
years = [1, 5, 10] # 投资年限
print("增长因子 ((1+r)^t):", realpow([1.05, 1.08, 1.1], [1, 5, 10]))
#增长因子 ((1+r)^t): [1.05 1.46932808 2.59374246]
示例4：工程学 - 材料强度计算
梁的弯曲强度与尺寸的幂次关系
beam_dimensions = [[2, 3], [4, 5], [6, 7]] # 梁的尺寸
power_exponents = [2, 4] # 不同情况下的幂指数
print("强度相关计算:", realpow(beam_dimensions, power_exponents))
#强度相关计算:
#[[4,81]
[16,625]
[36,2401]]
示例5：生物学 - 种群增长模型
种群增长：N(t) = N0 * r^t
growth_rates = [1.2, 1.5, 2.0] # 增长率
time_periods = [1, 2, 3] # 时间周期
print("增长倍数 (r^t):", realpow(growth_rates, time_periods))
#增长倍数 (r^t): [1.2 2.25 8. ]
示例6：声学 - 声压级计算
声压级：L = 20 * log10(P/P0)，反推压力比
pressure_ratios = [1, 10, 100, 1000]
print("声压级计算中的幂运算:", realpow(10, [0.05, 0.5, 1, 1.5]))
#声压级计算中的幂运算: [ 1.12201845 3.16227766 10. 31.6227766 ]
示例7：广播运算示例
测试NumPy广播功能
vector = [1, 2, 3, 4]
scalar_power = 2
matrix = [[1, 2], [3, 4]]
vector_powers = [1, 2]
print("向量^标量:", realpow(vector, scalar_power))
#向量^标量: [ 1. 4. 9. 16.]
print("矩阵^向量:", realpow(matrix, vector_powers))
#矩阵^向量:
#[[ 1. 4.]
[ 3. 16.]]
示例8：信号处理 - 滤波器设计
巴特沃斯滤波器响应：|H(ω)|² = 1 / (1 + (ω/ωc)^(2n))
omega_values = [0.1, 0.5, 1.0] # 频率比
filter_orders = [2, 4, 6] # 滤波器阶数
print("频率项的幂:", realpow(omega_values, [4, 8, 12])) # 2n次幂
#频率项的幂: [1.00000e-04 3.90625e-03 1.00000e+00]
示例9：统计学 - 方差计算
标准差与方差的关系：σ² = (x - μ)^2
deviations = [-2, -1, 0, 1, 2] # 偏差
print("方差项 (偏差^2):", realpow(deviations, 2))
#方差项 (偏差^2): [4. 1. 0. 1. 4.]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def real_pow_array(input_str):
"""
计算实数安全的逐元素幂运算（对标MATLAB的realpow函数）
参数:
input_str: 输入可以是以下类型：
- 包含两个元素的元组（X, Y）
- 元组的字符串表示（如"([1, 2], [3, 4])"）
- 直接传入Python列表/SymPy矩阵的元组
返回:
sp.Matrix: 计算结果矩阵
str: 错误时返回错误信息
注意:
- 要求X和Y必须同维度
- X中元素为负数时，对应Y元素必须为整数
- 零的负数次幂会触发错误
- 符号变量需明确声明属性（如非负性、整数性）
"""
try:
# 将输入转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
def pow_element_sym(a, b):
"""元素级安全性检查"""
# 检查零的负数次幂
if a.is_zero and b.is_negative:
return sp.nan
# 检查负数底数与非整数指数
if a.is_negative:
if not (b.is_integer or b.is_Integer):
return sp.nan
return sp.Pow(a, b)
def pow_element_num(a, b):
"""
数值标量安全性检查的幂运算
参数：
a, b: 数值标量（整数或浮点数）
返回：
float: 计算结果，若存在安全问题则返回NaN
安全性检查：
1. 零的负数次幂返回NaN
2. 负数底数与非整数指数返回NaN
3. 正确处理0^0情况（定义为1）
示例：
>>> eval_pow_element(4, 0.5) # 2.0
>>> eval_pow_element(-4, 0.5) # nan
>>> eval_pow_element(0, -1) # nan
>>> eval_pow_element(8, 1/3) # 2.0
>>> eval_pow_element(0, 0) # 1.0
"""
# 安全性检查1：零的负数次幂
if a == 0 and b < 0:
return float('nan')
# 安全性检查2：负数底数与非整数指数
if a < 0:
# 检查b是否为整数（考虑浮点精度）
if abs(b - round(b)) > 1e-10: # 容差处理浮点精度问题
return float('nan')
# 特殊处理：0^0定义为1
if a == 0 and b == 0:
return 1.0
# 安全计算幂
return a ** b
# 主处理逻辑
if isinstance(expr, tuple) and len(expr) == 2:
# 计算结果
if all(e.is_number for e in expr):
params = tuple(float(e) for e in expr)
result = pow_element_num(*params)
elif any(e.free_symbols for e in expr):
result = pow_element_sym(*expr)
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:")
print(real_pow_array("4,0.5"))
# 2.0
# 示例2：符号运算合法
x = sp.Symbol('x')
y = sp.Symbol('y')
print("示例2:", real_pow_array("x,y"))
# x**y
非负实数数组的平方根
Y = realsqrt(X) 返回数组X的每个元素的平方根. Y的大小与X的大小相同.
如果您希望负的复数返回复数结果而不是错误消息,请改用sqrt.
X - 输入矩阵,标量,向量,矩阵
示例1：几何学 - 距离计算
计算点到原点的距离（勾股定理）
squared_distances = [4, 9, 16, 25] # 平方距离
print("实际距离:", realsqrt(squared_distances))
#实际距离: [2.0, 3.0, 4.0, 5.0]
示例2：物理学 - 速度计算
动能公式：E = 1/2 * m * v^2，求速度 v = sqrt(2E/m)
kinetic_energies = [50, 100, 200] # 动能(J)
mass = 2 # 质量(kg)
print("速度:", realsqrt([2 * e / mass for e in kinetic_energies]))
#速度: [7.07106781186548, 10.0, 14.142135623731]
示例3：统计学 - 标准差计算
方差开方得到标准差
variances = [4, 9, 16, 25] # 方差
print("标准差:", realsqrt(variances))
#标准差: [2.0, 3.0, 4.0, 5.0]
示例4：金融学 - 波动率计算
波动率是收益率方差的开方
return_variances = [0.04, 0.09, 0.16] # 收益率方差
print("年化波动率:", realsqrt(return_variances))
#年化波动率: [0.2, 0.3, 0.4]
示例5：工程学 - 应力计算
应力与应变的平方根关系
strain_energy_densities = [100, 400, 900] # 应变能密度
print("相关应力项:", realsqrt(strain_energy_densities))
#相关应力项: [10.0, 20.0, 30.0]
示例6：信号处理 - RMS计算
计算信号的RMS（均方根）值
signal_squared_values = [[1, 4, 9], [16, 25, 36]] # 信号平方值
print("RMS值:", realsqrt(signal_squared_values))
#RMS值: [[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0]]
示例7：图像处理 - 亮度调整
对像素值的平方根变换，用于对比度增强
pixel_values = [[64, 144, 196], [100, 169, 225]] # 像素值（0-255）
print("平方根变换后:", realsqrt(pixel_values))
#平方根变换后: [[8.0, 12.0, 14.0],
[10.0, 13.0, 15.0]]
示例8：声学 - 声压计算
声强与声压的平方关系：I ∝ p²
sound_intensities = [1e-12, 1e-6, 1e-2] # 声强(W/m²)
print("相对声压:", realsqrt(sound_intensities))
#相对声压: [1.0e-6, 0.001, 0.1]
示例9：化学 - 浓度计算
反应速率与浓度的平方根关系
reaction_rates_squared = [0.25, 1.0, 2.25] # 反应速率的平方
print("反应速率:", realsqrt(reaction_rates_squared))
#反应速率: [0.5, 1.0, 1.5]
示例10：生物学 - 生长模型
生物体尺寸与表面积的平方根关系
surface_areas = [100, 400, 900] # 表面积
print("特征长度:", realsqrt(surface_areas))
#特征长度: [10.0, 20.0, 30.0]
示例11：符号计算示例
使用符号变量进行计算
symbolic_matrix = [[x, 4*y], [x**2, y+1]]
print("符号矩阵平方根:", realsqrt(symbolic_matrix))
#符号矩阵平方根: [[sqrt(x), 2*sqrt(y)],
[x, sqrt(y + 1)]]
示例12：机器学习 - 特征缩放
对特征进行平方根变换，减少偏度
skewed_features = [[1, 4, 9, 16], [25, 36, 49, 64]]
print("平方根变换后:", realsqrt(skewed_features))
#平方根变换后: [[1.0, 2.0, 3.0, 4.0],
[5.0, 6.0, 7.0, 8.0]]
示例13：质量控制 - 过程能力分析
过程能力指数中的标准差计算
variance_measurements = [0.01, 0.04, 0.09, 0.16]
print("标准差:", realsqrt(variance_measurements))
#标准差: [0.1, 0.2, 0.3, 0.4]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def real_sqrt_array(input_str):
"""
计算非负实数数组的平方根，对标MATLAB的realsqrt函数
参数:
input_str: 输入表达式（SymPy矩阵/列表/矩阵字符串）
返回:
sp.Matrix: 数值平方根矩阵（元素为浮点数）
str: 错误信息（若输入非法或存在负数元素）
示例:
>>> real_sqrt_array("[1, 4, 9]")
Matrix([[1.0], [2.0], [3.0]])
>>> real_sqrt_array("Matrix([[0, 2], [3, 4.0]])")
Matrix([[0.0, 1.4142135623730951], [1.7320508075688772, 2.0]])
>>> real_sqrt_array("[[-1, 4], [2, 3]]")
'错误：元素 -1 是负数'
>>> real_sqrt_array("[x, 4]")
'错误：元素 x 无法转换为数值'
"""
try:
expr = sp.sympify(input_str)
sp_matrix = sp.Matrix(expr) if isinstance(expr, list) else None
# 检查输入是否为有效矩阵
if sp_matrix is None:
return f"输入错误: {input_str} 无法转换为矩阵"
# 检查所有元素是否为非负实数
for element in sp_matrix:
try:
num = element.evalf() # 尝试转换为数值
if not num.is_real: # 确保元素为实数
return f"错误：元素 {element} 不是实数"
if num < 0: # 确保元素非负
return f"错误：元素 {element} 是负数"
float(num) # 强制转换为浮点数验证数值有效性
except:
return f"错误：元素 {element} 无法转换为数值"
# 计算平方根并数值化结果
return sp_matrix.applyfunc(sp.sqrt).evalf()
except Exception as e:
return f"运行时错误: {e}"
# 合法数值输入
print(real_sqrt_array("[1, 4, 9]"))
# Matrix([[1.00000000000000],
# [2.00000000000000],
# [3.00000000000000]])
print(real_sqrt_array("[[0, 2], [3, 4.0]]"))
# Matrix([[0, 1.41421356237310],
# [1.73205080756888, 2.00000000000000]])
参数区域图
RegionPlot(expression1, expression2)通过参数方程定义了一个二维区域，并可视化该区域在笛卡尔平面中的形状。
1. 机械振动分析 - 李萨如图形区域
可视化机械系统在特定频率比下的振动轨迹包络区域
RegionPlot(2*sin(3*theta+phi),1.5*sin(2*theta),phi=[0,@pi/4],theta=[0,2@pi])
2. 电磁场分析 - 偶极子辐射场
天线辐射方向图的3D可视化
RegionPlot(1/(1-0.5*cos(phi))*sin(theta)*cos(phi),1/(1-0.5*cos(phi))*sin(theta)*sin(phi),theta=[0,@pi],phi=[0,2@pi])
3. 流体力学 - 卡门涡街区域
圆柱绕流中涡旋脱落的周期性区域
RegionPlot(2t+0.3*sin(2*@pi*0.4*t)*exp(-0.2y),y,t=[0,5],y=[-1,1])
4. 材料科学 - 晶体滑移系统
六方晶系材料在特定取向下的滑移系统激活区域
RegionPlot(r*cos(3*theta),r*sin(3*theta),r=[0.8,1.2],theta=[0,@pi/3])
5. 化学工程 - 反应器停留时间分布
连续搅拌釜反应器系统的停留时间分布区域
RegionPlot(t,(t^(n-1)*exp(-t*n))/(n-1)!,t=[0,3],n=[1,10])
复变函数的实部和虚部关系图
ReImPlot(expression)展示了一个复变函数在实轴这个“切片”上的行为，让我们能像观察一元实函数一样，直观地看到它的实部和虚部如何随实数x变化。
1. 平方根函数
标记出了函数的分支切割位置。在切割线上，函数发生了“跳跃”，实部和虚部表现出不连续的特性。
这对于物理和工程中处理多值函数（如电磁场、流体力学）至关重要，因为它指导我们如何选择正确的分支以避免得到非物理的结果。
ReImPlot(sqrt((x^2-1)*(x^2-4)),x=[-3,3])
2. 倒数函数
在 x=0 点，函数趋向于无穷，这是一个极点。图清晰地显示了极点的位置和性质。
在电磁学中，这可以表示一个点电荷的电位在电荷所在位置是奇异的。
ReImPlot(1/x,x=[-3,3])
3. 指数函数
这展示了指数函数在实轴上是纯实数的、且没有任何奇点的良好性质。
它常用于描述增长和衰减过程，如人口增长、放射性衰变、电路暂态响应等。其导数仍是自身，这一特性在求解微分方程时极为有用。
ReImPlot(exp(x),x=[-2,2])
4. 对数函数
这个图完美地展示了对数函数的分支切割。在图上，当 x 从正方向接近0时，ln(x) 趋向 -∞。当 x 穿越0到负实轴时，虚部发生了一个π的跳跃，这正体现了函数在分支切割上的不连续性。
对数函数在复分析、数论（例如素数定理）和信号处理（复倒谱分析）中应用广泛。
ReImPlot(ln(abs(x)),x=[-3,3])
5. 三角函数
这验证了复三角函数是实三角函数的自然推广，在实轴上完全退化为我们熟悉的样子。
它在描述波动现象（声波、光波、交流电）时是核心工具。
ReImPlot(sin(x),x=[-2@pi,2@pi])
叠加复变函数的实部和虚部关系图
ReImPlot(expression1, expression2)在同一个坐标系中直观地比较多个函数的性质，揭示它们之间深刻的内在联系。
直接比较不同函数的行为, 发现函数之间的恒等关系、对称性或差异, 同时观察多个函数共有的特性（如奇点、分支点）。
1. 三角恒等式
在 [-1, 1] 区间，两条实部曲线关于 π/4 对称，像镜子一样，它们的和是常数。
超出这个区间，虽然函数变得复杂，但它们的实部之和仍然保持为常数 π/2。这验证了恒等式的强大和普适性。
ReImPlot(asin(x),acos(x),x=[-4,4])
2. 指数函数与双曲函数
悬链线问题：描述悬挂在两点的绳子形状。
狭义相对论：洛伦兹变换因子 γ 和快度（rapidity）的关系就用双曲函数描述，其本质是指数函数在闵可夫斯基时空的体现。
ReImPlot(exp(x),sinh(x),cosh(x),x=[-2,2])
3. 函数与其导数
导数即斜率：sec(x)^2 的曲线完美地描绘了 tan(x) 曲线上每一点的瞬时变化率（斜率）。
奇点行为：当 x 接近 ±π/2 时，tan(x) 和 sec(x)^2 都趋向于无穷，这直观地展示了函数及其导数在奇点附近的行为。
在光学中，正切函数可以描述某些边界条件，其导数则与场的强度相关。
ReImPlot(tan(x),sec(x)^2,x=[-1.5,1.5])
4. 恒等函数与平方函数
f(z)=z 只是一个旋转和平移。
f(z)=z^2 会将角度加倍（例如，第一象限会被映射到上半平面）。
ReImPlot(x,x^2,x=[-2,2])
复变函数的实部和虚部关系三维图
ReImPlot3D(expression1) 不再将自变量限制在实轴上，而是允许其在复平面的一个区域上自由变化。它通过两个三维曲面来完整地表示一个复变函数.
ReImPlot3D 将复变函数从抽象的公式变成了可视化，理解函数奇点（极点、分支点）、全局特性（解析性、周期性）和特殊行为（如ζ函数的零点）
1. 平方根函数
实部曲面 (u(x,y))： 形状像一个拧毛巾的曲面。
虚部曲面 (v(x,y))： 形状像一个螺旋桨或双曲抛物面。
在物理学中，sqrt(z) 经常出现在涉及平方关系的领域，例如波矢 k = sqrt(ω² - ω_c²) 在波导理论中的计算，正确处理其分支至关重要。
ReImPlot3D(sqrt(x),x=[-3-3@i,3+3@i])
2. 倒数函数
实部曲面：u(x,y) = x/(x²+y²)。形状像一个无限高的漏斗或尖峰。
虚部曲面：v(x,y) = -y/(x²+y²)。
极点可视化：这个图形将极点 (Pole) 的概念展现得淋漓尽致。
电磁场应用：在电磁学中，1/z 可以表示一个二维点电荷的复电位。
ReImPlot3D(1/x,x=[-3-3@i,3+3@i])
3. 指数函数
实部曲面：u(x,y) = e^x * cos(y)。这是一个随着 x 增大幅度不断增大的周期振荡波形。
虚部曲面：v(x,y) = e^x * sin(y)。同样是一个振幅增长的周期波形，但与实部曲面存在 π/2 的相位差。
调和函数：它的实部和虚部本身都是调和函数（满足拉普拉斯方程 ∇²u=0）。
信号处理：它描述了复指数信号的实部和虚部，是傅里叶分析和信号处理的基础。
ReImPlot3D(exp(x),x=[-2-2@i,2+2@i])
4. 对数函数
实部曲面： 这个曲面是一个关于原点对称的曲面，形状像一个“漏斗”或“火山口”。
虚部曲面： 虚部曲面是复平面 (x, y)(x,y) 本身（即高度始终为0的平面）。
电磁学：在二维静电学中，\ln(|z|)ln(∣z∣) 表示一个无限长线电荷的电势（泊松方程的解）。
流体力学：描述不可压缩流体的源或汇的 velocity potential（速度势），其中流体从原点流出或流入。
ReImPlot3D(ln(abs(x)),x=[-3-2@i,3+3@i])
5. 三角函数
实部曲面： u(x, y) = sin(x)*cosh(y) 整体形状像一系列“振荡的山脊”，在 y=0y=0 附近振幅较小，远离时振幅变大。
虚部曲面： 整体形状像一种“扭曲的 saddle 形状”，在 y=0y=0 处为0，随 |y|∣y∣ 增大而振幅增大。
ReImPlot3D(sin(x),x=[-2@pi-2@pi*@i,2@pi+2@pi*@i])
多曲面复变函数的实部和虚部关系三维图
ReImPlot3D(expression1, expression2, ...)用于同时比较多个复变函数在复平面上的行为。这种可视化方法特别有助于理解函数之间的关系、对称性和变换特性。
特别适合研究函数族、变换对和特殊函数之间的关系。
1. 多值函数，通过分支切割定义其主值
实部曲面: Re(asin(z)) 曲面, Re(acos(z)) 曲面, 两个实部曲面之间的关系反映了恒等式 Re(asin(z)) + Re(acos(z)) = π/2
虚部曲面: Im(asin(z)) 曲面, Im(acos(z)) 曲面, 在分支切割线上有不连续性
三维可视化直观展示了 asin(z) + acos(z) = π/2 在复平面上的表现, 清晰地显示了分支切割的位置和性质
ReImPlot3D(asin(x),acos(x),x=[-4-4@i,4+4@i])
2. 指数函数与对数函数
实际意义：
展示指数函数和对数函数作为互逆运算的关系
揭示指数函数的周期性和对数函数的多值性之间的关系
曲面特征：
exp(z) 的实部和虚部都是周期振荡的曲面
log(z) 的实部是径向对称的曲面，虚部是幅角函数（螺旋状曲面）
ReImPlot3D(exp(z),log(z),z=[-2-2@i,2+2@i])
3. 双曲函数与三角函数
实际意义：
展示双曲函数与三角函数之间的深刻联系：sinh(iz) = i sin(z) 和 cosh(iz) = cos(z)
在物理学中，这对理解振动问题（三角函数）和悬链线问题（双曲函数）之间的关系很重要
在相对论中，双曲函数描述快度，而三角函数描述旋转，这种对应关系有重要意义
曲面特征：
三角函数的实部和虚部是周期振荡的曲面
双曲函数的实部和虚部是指数增长/衰减的曲面
可以观察到通过虚数变换 z → iz 时，双曲函数曲面如何转换为三角函数曲面
ReImPlot3D(sinh(z),cosh(z),sin(z),cos(z),z=[-2-2@i,2+2@i])
4. 特殊函数
实际意义：
Gamma 函数是阶乘在复平面上的推广，在数论、概率论和物理学中有广泛应用
同时显示 Gamma 函数及其倒数有助于理解函数的极点和零点分布
在量子场论中，Gamma 函数出现在费曼积分计算中
曲面特征：
Gamma(z) 的曲面在负整数处有极点（趋向无穷）
1/Gamma(z) 的曲面在负整数处有零点（趋向零）
可以清晰地看到函数在正实轴上的增长行为以及在其他区域的振荡行为
ReImPlot3D(gamma(z),1/gamma(z), z=[-4-4@i,4+4@i])
除后的余数
r = rem(a,b)返回a除以b后的余数,其中a是被除数,b是除数.此函数通常称为求余运算,表达式为r=a-b*fix(a/b).
rem 函数遵从rem(a,0)是NaN的约定.
a — 被除数,标量,向量,矩阵
b — 除数,标量,向量,矩阵
示例1：时间计算 - 小时分钟转换
将总分钟数转换为小时和剩余分钟
total_minutes = [65, 130, 245, 380] # 总分钟数
hours_divisor = 60
print("剩余分钟数:", rem(total_minutes, hours_divisor))
#剩余分钟数: [ 5. 10. 5. 20.]
示例2：日历计算 - 星期几计算
计算某天是星期几（假设第1天是星期一）
days = [1, 8, 15, 22, 29] # 日期
week_days = 7
print("星期几 (0=周日, 1=周一, ...):", rem(days, week_days))
#星期几 (0=周日, 1=周一, ...): [1. 1. 1. 1. 1.]
示例3：工程学 - 角度归一化
将角度归一化到0-360度范围内
angles = [45, 370, -90, 720] # 角度
full_circle = 360
print("归一化角度:", rem(angles, full_circle))
#归一化角度: [ 45. 10. -90. 0.]
示例4：计算机科学 - 哈希函数
简单的哈希函数使用余数操作
keys = [123, 456, 789, 101112] # 键值
hash_table_size = 100
print("哈希位置:", rem(keys, hash_table_size))
#哈希位置: [23. 56. 89. 12.]
示例5：金融学 - 零钱计算
计算需要多少张整钞和剩余零钱
amounts = [123.45, 67.89, 256.78] # 金额
bill_denomination = 100 # 钞票面额
print("剩余零钱:", rem(amounts, bill_denomination))
#剩余零钱: [23.45 67.89 56.78]
示例6：信号处理 - 相位计算
计算信号的相位余数
phases = [3.5 * @pi, 5.2 * @pi, -2.1 * @pi] # 相位
two_pi = 2 * @pi
print("归一化相位:", rem(phases, two_pi))
#归一化相位: [ 4.71238898 3.76991118 -0.31415927]
示例7：数学 - 奇偶性判断
使用余数判断数字的奇偶性
numbers = [2, 3, 4, 5, 6, 7]
divisor = 2
print("余数 (0=偶数, 1=奇数):", rem(numbers, divisor))
#余数 (0=偶数, 1=奇数): [0. 1. 0. 1. 0. 1.]
示例8：物理学 - 周期性运动
计算简谐振动的相位余数
positions = [0.5, 1.2, 2.8, 3.5] # 位置
wavelength = 2.0 # 波长
print("相位余数:", rem(positions, wavelength))
#相位余数: [0.5 1.2 0.8 1.5]
示例9：游戏开发 - 地图坐标
计算游戏地图中的相对坐标（循环地图）
player_positions = [[15, 25], [48, 62], [105, 89]] # 玩家位置
map_size = [100, 100] # 地图尺寸
print("相对坐标:", rem(player_positions, map_size))
#相对坐标:
#[[15. 25.]
[48. 62.]
[ 5. 89.]]
示例10：密码学 - 模运算
简单的模运算应用
plaintext = [10, 23, 45, 67] # 明文
modulus = 26 # 模数（对应26个字母）
print("密文位置:", rem(plaintext, modulus))
#密文位置: [10. 23. 19. 15.]
示例11：矩阵运算
矩阵与标量的余数运算
matrix_a = [[10, 21], [32, 43]]
divisor = 10
print("矩阵余数:", rem(matrix_a, divisor))
#矩阵余数:
#[[0. 1.]
[2. 3.]]
示例12：广播运算
vector = [11, 22, 33, 44]
divisors = [3, 4, 5, 6]
print("广播余数:", rem(vector, divisors))
#广播余数: [2. 2. 3. 2.]
示例13：负数和浮点数
测试负数和浮点数的余数计算
mixed_numbers = [-10, -7.5, 8.3, 15.7]
divisor = 3
print("余数:", rem(mixed_numbers, divisor))
#余数: [-1. -1.5 2.3 0.7]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def remainder_after_division(input_expr):
"""
计算除后的余数，对标MATLAB的rem函数
参数:
input_expr: 输入表达式（元组/列表/字符串），格式为(被除数, 除数)
返回:
sp.Matrix 或 sp.Expr: 计算结果（支持符号变量）
str: 错误信息
示例:
>>> remainder_after_division((5, 3))
2
>>> remainder_after_division("(Matrix([5, -5, 5]), 3)")
Matrix([[2, -2, 2]])
>>> remainder_after_division("([x, 7], 3)")
Matrix([[Mod(x, 3)], [1]])
>>> remainder_after_division("(2, 0)")
'错误：除数不能为零'
"""
try:
expr = sp.sympify(input_expr)
error = False
result = None
def remainder_dividend_sym(a, b):
quotient = sp.sign(a) * sp.floor(sp.Abs(a) / sp.Abs(b))
return a - b * quotient
def remainder_dividend_num(a, b):
quotient = np.fix(a / b)
return a - b * quotient
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = remainder_dividend_num(*params)
elif any(e.free_symbols for e in expr):
result = remainder_dividend_sym(*expr)
else:
error = True
return result if not error else f"输入错误: {input_expr}"
except Exception as e:
return f"计算错误：{str(e)}"
# 标量与标量运算
print(remainder_after_division("-5.5, 3"))
# -2.5
# 符号与标量运算
print(remainder_after_division("x, 3"))
# x - 3*floor(Abs(x)/3)*sign(x)
重构数组
B = reshape(A,sz) 使用大小向量sz重构A以定义size(B).
A - 输入数组,向量,矩阵
sz - 输出大小, 由整数组成的行向量
示例1：图像处理 - 图像展平与重构
将图像像素矩阵展平为一维向量，然后重构
image_pixels = [[100, 150, 200],
[120, 180, 220],
[90, 130, 170]]
展平为9x1向量:
print(reshape(image_pixels, [9, 1]))
#[[100],
[120],
[90],
[150],
[180],
[130],
[200],
[220],
[170]]
重构为1x9行向量:
print(reshape(image_pixels, [1, 9]))
#[[100, 120, 90, 150, 180, 130, 200, 220, 170]]
示例2：机器学习 - 特征向量重塑
将特征向量重塑为不同的维度用于不同模型
feature_vector = list(range(1, 13)) # 12个特征
重塑为3x4特征矩阵:
print(reshape([feature_vector], [3, 4]))
#[[1, 4, 7, 10],
[2, 5, 8, 11],
[3, 6, 9, 12]]
重塑为4x3特征矩阵:
print(reshape([feature_vector}], [4, 3]))
#[[1, 5, 9],
[2, 6, 10],
[3, 7, 11],
[4, 8, 12]]
重塑为2x6特征矩阵:
print(reshape([feature_vector], [2, 6]))
#[[1, 3, 5, 7, 9, 11],
[2, 4, 6, 8, 10, 12]]
示例3：信号处理 - 时频转换
将时域信号重塑为频域分析格式
time_signal = [0.1, 0.5, 0.8, 0.3, 0.9, 0.2, 0.7, 0.4]
重塑为2x4矩阵进行分段分析:
print(reshape([time_signal], [2, 4]))
#[[0.1, 0.8, 0.9, 0.7],
[0.5, 0.3, 0.2, 0.4]]
重塑为4x2矩阵进行频域分析:
print(reshape([time_signal], [4, 2]))
#[[0.1, 0.9],
[0.5, 0.2],
[0.8, 0.7],
[0.3, 0.4]]
示例4：数值计算 - 矩阵分块
将大矩阵分块处理
large_matrix = [[i + j * 4 for i in range(1, 5)] for j in range(4)]
重塑为2x8矩阵进行分块计算:
print(reshape(large_matrix, [2, 8]))
#[[1, 9, 2, 10, 3, 11, 4, 12],
[5, 13, 6, 14, 7, 15, 8, 16]]
示例5：数据科学 - 数据表转换
将数据表重塑为不同的分析格式
sales_data = [[100, 150, 200], # 季度销售数据
[120, 180, 220],
[90, 130, 170],
[110, 160, 210]]
重塑为3x4矩阵 (3产品×4季度):
print(reshape(sales_data, [3, 4]))
#[[100, 110, 130, 220],
[120, 150, 160, 170],
[90, 180, 200, 210]]
重塑为6x2矩阵进行趋势分析:
print(reshape(sales_data, [6, 2]))
#[[100, 130],
[120, 160],
[90, 200],
[110, 220],
[150, 170],
[180, 210]]
示例6：物理学 - 坐标变换
将3D坐标重塑为不同的表示形式
coordinates_3d = [[1, 2, 3], # 点1
[4, 5, 6], # 点2
[7, 8, 9], # 点3
[10, 11, 12]] # 点4
重塑为3x4矩阵 (x,y,z坐标分别排列):
print(reshape(coordinates_3d, [3, 4]))
#[[1, 10, 8, 6],
[4, 2, 11, 9],
[7, 5, 3, 12]]
重塑为2x6矩阵进行平面投影:
print(reshape(coordinates_3d, [2, 6]))
#[[1, 7, 2, 8, 3, 9],
[4, 10, 5, 11, 6, 12]]
示例7：金融分析 - 投资组合重塑
重塑投资组合数据用于不同分析
portfolio_returns = [[0.05, 0.03, 0.07], # 资产A
[0.02, 0.04, 0.06], # 资产B
[0.08, 0.01, 0.09], # 资产C
[0.03, 0.05, 0.04]] # 资产D
重塑为3x4矩阵进行时间序列分析:
print(reshape(portfolio_returns, [3, 4]))
#[[0.05, 0.03, 0.01, 0.06],
[0.02, 0.03, 0.05, 0.09],
[0.08, 0.04, 0.07, 0.04]]
重塑为12x1向量进行整体分析:
print(reshape(portfolio_returns, [12, 1]))
#[[0.05],
[0.02],
[0.08],
[0.03],
[0.03],
[0.04],
[0.01],
[0.05],
[0.07],
[0.06],
[0.09],
[0.04]]
示例8：工程学 - 传感器数据重组
重组传感器阵列数据
sensor_readings = [[23.5, 24.1, 23.8], # 传感器1
[22.9, 23.7, 24.2], # 传感器2
[24.3, 23.2, 23.9], # 传感器3
[23.1, 24.0, 23.6]] # 传感器4
重塑为2x6矩阵进行空间分析:
print(reshape(sensor_readings, [2, 6]))
#[[23.5, 24.3, 24.1, 23.2, 23.8, 23.9],
[22.9, 23.1, 23.7, 24.0, 24.2, 23.6]]
示例9：符号计算示例
symbolic_matrix = [[x, y], [z, x + y], [y + z, x * z]]
重塑为3x2矩阵:
print(reshape(symbolic_matrix, [3, 2]))
#[[x, y],
[z, x + y],
[y + z, x*z]]
示例10：神经网络 - 权重矩阵重塑
神经网络权重矩阵的重塑
neural_weights = list(range(1, 25)) # 24个权重
重塑为4x6隐藏层权重矩阵:
print(reshape([neural_weights], [4, 6]))
#[[1, 5, 9, 13, 17, 21],
[2, 6, 10, 14, 18, 22],
[3, 7, 11, 15, 19, 23],
[4, 8, 12, 16, 20, 24]]
重塑为8x3输出层权重矩阵:
print(reshape([neural_weights], [8, 3]))
#[[1, 9, 17],
[2, 10, 18],
[3, 11, 19],
[4, 12, 20],
[5, 13, 21],
[6, 14, 22],
[7, 15, 23],
[8, 16, 24]]
示例11：控制系统 - 状态空间矩阵
控制系统状态空间表示的重塑
state_matrix = [[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12]]
重塑为2x6控制器输入格式:
print(reshape(state_matrix, [2, 6]))
#[[1, 9, 6, 3, 11, 8],
[5, 2, 10, 7, 4, 12]]
重塑为6x2观测器格式:
print(reshape(state_matrix, [6, 2]))
#[[1, 3],
[5, 7],
[9, 11],
[2, 4],
[6, 8],
[10, 12]]
示例12：量子计算 - 量子态重塑
量子态向量的重塑
quantum_state = [0.5, 0.5, 0.5, 0.5] # 2-qubit状态
重塑为2x2密度矩阵:
print(reshape([quantum_state], [2, 2]))
#[[0.5, 0.5],
[0.5, 0.5]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def reshape_matrix(input_str):
"""
对标MATLAB的reshape函数，按列优先顺序重塑矩阵维度
参数:
input_str (str): 输入字符串，格式为 "([[矩阵元素], [新形状])"，例如 "([[1,2],[3,4]], [2,2])"
返回:
sp.Matrix or str: 成功返回重塑后的SymPy矩阵，失败返回错误信息字符串
示例:
>>> reshape_matrix("([[1,2],[3,4]], [4, 1])")
Matrix([
[1],
[3],
[2],
[4]])
>>> reshape_matrix("([[1,2,3],[4,5,6]], [3,2])")
Matrix([
[1, 5],
[4, 3],
[2, 6]])
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
# 检查输入是否为 (矩阵, 新形状) 的元组
if not (isinstance(expr, tuple) and len(expr) == 2):
return "输入错误: 输入格式应为 (矩阵, 新形状)，例如 ([[1,2],[3,4]], [2,2])"
matrix_part, shape_part = expr
# 将矩阵部分转换为SymPy矩阵
sp_matrix = sp.Matrix(matrix_part) if isinstance(matrix_part, list) else None
if sp_matrix is None:
return f"输入错误: 无法解析矩阵部分 {matrix_part}"
# 解析新形状参数
if not isinstance(shape_part, (list, tuple)):
return f"输入错误: 新形状应为列表或元组，但得到 {type(shape_part).__name__}"
if len(shape_part) != 2:
return f"输入错误: 新形状需要2个维度，但得到 {len(shape_part)} 个"
try:
new_rows, new_cols = map(int, shape_part) # 确保形状为整数
except TypeError:
return "输入错误: 新形状维度必须为整数"
# 验证元素总数匹配
original_size = sp_matrix.rows * sp_matrix.cols
new_size = new_rows * new_cols
if original_size != new_size:
return f"错误: 元素总数不匹配 ({original_size} → {new_size})"
# 按列优先顺序展开元素
elements = []
for col in range(sp_matrix.cols):
for row in range(sp_matrix.rows):
elements.append(sp_matrix[row, col])
# 重塑为临时矩阵 (new_cols行, new_rows列)
temp_matrix = sp.Matrix(elements).reshape(new_cols, new_rows)
# 转置得到最终矩阵 (new_rows行, new_cols列)
result_matrix = temp_matrix.T
return result_matrix
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1: 2x2矩阵重塑为4x1
print("示例1:")
print(reshape_matrix("([[1, 2], [3, 4]], [4, 1])"))
# Matrix([[1],
# [3],
# [2],
# [4]])
# 示例2: 2x3矩阵重塑为3x2
print("\n示例2:")
print(reshape_matrix("([[1, 2, 3], [4, 5, 6]], [3, 2])"))
# Matrix([[1, 5],
# [4, 3],
# [2, 6]])
部分分式展开
[r,p,k] = residue(b,a) 计算以如下形式展开的两个多项式之比的部分分式展开式的留数,极点和直项.
b — 分子多项式的系数,数字向量
a — 分母多项式的系数,数字向量
示例1: 二阶系统的单位阶跃响应
二阶系统传递函数:H(s) = 1/(s² + 2s + 2)
print(residue([1], [1, 2, 2]))
#(1.0, 1.0), (-1.0 - 1.0*I, -1.0 + 1.0*I), [0.0]
示例2: 带零点的系统
带零点的二阶系统:H(s) = (2s+3)/(s²+5s+6)
print(residue([2, 3], [1, 5, 6]))
#(3.0, -1.0), (-3.0, -2.0), [0.0]
示例3: RLC串联电路阻抗
RLC串联电路阻抗:Z(s) = s²/(s²+2s+2)
print(residue([1, 0, 0], [1, 2, 2]))
#(2.0*I, -2.0*I), (-1.0 - 1.0*I, -1.0 + 1.0*I), [1.0]
示例4: 并联谐振电路
并联RLC电路导纳:Y(s) = (s+2)/(s²+3s+2)
print(residue([1, 2], [1, 3, 2]))
#(1.0,), (-1.0,), [0.0]
示例5: 低通滤波器
Butterworth低通滤波器:H(s) = 1/(s² + √2s + 1)
print(residue([1], [1, 1.414, 1]))
#(0.707213578500707, 0.707213578500707), (-0.707 - 0.707213546250353*I, -0.707 + 0.707213546250353*I), [0.0]
示例6: 带通滤波器
带通滤波器传递函数:H(s) = s/(s²+s+1)
print(residue([1, 0], [1, 1, 1]))
#(-0.5 - 0.866025403784439*I, -0.5 + 0.866025403784439*I), (-0.5 - 0.866025403784439*I, -0.5 + 0.866025403784439*I), [0.0]
示例7: 单自由度振动系统
质量-弹簧-阻尼系统:G(s) = 1/(s²+0.5s+4)
print(residue([1], [1, 0.5, 4]))
#(0.25, 0.25), (-0.25 - 1.98431348329844*I, -0.25 + 1.98431348329844*I), [0.0]
示例8: 一阶热传导系统
一阶热系统:T(s) = 2/(s+3)
print(residue([2], [1, 3]))
#(2.0,), (-3.0,), [0.0]
示例9: 具有延迟补偿的系统
带补偿的热系统:G(s) = (s+1)/(s²+4s+3)
print(residue([1, 1], [1, 4, 3]))
#(1.0,), (-3.0,), [0.0]
示例10: 三阶系统
三阶控制系统:H(s) = (s+2)/(s³+6s²+11s+6)
print(residue([1, 2], [1, 6, 11, 6]))
#(-0.5, 0.5), (-3.0, -1.0), [0.0]
示例11: 具有复数极点的情况
振荡系统（复数极点）:H(s) = 1/(s²+1)
print(residue([1], [1, 0, 1]))
#(1.0, 1.0), (-1.0*I, 1.0*I), [0.0]
示例12: 分子阶次高于分母（直项不为空）
非真分式系统:H(s) = (s²+2)/(s+1) = s-1 + 3/(s+1)
print(residue([1, 0, 2], [1, 1]))
#(3.0,), (-1.0,), [1.0, -1.0]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def remain_partial_fraction(input_str):
"""
对标MATLAB的residue函数，执行部分分式分解，返回留数、极点和直项
参数:
input_str (str): 输入字符串，格式为 "([分子系数列表], [分母系数列表])"，例如 "([1,2], [1,4,3])"
返回:
tuple or str: 成功时返回 (留数列表, 极点列表, 直项系数列表)，失败时返回错误信息字符串
示例:
>>> remain_partial_fraction("([1, 2], [1, 4, 3])")
([1.5, -0.5], [-3.0, -1.0], [])
>>> remain_partial_fraction("([1], [1, -2, 1])")
([0.0, 1.0], [1.0, 1.0], [])
"""
try:
x = sp.symbols('x')
expr = sp.sympify(input_str)
# 验证输入格式
if not (isinstance(expr, tuple) and len(expr) == 2):
return "输入错误：输入应为包含两个系数列表的元组，例如 ([1,2], [3,4,5])"
num_coeffs, den_coeffs = expr
# 构建分子和分母多项式
numerator = sp.Poly(num_coeffs, x)
denominator = sp.Poly(den_coeffs, x)
# 分离直项和真分式部分
quotient, remainder = sp.div(numerator, denominator)
direct_terms = [float(coeff.evalf()) for coeff in quotient.all_coeffs()]
# 对真分式进行部分分式分解
f_true = remainder / denominator.as_expr()
partial = sp.apart(f_true, x, full=True)
# 初始化结果容器
residues = []
poles = []
# 处理每个分式项
for term in partial.as_ordered_terms():
term = term.simplify()
if term == 0:
continue
if term.is_polynomial():
continue # 直项已处理
# 分离分子分母并因式分解
numer, denom = term.as_numer_denom()
denom_factorized = sp.factor(denom)
# 遍历所有分母因子
for factor, power in sp.factor_list(denom_factorized)[1]:
# 处理幂次因子
if factor.is_Pow:
base, exp = factor.args
else:
base = factor
exp = 1
# 处理线性因子 (x - a)
if sp.Poly(base).is_linear:
a = sp.solve(base, x)[0]
other_denom = denom_factorized / (base ** exp)
residue_expr = numer / other_denom
# 添加极点并计算对应留数
for i in range(exp):
pole = a.evalf()
residue_val = residue_expr.subs(x, a).evalf() if i == exp - 1 else 0.0
poles.append(float(pole))
residues.append(float(residue_val))
# 处理不可约二次因子（复数根）
else:
# 提取复数根
roots = sp.roots(base, multiple=True)
for root in roots:
pole = root.evalf()
residue_val = (numer / (denom_factorized / base ** exp)).subs(x, pole).evalf()
poles.append(complex(pole))
residues.append(complex(residue_val))
# MATLAB式排序：按实部升序，虚部升序
sorted_pairs = sorted(zip(poles, residues),
key=lambda x: (x[0].real if isinstance(x[0], complex) else x[0],
x[0].imag if isinstance(x[0], complex) else 0))
if sorted_pairs:
poles, residues = zip(*sorted_pairs)
return (residues, poles, direct_terms)
except Exception as e:
return f"错误：{str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1: 单极点情况
print("示例1:")
print(remain_partial_fraction("([-4,8], [1,6,8])"))
# ((-12.0, 8.0), (-4.0, -2.0), [0.0])
# 示例2: 重复极点
print("\n示例2:")
print(remain_partial_fraction("([2,0,0,1,0], [1, 0, 1])"))
# (((2-1j), (2+1j)), (-1j, 1j), [2.0, 0.0, -2.0])
# 示例3: 复数极点
print("\n示例3:")
print(remain_partial_fraction("([1], [1, 0, 1])"))
# (((1+0j), (1+0j)), (-1j, 1j), [0.0])
旋转曲面
RevolutionPlot3D(expr1,expr2)通过将平面曲线绕固定轴旋转而生成的 3D 曲面。这种曲面在数学、工程和自然界中广泛存在，具有优美的对称性和实用性。
1. 圆环面 (Torus)
几何意义：将小圆绕大圆旋转
实际应用： 甜甜圈形状, 管道连接件, 粒子加速器的环形轨道
RevolutionPlot3D(3+cos(t),sin(t),t=[0,2@pi])
2. 球面 (Sphere)
几何意义：半圆绕直径旋转
实际应用：球形储罐, 天文馆穹顶, 轴承滚珠
RevolutionPlot3D(cos(t),sin(t),t=[0,2@pi]
3. 圆锥面 (Cone)
几何意义：直线绕 z 轴旋转
实际应用：圆锥形屋顶, 离心分离器, 交通锥标
RevolutionPlot3D(t,t=[0,5])
4. 双曲面 (Hyperboloid)
几何意义：双曲线绕轴旋转
实际应用：核电站冷却塔, 太空望远镜支撑结构, 建筑薄壳结构
RevolutionPlot3D(cosh(t),t,t=[-1,1])
5. 花瓶曲面 (Vase)
几何意义：波动曲线旋转
实际应用：陶瓷器皿设计, 建筑柱式装饰, 容器造型优化
RevolutionPlot3D(2+sin(t)+t/5,t,t=[0,4@pi])
6. 螺旋曲面 (Spiral Surface)
几何意义：螺旋线绕轴旋转
实际应用：螺旋输送机, DNA分子模型, 涡轮叶片设计
RevolutionPlot3D(t,sin(2t),t=[0,4@pi])
rewrite（expr，target）根据目标函数target重写符号表达式expr.
重写后的表达式在数学上等同于原始表达式,如果expr是一个向量或矩阵,rewrite对expr执行元素式操作.
示例 1: 将三角函数用指数函数重写
print("示例1输出:", rewrite(sin(x), exp))
#示例1输出: -I*(exp(I*x) - exp(-I*x))/2
示例 2: 将矩阵中的三角函数用指数函数重写
print("示例2输出:", rewrite([[sin(x), cos(x), tan(x)]], exp))
#示例2输出: [[-I*(exp(I*x) - exp(-I*x))/2, exp(I*x)/2 + exp(-I*x)/2, I*(-exp(I*x) + exp(-I*x))/(exp(I*x) + exp(-I*x))]]
示例 3: 将双曲函数用指数函数重写
print("示例3输出:", rewrite(sinh(x), exp))
#示例3输出: exp(x)/2 - exp(-x)/2
示例 4: 将反三角函数用对数函数重写
print("示例4输出:", rewrite(asin(x), log))
#示例4输出: -I*log(I*x + sqrt(1 - x**2))
示例 5: 将矩阵中的各种函数用指数重写
print("示例5输出:", rewrite([[sin(x), cos(x), asin(x), acos(x)]], exp))
#示例5输出: [[-I*(exp(I*x) - exp(-I*x))/2, exp(I*x)/2 + exp(-I*x)/2, asin(x), acos(x)]]
示例 6: 将三角函数用双曲函数重写
print("示例6输出:", rewrite(sin(x), sinh))
#示例6输出: sin(x)
示例 7: 将伽马函数用其他特殊函数重写
print("示例7输出:", rewrite(gamma(x), factorial))
#示例7输出: factorial(x - 1)
示例 8: 将误差函数用其他函数重写
print("示例8输出:", rewrite(erf(x), exp))
#示例8输出: erf(x)
示例 9: 将矩阵中的双曲函数用指数重写
print("示例9输出:", rewrite([[sinh(x), cosh(x), tanh(x)]], exp))
#示例9输出: [[exp(x)/2 - exp(-x)/2, exp(x)/2 + exp(-x)/2, (exp(x) - exp(-x))/(exp(x) + exp(-x))]]
示例 10: 将正切函数用正弦和余弦重写
print("示例10输出:", rewrite(tan(x), sin))
#示例10输出: 2*sin(x)**2/sin(2*x)
示例 11: 将分段函数用其他形式重写
print("示例11输出:", rewrite(Piecewise((x, x>0), (0, True)), exp))
#示例11输出: Piecewise((x, x > 0), (0, True))
示例 12: 复杂的矩阵表达式重写
print("示例12输出:", rewite([[sin(x)*cos(y), tan(x+y)], [asin(z), acos(w)]], exp))
#示例12输出: [[-I*(exp(I*x) - exp(-I*x))*(exp(I*y)/2 + exp(-I*y)/2)/2, I*(exp(I*(-x - y)) - exp(I*(x + y)))/(exp(I*(-x - y)) + exp(I*(x + y)))],
[asin(z), acos(w)]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def rewrite_another_expression(input_str):
"""
将输入表达式或矩阵中的元素用另一种函数重写。
参数:
input_str: 输入字符串，应为一个元组，格式为 (表达式或矩阵, 目标函数)。
例如："(sin(x), exp)" 或 "([[sin(x), cos(x)]], exp)"
返回:
重写后的表达式或矩阵，若输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为元组且包含两个元素
if isinstance(expr, tuple) and len(expr) == 2:
input_expr, target = expr[0], expr[1]
# 检查是否为矩阵
sp_matrix = sp.Matrix(input_expr) if isinstance(input_expr, list) else None
if sp_matrix is not None:
# 对矩阵每个元素应用 rewrite
result = sp_matrix.applyfunc(lambda x: x.rewrite(target))
else:
# 普通表达式应用 rewrite
result = input_expr.rewrite(target)
else:
error = True
return result if not error else f"输入格式错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例 1: 将 sin(x) 用 exp 重写
input_str1 = "(sin(x), exp)"
print("示例1输入:", input_str1)
print("示例1输出:", rewrite_another_expression(input_str1))
# -I*(exp(I*x) - exp(-I*x))/2
# 示例 2: 将矩阵中的三角函数用 exp 重写
input_str2 = "([[sin(x), cos(x)]], exp)"
print("\n示例2输入:", input_str2)
print("示例2输出:", rewrite_another_expression(input_str2))
# Matrix([[-I*(exp(I*x) - exp(-I*x))/2, exp(I*x)/2 + exp(-I*x)/2]])
黎曼球复变函数图
RiemannSphereComplexPlot(exp)提供复变函数在整个扩展复平面（包括无穷远）上的全景视图。
在标准复平面上，函数在无穷远的行为很难可视化（通常用箭头表示趋于无穷的方向）。在球面上，这直接对应函数值在北极N附近的行为。
许多复变变换（特别是分式线性变换/Möbius 变换）在黎曼球面上有非常简洁和优美的几何解释（旋转、缩放、平移）。
函数在无穷远点的性质（极点、本性奇点、可去奇点）可以通过分析其在北极 N 附近的行为来定义和研究。
1. 分式线性变换 (Möbius 变换) 之一
清晰展示了分式线性变换在黎曼球面上的全局等距（保角）性质，以及无穷远点与有限点的互换
RiemannSphereComplexPlot(1/z, z=[-2-2@i, 2+2@i])
2. 正弦函数是整函数（全平面解析），但超越整函数，有无穷多个零点
生动展示了本性奇点在无穷远处的复杂行为，这是标准复平面绘图难以清晰表现的（通常只能看到模的爆炸增长或振荡，但无法同时体现其“填满”整个邻域的特性）。
RiemannSphereComplexPlot(sin(z), z=[-4-4@i,4+4@i])
3. 指数函数是整函数，没有零点
直观地展示了指数函数的周期性、模随实部的指数变化规律，以及本性奇点在无穷远处的行为（不同路径趋向无穷导致函数值行为迥异：趋向北极、趋向南极、或绕赤道旋转）。
RiemannSphereComplexPlot(exp(z), z=[-4-4@i,4+4@i])
4. 这是一个有理函数（亚纯函数），在 z = ±i 处有单极点
清晰地展示了一个简单亚纯函数的全局结构：两个有限极点（映射到 N），一个在无穷远的零点（映射到 S）。
这是理解更复杂有理函数零极点分布的基础。
在系统理论（如控制论、信号处理）中，这种零极点图（通常在复平面上画）对分析系统稳定性至关重要，黎曼球面提供了一种包含无穷远点的完整视角。
RiemannSphereComplexPlot(1/(z^2+1), z=[-3-3@i,3+3@i])
黎曼球面
RiemannSphereComplexPlot3D(expression) 使得复平面上的每个点都唯一对应球面上的一个点，从而允许我们可视化复函数在无穷远点的行为，并简化某些复分析问题的处理。
1. 分式函数
展示极点和零点的对称性。
RiemannSphereComplexPlot(1/z, z=[-2-2@i, 2+2@i])
2. 简单多项式
显示二次映射的特性。
RiemannSphereComplexPlot(z^5+1/10,z=[-2-2@i,2+2@i])
3. 三角函数
在黎曼球面上显示其周期性和极点分布。
RiemannSphereComplexPlot(sin(z), z=[-4-4@i,4+4@i])
删除缺失的条目
R = rmmissing(A) 从数组或表中删除缺失的条目. 如果A是向量,则rmmissing会删除包含缺失数据的所有条目. 如果A是矩阵,则rmmissing会删除包含缺失数据的所有行.
R = rmmissing(A,dim) 指定要沿其进行运算的 A 的维度. 默认情况下, rmmissing 沿其大小不为1的第一个维度进行运算
A — 输入数据,向量,矩阵
dim — 沿其运算的维度,1或2
示例 1: 处理包含缺失值的向量
print(f"输出: {rmmissing([1, nan, 3, nan, 5])}")
#输出: [1.0, 3.0, 5.0]
示例 2: 删除包含缺失值的行
print(f"输出: {rmmissing([[1, 2, 3], [4, nan, 6], [7, 8, 9]], 1)}")
#输出: [[1, 2, 3],
[7, 8, 9]]
示例 3: 删除包含缺失值的列
print(f"输出: {rmmissing([[1, nan, 3], [4, 5, 6], [7, 8, nan]], 2)}")
#输出: [[1],
[4],
[7]]
示例 4: 传感器数据处理
时间序列数据，某些传感器读数缺失
print(f"输出: {rmmissing([[1, 23.5, 45.2, 101.3], [2, nan, 46.1, 102.1], [3, 24.1, nan, 100.8]], 1)}")
#输出: [[1, 23.5, 45.2, 101.30]]
示例 5: 实验数据预处理
不同实验条件下的测量值
print(f"输出: {rmmissing([[1.2, 2.3, nan], [3.4, nan, 5.6], [7.8, 9.1, 10.2]], 2)}")
#输出: [[1.2],
[3.4],
[7.8]]
示例 6: 调查问卷数据处理
受访者ID，年龄，收入，满意度
print(f"输出: {rmmissing([[1, 25, 50000, 4], [2, nan, 60000, 5], [3, 30, nan, 3], [4, 35, 70000, nan]], 1)}")
#输出: [[1, 25, 50000, 4]]
示例 7: 多变量时间序列
温度，湿度，气压时间序列
print(f"输出: {rmmissing([[20.5, 65, 1013], [nan, 70, 1012], [22.1, nan, 1011], [21.8, 68, nan]], 1)}")
#输出: [[20.5, 65, 1013]]
示例 8: 边缘情况 - 所有数据都缺失
print(f"输出: {rmmissing([nan, nan, nan])}")
#输出: []
示例 9: 边缘情况 - 没有缺失值
print(f"输出: {rmmissing([[1, 2, 3], [4, 5, 6]], 1)}")
#输出: [[1, 2, 3],
[4, 5, 6]]
示例 10: 大规模数据抽样
模拟从大数据集中抽取的样本
print(f"输出: {rmmissing([[1, 2.5, 3], [4, nan, 6], [7, 8.2, 9], [10, 11.1, nan]], 1)}")
#输出: [[1, 2.5, 3],
[7, 8.2, 9]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.matrices import MatrixBase
import numpy as np
def remove_missing(input_str):
"""
删除矩阵或向量中的缺失值（NaN）
参数:
input_str: 输入字符串，支持以下格式：
- 单个矩阵：例如 "[[1, nan], [3, 4]]"
- 带参数的元组：例如 "([[1, nan], [3, 4]], 1)"
参数说明：1=删除包含NaN的行，2=删除包含NaN的列
返回:
处理后的矩阵（保持矩阵结构）或向量（返回列表形式）
如果输入错误，返回错误信息字符串
"""
try:
# 将字符串中的 nan 转换为 SymPy 的 nan 符号
expr = sp.sympify(input_str, locals={'nan': sp.nan})
result = None
error = False
# 处理带参数的元组输入
if isinstance(expr, tuple):
if len(expr) != 2:
error = True
else:
matrix_part, direction = expr
sp_matrix = sp.Matrix(matrix_part) if isinstance(matrix_part, list) else None
if sp_matrix is None or direction not in (1, 2):
error = True
else:
# 转换为 NumPy 数组处理
np_array = np.array(sp_matrix, dtype=float)
if direction == 1: # 删除包含NaN的行
mask = ~np.isnan(np_array).any(axis=1)
result_array = np_array[mask]
else: # 删除包含NaN的列
mask = ~np.isnan(np_array).any(axis=0)
result_array = np_array[:, mask]
# 转换回 SymPy 矩阵
result = sp.Matrix(result_array.tolist())
# 处理单个矩阵输入
else:
sp_matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if sp_matrix is None:
error = True
else:
np_array = np.array(sp_matrix, dtype=float)
# 判断是否为向量（行数或列数为1）
if 1 in sp_matrix.shape:
# 展平并过滤NaN
flat_array = np_array.flatten()
filtered = flat_array[~np.isnan(flat_array)]
result = filtered.tolist() # 向量返回列表形式
else:
# 默认删除包含NaN的行
mask = ~np.isnan(np_array).any(axis=1)
result_array = np_array[mask]
result = sp.Matrix(result_array.tolist())
return result if not error else f"输入格式错误: {input_str}"
except Exception as e:
return f"处理错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：删除包含NaN的行
print("示例1 输入:", "([1,3,nan,6,nan])")
print("示例1 输出:", remove_missing("([1,3,nan,6,nan])"))
# [1.0, 3.0, 6.0]
# 示例2：删除包含NaN的列
print("\n示例2 输入:", "([[1, nan], [3, 4]], 2)")
print("示例2 输出:", remove_missing("([[1, nan], [3, 4]], 2)"))
# Matrix([[1.00000000000000],
# [3.00000000000000]])
# 示例3：处理向量输入
print("\n示例3 输入:", "[nan, 2, 3]")
print("示例3 输出:", remove_missing("[nan, 2, 3]"))
# [2.0, 3.0]
检测并删除数据中的离群值
B = rmoutliers(A) 在 A 的数据中检测并删除离群值。如果 A 是矩阵，则 rmoutliers 会分别检测 A 的每列中的离群值，并删除整行。
B = rmoutliers(A,method) 指定检测离群值的方法。例如，rmoutliers(A,"mean") 将 A 中与均值相差超过三倍标准差的元素定义为离群值。
median(默认) — 离群值定义为与中位数相差超过三倍换算 MAD 的元素。换算 MAD 定义为 c*median(abs(A-median(A)))，其中 c=-1/(sqrt(2)*erfcinv(3/2))。
mean - 离群值定义为与均值相差超过三倍标准差的元素。此方法比 "median" 快，但没有它可靠。
quartiles - 离群值定义为比上四分位数 (75%) 大 1.5 个四分位差以上或比下四分位数 (25%) 小 1.5 个四分位差以上的元素。当 A 中的数据不是正态分布时，此方法很有用。
grubbs - 使用针对离群值的格拉布斯检验检测离群值，即基于假设检验每次迭代删除一个离群值。此方法假设 A 中的数据呈正态分布。
A — 数值向量，矩阵。
示例 1: 传感器数据清理 - 均值法
模拟温度传感器读数，有一个异常高温值
sensor_data = [23.5, 24.1, 23.8, 25.2, 23.9, 45.6, 24.3, 23.7]
print(f"使用均值法: {rmoutliers(sensor_data, mean)}")
#使用均值法: [[23.5, 24.1, 23.8, 25.2, 23.9, 45.6, 24.3, 23.7]]
#应用场景: 工业传感器数据质量控制
示例 2: 股票收益率异常检测 - 中位数法
模拟股票日收益率，有一个极端异常值
returns = [0.02, -0.01, 0.015, -0.005, 0.03, 0.25, -0.008, 0.012, -0.002]
print(f"使用中位数法: {rmoutliers(returns, median)}")
#使用中位数法: [[0.02, -0.01, 0.015, -0.005, 0.03, -0.008, 0.012, -0.002]]
#应用场景: 金融风险管理，识别异常交易
示例 3: 学生成绩分析 - 四分位距法
模拟班级数学成绩，识别异常高/低分
math_scores = [65, 72, 68, 85, 90, 92, 88, 45, 95, 87, 82, 78, 96, 40, 89]
print(f"使用四分位距法: {rmoutliers(math_scores, quartiles)}")
#使用四分位距法: [[65, 72, 68, 85, 90, 92, 88, 45, 95, 87, 82, 78, 96, 89]]
#应用场景: 教育数据分析，识别异常成绩
示例 4: 实验测量数据 - 格拉布斯检验
模拟物理实验测量值，有一个明显离群值
measurements = [9.81, 9.79, 9.83, 9.80, 9.82, 8.95, 9.81, 9.79, 9.84]
print(f"使用格拉布斯检验: {rmoutliers(measurements, grubbs)}")
#使用格拉布斯检验: [[9.81, 9.79, 9.83, 9.8, 9.82, 8.95, 9.81, 9.79, 9.84]]
#应用场景: 科学研究，确保数据可靠性
示例 5: 二维数据 - 产品质量多维检测数据 (中位数法)
模拟产品尺寸测量（长度、宽度、高度）
quality_data = [[10.1, 5.2, 3.1], [10.0, 5.1, 3.0], [9.9, 5.3, 3.2], [12.5, 5.0, 3.1], [10.2, 5.2, 3.0]]
print(f"使用中位数法:\n{rmoutliers(quality_data, median)}")
#使用中位数法:
#[[10.1, 5.2, 3.1],
[10, 5.1, 3],
[9.9, 5.3, 3.2],
[10.2, 5.2, 3]]
#应用场景: 制造业质量控制，识别缺陷产品
示例 6: 时间序列异常点检测 - 均值法
模拟日访问量数据，包含异常峰值
traffic_data = [1250, 1320, 1280, 1290, 1350, 4500, 1310, 1270, 1330, 1295]
print(f"使用均值法: {rmoutliers(traffic_data, mean)}")
#使用均值法:
#[[1250, 1320, 1280, 1290, 1350, 4500, 1310, 1270, 1330, 1295]]
#应用场景: 网络分析，识别异常流量
示例 7: 医学数据 - 血压测量数据分析 (四分位距法)
模拟收缩压测量值
blood_pressure = [120, 118, 122, 119, 121, 185, 117, 123, 119, 124]
print(f"使用四分位距法: {rmoutliers(blood_pressure, quartiles)}")
#使用四分位距法: [[120, 118, 122, 119, 121, 117, 123, 119, 124]]
#应用场景: 医疗数据分析，确保测量准确性
示例 8: 环境监测数据 - 多变量异常检测(均值法)
模拟空气质量数据（PM2.5, CO2, 温度）
environment_data = [[35, 420, 22.5], [38, 415, 23.0], [32, 430, 22.8], [150, 410, 22.7], [36, 425, 22.9]]
print(f"使用均值法:\n{rmoutliers(environment_data, mean)}")
#使用均值法:
#[[35, 420, 22.5],
[38, 415, 23],
[32, 430, 22.8],
[150, 410, 22.7],
[36, 425, 22.9]]
#应用场景: 环境监测，识别传感器故障
示例 9: 电商交易金额分析 (格拉布斯检验)
模拟客户交易金额，识别可能的欺诈交易
transaction_amounts = [85, 120, 95, 110, 105, 2500, 90, 115, 100, 125]
print(f"使用格拉布斯检验: {rmoutliers(transaction_amounts, grubbs)}")
#使用格拉布斯检验: [[85, 120, 95, 110, 105, 2500, 90, 115, 100, 125]]
#应用场景: 金融反欺诈，检测异常交易
示例 10: 生产设备运行数据 (中位数法)
模拟机器运行温度数据
machine_temp = [72.5, 73.0, 72.8, 73.2, 72.9, 85.5, 73.1, 72.7, 73.3]
print(f"使用中位数法: {rmoutliers(machine_temp, median)}")
#使用中位数法: [[72.5, 73, 72.8, 73.2, 72.9, 73.1, 72.7, 73.3]]
#应用场景: 工业物联网，预测性维护
示例 11: 投资组合收益率分析 (四分位距法)
模拟不同资产的收益率数据
portfolio_returns = [[0.02, 0.015, 0.025], [-0.01, -0.008, -0.005], [0.03, 0.028, 0.032], [0.12, -0.02, 0.035], [0.018, 0.022, 0.02]]
print(f"使用四分位距法:\n{rmoutliers(portfolio_returns, quartiles)}")
#使用四分位距法:
#[[0.02, 0.015, 0.025],
[0.03, 0.028, 0.032],
[0.018, 0.022, 0.02]]
#应用场景: 投资管理，风险控制
示例 12: 边缘情况 - 无离群值数据
normal_data = [10, 11, 9, 10.5, 11.2, 9.8, 10.3]
print(f"使用均值法: {rmoutliers(normal_data, mean)}")
#使用均值法:
#[[10, 11, 9, 10.5, 11.2, 9.8, 10.3]]
#应用场景: 验证算法在正常数据上的行为
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.special import erfcinv
def remove_outliers_data(input_str):
"""
对标 MATLAB 的 rmoutliers 函数，检测并删除数据中的离群值
参数:
input_str: 输入数据的字符串表达式，可以是矩阵、列表或元组（矩阵, 方法）
返回:
处理后的 SymPy 矩阵或错误信息字符串
支持方法:
'mean' - 三倍标准差法（默认）
'median' - 中位数绝对偏差法
'quartiles' - 四分位距法
'grubbs' - 格拉布斯检验法
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def evaluation_rmoutliers(A, method="median"):
"""
函数功能：在输入数据A中检测并删除离群值
参数说明：
A: 输入的数据，可以是一维数组或者二维矩阵形式
method: 检测离群值的方法，支持'mean'（与均值相差超过三倍标准差）、
'median'（与中位数相差超过三倍换算MAD的元素）、
'quartiles'（基于四分位数判定）、
'grubbs'（格拉布斯检验）、
'gesd'（广义极端Student化偏差检验），默认是'mean'
返回值：
删除离群值后的结果数据
"""
if method == "mean":
if A.ndim == 1:
mean_val = np.mean(A)
std_val = np.std(A)
mask = np.abs(A - mean_val) <= 3 * std_val
return A[mask]
elif A.ndim == 2:
row_num, col_num = A.shape
keep_rows = []
for col_idx in range(col_num):
col_data = A[:, col_idx]
mean_val = np.mean(col_data)
std_val = np.std(col_data)
col_mask = np.abs(col_data - mean_val) <= 3 * std_val
keep_rows.append(col_mask)
combined_mask = np.all(keep_rows, axis=0)
return A[combined_mask, :]
elif method == "median":
if A.ndim == 1:
median_val = np.median(A)
mad = np.median(np.abs(A - median_val))
c = -1 / (np.sqrt(2) * erfcinv(3 / 2))
diff = np.abs(A - median_val)
mask = diff <= 3 * c * mad
return A[mask]
elif A.ndim == 2:
row_num, col_num = A.shape
keep_rows = []
for col_idx in range(col_num):
col_data = A[:, col_idx]
median_val = np.median(col_data)
mad = np.median(np.abs(col_data - median_val))
c = -1 / (np.sqrt(2) * erfcinv(3 / 2))
diff = np.abs(col_data - median_val)
col_mask = diff <= 2 * c * mad
keep_rows.append(col_mask)
combined_mask = np.all(keep_rows, axis=0)
return A[combined_mask, :]
elif method == "quartiles":
if A.ndim == 1:
q1 = np.percentile(A, 25)
q3 = np.percentile(A, 75)
iqr = q3 - q1
lower_bound = q1 - 1.5 * iqr
upper_bound = q3 + 1.5 * iqr
mask = (A >= lower_bound) & (A <= upper_bound)
return A[mask]
elif A.ndim == 2:
row_num, col_num = A.shape
keep_rows = []
for col_idx in range(col_num):
col_data = A[:, col_idx]
q1 = np.percentile(col_data, 25)
q3 = np.percentile(col_data, 75)
iqr = q3 - q1
lower_bound = q1 - 1.5 * iqr
upper_bound = q3 + 1.5 * iqr
col_mask = (col_data >= lower_bound) & (col_data <= upper_bound)
keep_rows.append(col_mask)
combined_mask = np.all(keep_rows, axis=0)
return A[combined_mask, :]
elif method == "grubbs":
from scipy.stats import zscore
if A.ndim == 1:
while True:
z_scores = zscore(A)
max_abs_z = np.max(np.abs(z_scores))
if max_abs_z < 3: # 可根据实际调整临界值，这里暂用3作为参考
break
index_to_remove = np.argmax(np.abs(z_scores))
A = np.delete(A, index_to_remove)
return A
elif A.ndim == 2:
row_num, col_num = A.shape
keep_rows = []
for col_idx in range(col_num):
col_data = A[:, col_idx]
while True:
z_scores = zscore(col_data)
max_abs_z = np.max(np.abs(z_scores))
if max_abs_z < 3:
break
index_to_remove = np.argmax(np.abs(z_scores))
col_data = np.delete(col_data, index_to_remove)
keep_rows.append(np.isin(A[:, col_idx], col_data))
combined_mask = np.all(keep_rows, axis=0)
return A[combined_mask, :]
elif method == "gesd":
# 这里只是个简单框架示例，实际的gesd算法较复杂，可能需进一步完善和优化
if A.ndim == 1:
# 假设这里有合适的gesd实现逻辑，暂用简单逻辑替代
return A
elif A.ndim == 2:
# 同理，对于二维情况暂简单处理
return A
else:
return None
# 解析输入参数
if isinstance(expr, tuple) and len(expr) == 2:
matrix_part = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
method_part = expr[1]
if matrix_part is None or not str(method_part) in ["mean", "median", "quartiles", "grubbs"]:
error = True
else:
# 转换为 numpy 数组
if matrix_part.shape[1] == 1:
A = np.ravel(np.array(matrix_part, dtype=float))
else:
A = np.array(matrix_part, dtype=float)
result = evaluation_rmoutliers(A, method=str(method_part))
# 无方法参数的情况
elif isinstance(expr, list):
matrix_part = sp.Matrix(expr)
if matrix_part.shape[1] == 1:
A = np.ravel(np.array(matrix_part, dtype=float))
else:
A = np.array(matrix_part, dtype=float)
result = evaluation_rmoutliers(A, method="mean") # 默认使用均值法
else:
error = True
if error:
return f"输入格式错误: {input_str}"
else:
# 转换回 SymPy 矩阵并调整形状
result = sp.Matrix(result)
if result.shape[1] == 1 and A.ndim == 1:
return result.T # 一维数据返回行向量
return result
except Exception as e:
return f"处理过程中发生错误: {str(e)}"
if __name__ == "__main__":
# 示例代码
# 示例1：一维数据均值法
print("示例1 输出:", remove_outliers_data("[1, 2, 3, 100]"))
# Matrix([[1.00000000000000, 2.00000000000000, 3.00000000000000, 100.000000000000]])
# 示例2：二维数据中位数法
data = "[[1,200],[2,3],[3,4],[5,6]]"
print("示例2 输出:\n", remove_outliers_data(f"({data}, 'median')"))
# Matrix([[2.00000000000000, 3.00000000000000],
# [3.00000000000000, 4.00000000000000],
# [5.00000000000000, 6.00000000000000]])
# 示例3：四分位距法处理二维数据
data = "[[1, 2], [3, 4], [5, 6], [100, 200]]"
print("示例3 输出:\n", remove_outliers_data(f"({data}, 'quartiles')"))
# Matrix([[1.00000000000000, 2.00000000000000],
# [3.00000000000000, 4.00000000000000],
# [5.00000000000000, 6.00000000000000]])
# 示例4：格拉布斯检验法处理一维数据
print("示例4 输出:", remove_outliers_data("([10, 12, 13, 14, 100], 'grubbs')"))
# Matrix([[10.0000000000000, 12.0000000000000, 13.0000000000000, 14.0000000000000, 100.000000000000]])
多项式根
r = roots(p) 以列向量r的形式返回p中的系数表示的多项式的根.输入p是一个包含n+1多项式系数的向量,以x^n系数开头.例如：p = [3 2 -2]表示多项式3x^2+2x−2.0系数表示方程中不存在的中间幂.
p — 多项式系数,向量
示例1. 控制系统稳定性分析
control_system = "s**3 + 5*s**2 + 8*s + 6"
print(f"特征方程: {control_system}")
roots = roots(control_system.replace('s', 'x'))
print(f"特征根: {roots}")
#特征方程: s**3 + 5*s**2 + 8*s + 6
#特征根: [-3, -1.0 - 1.0*I, -1.0 + 1.0*I]
示例2. 机械振动系统固有频率
vibration_eq = "m*s**2 + c*s + k"
代入典型参数值：m=1, c=0.5, k=2
vibration_eq_numeric = vibration_eq.replace('m', '1').replace('c', '0.5').replace('k', '2').replace('s', 'x')
print(f"特征方程: {vibration_eq_numeric}")
roots = roots(vibration_eq_numeric)
print(f"特征根: {roots}")
#特征方程: 1*x**2 + 0.5*x + 2
#特征根: [-0.25 - 1.391941091*I, -0.25 + 1.391941091*I]
示例3. RLC电路谐振分析
rlc_circuit = "L*C*s**2 + R*C*s + 1"
代入参数：L=0.1, C=0.01, R=5
rlc_numeric = rlc_circuit.replace('L', '0.1').replace('C', '0.01').replace('R', '5').replace('s', 'x')
print(f"特征方程: {rlc_numeric}")
roots = roots(rlc_numeric)
print(f"特征根: {roots}")
#特征方程: 0.1*0.01*x**2 + 5*0.01*x + 1
#特征根: [-25.0 - 19.36491673*I, -25.0 + 19.36491673*I]
示例4. 量子谐振子能级
quantum_oscillator = "x**2 - 2*x + 2"
print(f"特征方程: {quantum_oscillator}")
roots = roots(quantum_oscillator)
print(f"特征根: {roots}")
#特征方程: x**2 - 2*x + 2
#特征根: [1.0 - 1.0*I, 1.0 + 1.0*I]
示例5. 经济增长模型稳定性
economic_growth = "r**2 - 1.5*r + 0.5"
print(f"特征方程: {economic_growth.replace('r', 'x')}")
roots = roots(economic_growth.replace('r', 'x'))
print(f"特征根: {roots}")
#特征方程: x**2 - 1.5*x + 0.5
#特征根: [0.5, 1]
示例6. 数字滤波器设计
digital_filter = "z**2 - 1.8*z + 0.81"
print(f"特征方程: {digital_filter.replace('z', 'x')}")
roots = roots(digital_filter.replace('z', 'x'))
print(f"特征根: {roots}")
#特征方程: x**2 - 1.8*x + 0.81
#特征根: [0.9 - 3.650017528e-9*I, 0.9 + 3.650025035e-9*I]
示例7. 结构力学振动模态
beam_vibration = "lambda**4 - 12.36*lambda**2 + 12.36"
print(f"特征方程: {beam_vibration.replace('lambda', 'x')}")
roots = roots(beam_vibration.replace('lambda', 'x'))
print(f"特征根: {roots}")
#特征方程: x**4 - 12.36*x**2 + 12.36
#特征根: [-3.355973541, -1.047588464, 1.047588464, 3.355973541]
示例8. 化学反应动力学
reaction_kinetics = "k**2 + 0.5*k - 2"
print(f"特征方程: {reaction_kinetics.replace('k', 'x')}")
roots = roots(reaction_kinetics.replace('k', 'x'))
print(f"特征根: {roots}")
#特征方程: x**2 + 0.5*x - 2
#特征根: [-1.686140662, 1.186140662]
示例9. 特殊数学多项式
chebyshev = "8*x**3 - 4*x"
print(f"多项式: {chebyshev}")
roots = roots(chebyshev)
print(f"根: {roots}")
#多项式: 8*x**3 - 4*x
#根: [-0.7071067812, 0, 0.7071067812]
示例10. 复杂高阶多项式
high_order = "[1, 0, 0, 0, 0, 0, -1]" # x⁶ - 1 = 0
print(f"多项式系数: {high_order}")
roots = roots(high_order)
print(f"根: {roots}")
#多项式系数: [1, 0, 0, 0, 0, 0, -1]
#根: [-1, 1, -0.5 - 0.8660254038*I, -0.5 + 0.8660254038*I, 0.5 - 0.8660254038*I, 0.5 + 0.8660254038*I]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def roots_of_poly(input_str):
"""
对标 MATLAB 的 roots 函数，计算多项式的根
参数:
input_str (str): 多项式的字符串表示，可以是：
- 系数列表（如 '[1, -3, 2]'，对应 x²-3x+2）
- 多项式表达式（如 'x**2 - 3*x + 2'）
返回:
list[complex]: 包含所有根的列表（复数形式），数值解经过四舍五入处理
str: 错误信息字符串（当输入无效时）
实现特性:
- 支持系数列表和多项式表达式两种输入形式
- 自动处理数值计算误差（四舍五入到小数点后10位）
- 返回包含复数根的完整解集
- 错误处理包含详细的异常信息
"""
try:
# 符号变量定义
x = sp.symbols('x')
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
# 多项式对象生成逻辑
if isinstance(expr, list):
# 系数列表输入模式（MATLAB roots函数风格）
# 示例：输入 "[1, 2, 3]" 对应 x² + 2x + 3
poly = sp.Poly(expr, x)
else:
# 多项式表达式输入模式
# 示例：输入 "x**3 - 1" 对应 x³ - 1
poly = sp.Poly(expr, x)
# 计算数值根（包含复数根）
numerical_roots = poly.nroots()
# 处理数值精度问题
processed_roots = []
for root in numerical_roots:
# 改进后的SymPy原生处理
exact_root = sp.N(root, 10) # SymPy精度控制
r = sp.re(exact_root).evalf(10) # 精确实部
i = sp.im(exact_root).evalf(10) # 精确虚部
if i == 0:
processed_root = r # 纯实数输出
else:
processed_root = r + i * sp.I # SymPy复数
processed_roots.append(processed_root)
return processed_roots
except Exception as e:
# 错误处理（包含常见错误类型）
return f"错误: {str(e)}"
if __name__ == "__main__":
# 示例测试用例
test_cases = [
("[1, -3, 2]", "二次多项式（实数根）"),
# [1.000000000, 2.000000000]
("x**3 - 1", "三次单位根"),
# [1.000000000, -0.5 - 0.8660254038*I, -0.5 + 0.8660254038*I]
("x**2 + 4", "纯虚数根"),
# [-2.0*I, 2.0*I]
("[1, 0, 0, -1]", "四次多项式（x⁴-1）"),
# [1.000000000, -0.5 - 0.8660254038*I, -0.5 + 0.8660254038*I]
]
for input_str, description in test_cases:
print(f"测试用例: {description}")
print(f"输入: {input_str}")
result = roots_of_poly(input_str)
print("结果:", result)
print("-" * 60)
将数组旋转90度
B = rot90(A) 将数组A逆时针旋转90度, rot90在由第一个和第二个维度构成的平面中旋转.
B = rot90(A,k) 将数组A按逆时针方向旋转k*90度,其中k是一个整数.
A — 输入数组,向量,矩阵或多维数组
k是旋转常量,整数
示例1. 图像处理 - 图像旋转
模拟一个3x3像素块的图像矩阵
image_matrix = [[255, 128, 64], [192, 96, 32], [128, 64, 0]]
原始图像像素矩阵:
original = rot90(image_matrix)
print(original)
#[[128, 192, 255],
[64, 96, 128],
[0, 32, 64]]
逆时针旋转90度:
rotated_90 = rot90(image_matrix, 1)
print(rotated_90)
#[[128, 192, 255],
[64, 96, 128],
[0, 32, 64]]
旋转180度 (上下颠倒):
rotated_180 = rot90(image_matrix, 2)
print(rotated_180)
#[[0, 64, 128],
[32, 96, 192],
[64, 128, 255]]
顺时针旋转90度:
rotated_270 = rot90(image_matrix, -1)
print(rotated_270)
#[[64, 32, 0],
[128, 96, 64],
[255, 192, 128]]
示例2. 游戏开发 - 俄罗斯方块旋转
L形方块
tetris_L = [[1, 0], [1, 0], [1, 1]]
L形俄罗斯方块原始状态:
original_L = rot90(tetris_L)
print(original_L)
#[[1, 1, 1],
[1, 0, 0]]
旋转90度:
rotated_L = rot90(tetris_L, 1)
print(rotated_L)
#[[1, 1, 1],
[1, 0, 0]]
T形方块
tetris_T = [[0, 1, 0], [1, 1, 1]]
T形俄罗斯方块原始状态:
original_T = rot90(tetris_T)
print(original_T)
#[[1, 0],
[1, 1],
[1, 0]]
旋转180度:
rotated_T = rot90(tetris_T, 2)
print(rotated_T)
#[[1, 1, 1],
[0, 1, 0]]
示例3. 数学计算 - 矩阵变换
创建希尔伯特矩阵（数值分析中常用）
hilbert_matrix = [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]]
3x3希尔伯特矩阵:
original_hilbert = rot90(hilbert_matrix)
print(original_hilbert)
#[[1/3, 1/2, 1],
[1/4, 1/3, 1/2],
[1/5, 1/4, 1/3]]
旋转90度后的希尔伯特矩阵:
rotated_hilbert = rot90(hilbert_matrix, 1)
print(rotated_hilbert)
#[[1/3, 1/2, 1],
[1/4, 1/3, 1/2],
[1/5, 1/4, 1/3]]
特殊矩阵：帕斯卡矩阵
pascal_matrix = [[1, 1, 1, 1], [1, 2, 3, 4], [1, 3, 6, 10], [1, 4, 10, 20]]
帕斯卡矩阵:
original_pascal = rot90(pascal_matrix)
print(original_pascal)
#[[1, 1, 1, 1],
[4, 3, 2, 1],
[10, 6, 3, 1],
[20, 10, 4, 1]]
旋转270度后的帕斯卡矩阵:
rotated_pascal = rot90(pascal_matrix, 3)
print(rotated_pascal)
#[[1, 4, 10, 20],
[1, 3, 6, 10],
[1, 2, 3, 4],
[1, 1, 1, 1]]
示例4. 数据分析 - 数据透视表旋转
季度销售数据：行代表产品，列代表季度
sales_data = [[120, 150, 130, 180], [80, 90, 95, 110], [200, 220, 210, 240]]
原始销售数据 (产品×季度):
original_sales = rot90(sales_data)
print(original_sales)
#[[200, 80, 120],
[220, 90, 150],
[210, 95, 130],
[240, 110, 180]]
旋转90度 (季度×产品):
rotated_sales = rot90(sales_data, 1)
print(rotated_sales)
#[[200, 80, 120],
[220, 90, 150],
[210, 95, 130],
[240, 110, 180]]
旋转180度 (数据完全转置):
full_rotate_sales = rot90(sales_data, 2)
print(full_rotate_sales)
#[[240, 210, 220, 200],
[110, 95, 90, 80],
[180, 130, 150, 120]]
示例5. 计算机图形学 - 坐标系变换
2D图形顶点坐标
triangle_vertices = [[0, 0], [1, 0], [0.5, 1]]
三角形原始顶点坐标:
original_triangle = rot90(triangle_vertices)
print(original_triangle)
#[[0.5, 1, 0],
[1, 0, 0]]
逆时针旋转90度后的顶点:
rotated_triangle = rot90(triangle_vertices, 1)
print(rotated_triangle)
#[[0.5, 1, 0],
[1, 0, 0]]
顺时针旋转90度后的顶点:
clockwise_triangle = rot90(triangle_vertices, -1)
print(clockwise_triangle)
#[[0, 0, 1],
[0, 1, 0.5]]
示例6. 密码学 - 矩阵编码
简单的替换密码矩阵
cipher_matrix = [[7, 4, 2], [1, 5, 8], [3, 9, 6]]
原始密码矩阵:
original_cipher = rot90(cipher_matrix)
print(original_cipher)
#[[3, 1, 7],
[9, 5, 4],
[6, 8, 2]]
旋转90度后的密码矩阵:
rotated_cipher = rot90(cipher_matrix, 1)
print(rotated_cipher)
#[[3, 1, 7],
[9, 5, 4],
[6, 8, 2]]
示例7. 物理模拟 - 应力张量旋转
2D应力张量
stress_tensor = [[100, -50], [-50, 80]]
原始应力张量 (MPa):
original_stress = rot90(stress_tensor)
print(original_stress)
#[[-50, 100],
[80, -50]]
坐标系旋转45度后的应力张量:
rotated_stress = rot90(stress_tensor, 1)
print(rotated_stress)
#[[-50, 100],
[80, -50]]
坐标系旋转90度后的应力张量:
stress_90 = rot90(stress_tensor, 2)
print(stress_90)
#[[80, -50],
[-50, 100]]
示例8. 机器学习 - 特征矩阵变换
特征矩阵：行代表样本，列代表特征
feature_matrix = [[1.2, 3.4, 5.6], [2.1, 4.3, 6.5], [0.9, 2.8, 4.7]]
原始特征矩阵 (样本×特征):
original_features = rot90(feature_matrix)
print(original_features)
#[[0.9, 2.1, 1.2],
[2.8, 4.3, 3.4],
[4.7, 6.5, 5.6]]
旋转90度后的特征矩阵 (特征×样本):
rotated_features = rot90(feature_matrix, 1)
print(rotated_features)
#[[0.9, 2.1, 1.2],
[2.8, 4.3, 3.4],
[4.7, 6.5, 5.6]]
#这种旋转在某些数据预处理和特征工程中有用
示例9. 电子表格处理 - 表格转置
学生成绩表：行代表学生，列代表科目
grade_sheet = [[85, 92, 78], [76, 88, 95], [92, 85, 90], [68, 75, 82]]
原始成绩表 (学生×科目):
original_grades = rot90(grade_sheet)
print(original_grades)
#[[68, 92, 76, 85],
[75, 85, 88, 92],
[82, 90, 95, 78]]
旋转180度 (科目×学生):
rotated_grades = rot90(grade_sheet, 2)
print(rotated_grades)
#[[82, 75, 68],
[90, 85, 92],
[95, 88, 76],
[78, 92, 85]]
#这种变换在数据可视化和报告生成中很有用
示例10. 棋盘游戏 - 游戏状态旋转
井字棋游戏状态 (0:空, 1:X, 2:O)
tic_tac_toe = [[1, 0, 2], [0, 1, 0], [2, 0, 1]]
井字棋当前状态:
original_board = rot90(tic_tac_toe)
print(original_board)
#[[2, 0, 1],
[0, 1, 0],
[1, 0, 2]]
旋转90度后的状态:
rotated_board = rot90(tic_tac_toe, 1)
print(rotated_board)
#[[2, 0, 1],
[0, 1, 0],
[1, 0, 2]]
#这种旋转在游戏AI中用于状态对称性检测
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def rot90_matrix(input_str):
"""
对标 MATLAB 的 rot90 函数，实现矩阵旋转功能
参数:
input_str (str): 输入字符串表达式，支持格式：
- 矩阵 (如 "[[1,2],[3,4]]")
- 元组格式 (矩阵, 旋转次数) (如 "([[1,2],[3,4]], 2)")
返回:
sp.Matrix: 旋转后的矩阵 (成功时)
str: 错误信息 (失败时)
功能特性:
- 支持逆时针旋转 (k>0)
- 支持顺时针旋转 (自动转换负数为等效正数)
- 自动处理整数溢出 (k 取模 4)
- 输入格式自动验证
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
# 处理元组输入 (矩阵, 旋转次数)
if isinstance(expr, tuple) and len(expr) == 2:
matrix = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
k_expr = expr[1]
# 验证旋转次数是否为整数
if matrix is not None and k_expr.is_Integer:
k = int(k_expr) % 4 # 规范化旋转次数
result = matrix.rot90(k)
else:
error = True
# 处理单个矩阵输入 (默认旋转1次)
else:
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is not None:
result = matrix.rot90(1)
else:
error = True
return result if not error else f"输入格式错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
if __name__ == "__main__":
# 示例测试用例
test_cases = [
("[[1,2],[3,4]]", "默认旋转1次 (逆时针90度)"),
# Matrix([[3, 1],
# [4, 2]])
("([[1,2],[3,4]], 2)", "旋转2次 (180度)"),
# Matrix([[4, 3],
# [2, 1]])
("([[1,2,3],[4,5,6]], -1)", "负数旋转 (顺时针90度)"),
# Matrix([[3, 6],
# [2, 5],
# [1, 4]])
("([1,2,3], 1)", "列向量旋转"),
# Matrix([[3, 2, 1]])
]
for input_str, desc in test_cases:
print(f"测试用例: {desc}")
print(f"输入: {input_str}")
result = rot90_matrix(input_str)
print("输出:" if isinstance(result, sp.Matrix) else "结果:", result)
print("=" * 60)
舍入至最近的小数或整数
Y = round(X)将X的每个元素四舍五入为最近的整数.在舍入机会均等的情况下,即有元素的十进制小数部分为0.5(在舍入误差内)时,round函数会偏离零四舍五入到最接近的具有更大幅值的整数.
Y = round(X,N) 四舍五入到N位数:
N > 0：舍入到小数点右侧的第N位数.
N = 0：四舍五入到最接近的整数.
N < 0：舍入到小数点左侧的第N位数.
X — 输入数组,标量,向量,矩阵
N — 位数,整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def round_decimal_place(input_str):
"""
对标 MATLAB 的 round 函数，实现数值舍入功能
参数:
input_str (str): 输入字符串表达式，支持格式：
- 单个数值/表达式 (如 "3.1415" 或 "sqrt(2)")
- 矩阵 (如 "[[1.2, 3.4], [5.6, 7.8]]")
- 元组格式 (数值/矩阵, 小数位数) (如 "(3.1415, 2)")
返回:
float/sp.Matrix: 舍入后的数值或矩阵 (成功时)
str: 错误信息 (失败时)
功能特性:
- 支持矩阵元素级舍入
- 自动处理符号表达式数值化
- 支持正/负数舍入位数
- 自动验证输入格式有效性
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
def evaluation_round(x, n=None):
"""核心舍入函数，处理数值和符号表达式"""
# 处理纯数值
return x.round(n) if n else x.round()
# 处理元组输入 (数值/矩阵, 小数位数)
if isinstance(expr, tuple) and len(expr) == 2:
x_part, n_part = expr
# 验证小数位数是否为整数
if not n_part.is_Integer:
error = True
else:
n = int(n_part)
# 处理矩阵输入
if isinstance(x_part, list):
matrix = sp.Matrix(x_part)
result = sp.Matrix(matrix.shape[0], matrix.shape[1],
lambda i, j: evaluation_round(matrix[i, j], n))
# 处理数值/符号输入
elif x_part.is_number or x_part.free_symbols:
result = evaluation_round(x_part, n)
else:
error = True
# 处理单个元素输入
else:
# 处理矩阵输入
if isinstance(expr, list):
matrix = sp.Matrix(expr)
result = sp.Matrix(matrix.shape[0], matrix.shape[1],
lambda i, j: evaluation_round(matrix[i, j]))
# 处理数值/符号输入
elif expr.is_number:
result = evaluation_round(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
if __name__ == "__main__":
# 示例测试用例
test_cases = [
("3.1415926", "基本数值舍入"),
# 3
("([[1.234, 5.678], [9.012, 3.456]], 1)", "矩阵保留1位小数"),
# Matrix([[1.2, 5.7],
# [9.00000000000000, 3.5]])
("([2.345, 6.789], 2)", "列向量保留2位小数"),
# Matrix([[2.34],
# [6.79]])
("(-3.5, 0)", "负数舍入到整数"),
# -4
]
for input_str, desc in test_cases:
print(f"测试用例: {desc}")
print(f"输入: {input_str}")
result = round_decimal_place(input_str)
if isinstance(result, sp.Matrix):
print("输出矩阵:")
print(result)
else:
print("结果:", result)
print("=" * 60)
简化的行阶梯形矩阵(高斯消去法)
R = rref(A) 使用 Gauss-Jordan 消去法和部分主元消去法返回简化行阶梯形的A
A — 输入矩阵, 矩阵
示例1. 线性方程组求解
方程组:
2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3
增广矩阵
augmented_matrix = [[2, 1, -1, 8], [-3, -1, 2, -11], [-2, 1, 2, -3]]
result = rref(augmented_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 0, 2],
[0, 1, 0, 3],
[0, 0, 1, -1]]
解:
print(f"x = {result[0, 3]}")
print(f"y = {result[1, 3]}")
print(f"z = {result[2, 3]}")
#x = 2
#y = 3
#z = -1
示例2. 矩阵秩的计算
matrix_a = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
result = rref(matrix_a)
简化行阶梯形:
print(result)
#[[1, 0, -1],
[0, 1, 2],
[0, 0, 0]]
计算非零行的数量（矩阵的秩）
rank = sum(1 for i in range(result.shape[0]) if any(result[i, j] != 0 for j in range(result.shape[1])))
print(f"矩阵的秩: {rank}")
#矩阵的秩: 2
另一个例子 - 满秩矩阵
matrix_b = [[1, 0, 2], [-1, 1, 1], [0, 1, 3]]
result_b = rref(matrix_b)
简化行阶梯形:
print(result_b)
#[[1, 0, 2],
[0, 1, 3],
[0, 0, 0]]
rank_b = sum(1 for i in range(result_b.shape[0]) if any(result_b[i, j] != 0 for j in range(result_b.shape[1])))
print(f"矩阵的秩: {rank_b}")
#矩阵的秩: 2
示例3. 向量组的线性相关性判断
vectors = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
result = rref(vectors)
简化行阶梯形:
print(result)
#[[1, 0, -1],
[0, 1, 2],
[0, 0, 0]]
判断线性相关性
if result.shape[0] == result.shape[1] and result == sp.eye(result.shape[0]):
print("结论: 向量组线性无关")
else:
print("结论: 向量组线性相关")
print("存在非平凡线性组合使得向量和为0")
#结论: 向量组线性相关
#存在非平凡线性组合使得向量和为0
示例4. 齐次线性方程组的解空间
齐次方程组:
x + 2y + 3z = 0
2x + 4y + 6z = 0
3x + 6y + 9z = 0
homogeneous_matrix = [[1, 2, 3], [2, 4, 6], [3, 6, 9]]
result = rref(homogeneous_matrix)
简化行阶梯形:
print(result)
#[[1, 2, 3],
[0, 0, 0],
[0, 0, 0]]
分析解空间
pivot_columns = []
for i in range(min(result.shape)):
for j in range(result.shape[1]):
if result[i, j] != 0:
pivot_columns.append(j)
break
free_vars = set(range(result.shape[1])) - set(pivot_columns)
print(f"主元列: {pivot_columns}")
print(f"自由变量: {free_vars}")
print(f"解空间的维数: {len(free_vars)}")
#主元列: [0]
#自由变量: {1, 2}
#解空间的维数: 2
示例5. 电路分析 - 基尔霍夫定律
电路节点电流方程:
I1 - I2 - I3 = 0
I2 - I4 - I5 = 0
10I1 + 5I2 = 20
5I2 + 10I3 + 5I4 = 0
5I4 + 10I5 = 10
circuit_matrix = [[1, -1, -1, 0, 0, 0], [0, 1, 0, -1, -1, 0], [10, 5, 0, 0, 0, 20], [0, 5, 10, 5, 0, 0], [0, 0, 0, 5, 10, 10]]
result = rref(circuit_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 0, 0, -1/2, 0],
[0, 1, 0, 0, 1, 0],
[0, 0, 1, 0, -3/2, 0],
[0, 0, 0, 1, 2, 0],
[0, 0, 0, 0, 0, 1]]
各支路电流:
currents = [I1, I2, I3, I4, I5]
for i in range(min(result.shape[0], len(currents))):
print(f"{currents[i]} = {result[i, -1]:.4f} A")
#I1 = 0.0000 A
#I2 = 0.0000 A
#I3 = 0.0000 A
#I4 = 0.0000 A
#I5 = 1.0000 A
示例6. 经济学 - 投入产出分析
三部门经济模型:
农业: x1 = 0.2x1 + 0.3x2 + 0.2x3 + 50
工业: x2 = 0.4x1 + 0.1x2 + 0.2x3 + 100
服务业: x3 = 0.3x1 + 0.2x2 + 0.1x3 + 80
重写为齐次形式: (I-A)x = d
input_output_matrix = "[[0.8, -0.3, -0.2, 50], [-0.4, 0.9, -0.2, 100], [-0.3, -0.2, 0.9, 80]]"
result = rref(input_output_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 0, 211.190476190476],
[0, 1, 0, 252.857142857143],
[0, 0, 1, 215.476190476190]]
各部门最优产出:
sectors = [农业, 工业, 服务业]
for i in range(min(result.shape[0], len(sectors))):
print(f"{sectors[i]}: {result[i, -1]:.2f}")
#农业: 211.19
#工业: 252.86
#服务业: 215.48
示例7. 计算机图形学 - 齐次坐标变换
已知4个点的投影前后坐标，求解投影矩阵
projection_matrix = [[100, 200, 1, 0, 0, 0, -150, -300],
[0, 0, 0, 100, 200, 1, -50, -100],
[300, 100, 1, 0, 0, 0, -450, -150],
[0, 0, 0, 300, 100, 1, -150, -50],
[200, 300, 1, 0, 0, 0, -200, -300],
[0, 0, 0, 200, 300, 1, -200, -300]]
result = rref(projection_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 0, 0, 0, 0, -7/6, 1/2],
[0, 1, 0, 0, 0, 0, 2/3, -1/2],
[0, 0, 1, 0, 0, 0, -500/3, -250],
[0, 0, 0, 1, 0, 0, -5/6, -1/2],
[0, 0, 0, 0, 1, 0, -2/3, -3/2],
[0, 0, 0, 0, 0, 1, 500/3, 250]]
投影变换矩阵参数:
params = [a, b, c, d, e, f, g, h]
for i in range(min(result.shape[0], len(params))):
print(f"{params[i]} = {result[i, -1]:.4f}")
#a = 0.5
#b = -0.5
#c = -250
#d = -0.5
#e = -1.5
#f = 250
示例8. 机器学习 - 线性回归参数估计
多元线性回归: y = β₀ + β₁x₁ + β₂x₂
正规方程: (XᵀX)β = Xᵀy
regression_matrix = [[5, 10, 20, 30], [10, 30, 60, 100], [20, 60, 140, 210], [30, 100, 210, 350]]
result = rref(regression_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]]
回归系数:
coefficients = [β₀, β₁, β₂]
for i in range(min(result.shape[0], len(coefficients))):
print(f"{coefficients[i]} = {result[i, -1]:.4f}")
#β₀ = 0.0000
#β₁ = 0.0000
#β₂ = 0.0000
示例9. 密码学 - 线性码的解码
假设我们有一个(7,4)汉明码的校验矩阵
hamming_matrix = [[1, 0, 1, 0, 1, 0, 1, 0], [0, 1, 1, 0, 0, 1, 1, 0], [0, 0, 0, 1, 1, 1, 1, 0]]
result = rref(hamming_matrix)
简化行阶梯形:
print(result)
#[[1, 0, 1, 0, 1, 0, 1, 0],
[0, 1, 1, 0, 0, 1, 1, 0],
[0, 0, 0, 1, 1, 1, 1, 0]]
校验结果分析:
if all(result[i, -1] == 0 for i in range(result.shape[0])):
print("接收到的码字正确，无错误")
else:
print("检测到错误，需要纠错")
error_position = None
for i in range(result.shape[0]):
if result[i, -1] != 0:
# 查找错误模式
for j in range(result.shape[1] - 1):
if result[i, j] != 0:
error_position = j
break
if error_position is not None:
break
print(f"错误可能出现在第 {error_position + 1} 位")
#校验结果分析:
#接收到的码字正确，无错误
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def row_reduce_gauss(input_str):
"""
对标 MATLAB 的 rref 函数，计算简化行阶梯形矩阵
参数:
input_str (str): 输入字符串表达式，支持格式：
- 数值矩阵 (如 "[[1,2],[3,4]]")
- 符号矩阵 (如 "[[x, 1], [2, y]]")
- 增广矩阵 (如 "[[1,2,5],[3,4,6]]")
返回:
sp.Matrix: 简化行阶梯形矩阵 (成功时)
str: 错误信息 (失败时)
功能特性:
- 支持符号运算
- 自动处理分数形式
- 保留行交换信息
- 完整错误处理链
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str)
# 转换为矩阵对象
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is None:
return f"输入错误: 无法转换为矩阵 - {input_str}"
# 执行高斯消去法
reduced_matrix, pivots = matrix.rref()
return reduced_matrix
except Exception as e:
return f"错误: {str(e)}"
if __name__ == "__main__":
# 示例测试用例
test_cases = [
("[[1,2],[3,4]]", "数值矩阵"),
# Matrix([[1, 0],
# [0, 1]])
("[[1,2,5],[3,4,6]]", "增广矩阵"),
# Matrix([[1, 0, -4],
# [0, 1, 9/2]])
("[[x,1],[2,y]]", "符号矩阵"),
# Matrix([[1, 0],
# [0, 1]])
("[1,2,3]", "行向量"),
# Matrix([[1],
# [0],
# [0]])
]
for input_str, desc in test_cases:
print(f"测试用例: {desc}")
print(f"输入: {input_str}")
result = row_reduce_gauss(input_str)
if isinstance(result, sp.Matrix):
print("简化行阶梯形:")
print(result)
else:
print("结果:", result)
print("=" * 60)
将实数舒尔形式转换为复数舒尔形式
[Unew,Tnew] = rsf2csf(U,T) 将实矩阵 X 的 [U,T] = schur(X) 的输出从实数舒尔形式变换为复数舒尔形式。
此操作变换 X 的特征值在 T 中的表示方式，并变换 U 以使 X = Unew*Tnew*Unew' 且 Unew'*Unew = eye(size(X))。.
U — 酉矩阵, 矩阵
T — 舒尔形式, 矩阵
Unew — 变换后的酉矩阵, 矩阵, 矩阵 Unew 满足 Unew'*Unew = eye(size(X))。
Tnew — 变换后的舒尔形式, 矩阵, Tnew 是上三角矩阵，X 的特征值在对角线上，且它满足 X = Unew*Tnew*Unew'。
示例1：机械振动系统的质量-弹簧-阻尼系统
自由度振动系统的实数舒尔分解结果
T_vibration = [
[-0.25, 8.66, 0, 0],
[-8.66, -0.25, 0, 0],
[0, 0, -0.5, 5.0],
[0, 0, -5.0, -0.5]
]
U_vibration = [
[0.7071, 0, 0.7071, 0],
[0, 0.7071, 0, 0.7071],
[0.7071, 0, -0.7071, 0],
[0, 0.7071, 0, -0.7071]
]
U_csf1, T_csf1 = rsf2csf(U_vibration, T_vibration)
复数舒尔形式 T_csf (振动系统):
print(T_csf1)
#[[0.499995204977008*I, 0.499995204977008, 0.499995204977008*I, 0.499995204977008],
[-0.499995204977008, -0.499995204977008*I, -0.499995204977008, -0.499995204977008*I],
[0.499995204977008*I, 0.499995204977008, -0.499995204977008*I, -0.499995204977008],
[-0.499995204977008, -0.499995204977008*I, 0.499995204977008, 0.499995204977008*I]]
酉矩阵 U_csf (振动系统):
print(U_csf1)
#[[-0.25 + 8.66*I, 2.66453525910038e-15, 0, 0],
[0, -0.25 - 8.66*I, 0, 0],
[0, 0, -0.5 + 5.0*I, -8.88178419700125e-16],
[0, 0, 0, -0.5 - 5.0*I]]
示例2：RLC电路系统的状态空间表示
RLC电路的实数舒尔形式
T_circuit = [
[-100, 314.16, 0, 0],
[-314.16, -100, 0, 0],
[0, 0, -50, 628.32],
[0, 0, -628.32, -50]
]
U_circuit = [
[0.8660, -0.5, 0, 0],
[0.5, 0.8660, 0, 0],
[0, 0, 0.8660, -0.5],
[0, 0, 0.5, 0.8660]
]
U_csf2, T_csf2 = rsf2csf(U_circuit, T_circuit)
复数舒尔形式 T_csf (电路系统):
print(T_csf2)
#[[0.353553390593274 + 0.61235447250755*I, 0.61235447250755 + 0.353553390593274*I, 0, 0],
[-0.61235447250755 + 0.353553390593274*I, 0.353553390593274 - 0.61235447250755*I, 0, 0],
[0, 0, 0.353553390593274 + 0.61235447250755*I, 0.61235447250755 + 0.353553390593274*I],
[0, 0, -0.61235447250755 + 0.353553390593274*I, 0.353553390593274 - 0.61235447250755*I]]
酉矩阵 U_csf (电路系统):
print(U_csf2)
#[[-100.0 + 314.16*I, 1.70530256582424e-13, 0, 0],
[0, -100.0 - 314.16*I, 0, 0],
[0, 0, -50.0 + 628.32*I, 2.27373675443232e-13 - 5.6843418860808e-14*I],
[0, 0, 0, -50.0 - 628.32*I]]
示例3：控制系统极点配置
控制系统的实数舒尔形式
T_control = [
[-2.0, 3.0, 0, 0, 0],
[-3.0, -2.0, 0, 0, 0],
[0, 0, -1.0, 1.5, 0],
[0, 0, -1.5, -1.0, 0],
[0, 0, 0, 0, -0.5]
]
U_control = [
[0.8944, 0.4472, 0, 0, 0],
[-0.4472, 0.8944, 0, 0, 0],
[0, 0, 0.8321, 0.5547, 0],
[0, 0, -0.5547, 0.8321, 0],
[0, 0, 0, 0, 1.0]
]
U_csf3, T_csf3 = rsf2csf(U_control, T_control)
复数舒尔形式 T_csf (控制系统):
print(T_csf3)
#[[-0.316218152546624 + 0.632436305093248*I, 0.632436305093248 - 0.316218152546624*I, 0, 0, 0],
[-0.632436305093248 - 0.316218152546624*I, -0.316218152546624 - 0.632436305093248*I, 0, 0, 0],
[0, 0, -0.392232131524178 + 0.588383552625326*I, 0.588383552625326 - 0.392232131524178*I, 0],
[0, 0, -0.588383552625326 - 0.392232131524178*I, -0.392232131524178 - 0.588383552625326*I, 0],
[0, 0, 0, 0, 1.0]]
示例4：建筑结构动力学
3层建筑结构的实数舒尔形式
T_structure = [
[-0.1, 12.57, 0, 0, 0, 0],
[-12.57, -0.1, 0, 0, 0, 0],
[0, 0, -0.2, 25.13, 0, 0],
[0, 0, -25.13, -0.2, 0, 0],
[0, 0, 0, 0, -0.3, 37.70],
[0, 0, 0, 0, -37.70, -0.3]
]
U_structure = [
[0.5, 0.8660, 0, 0, 0, 0],
[-0.8660, 0.5, 0, 0, 0, 0],
[0, 0, 0.5, 0.8660, 0, 0],
[0, 0, -0.8660, 0.5, 0, 0],
[0, 0, 0, 0, 0.5, 0.8660],
[0, 0, 0, 0, -0.8660, 0.5]
]
U_csf4, T_csf4 = rsf2csf(U_structure, T_structure)
复数舒尔形式 T_csf (结构动力学):
print(T_csf4)
#[[-0.61235447250755 + 0.353553390593274*I, 0.353553390593274 - 0.61235447250755*I, 0, 0, 0, 0],
[-0.353553390593274 - 0.61235447250755*I, -0.61235447250755 - 0.353553390593274*I, 0, 0, 0, 0],
[0, 0, -0.61235447250755 + 0.353553390593274*I, 0.353553390593274 - 0.61235447250755*I, 0, 0],
[0, 0, -0.353553390593274 - 0.61235447250755*I, -0.61235447250755 - 0.353553390593274*I, 0, 0],
[0, 0, 0, 0, -0.61235447250755 + 0.353553390593274*I, 0.353553390593274 - 0.61235447250755*I],
[0, 0, 0, 0, -0.353553390593274 - 0.61235447250755*I, -0.61235447250755 - 0.353553390593274*I]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy import linalg
import numpy as np
def real_schur_complex_schur(input_str):
"""
将实数舒尔形式转换为复数舒尔形式（参数顺序调整为MATLAB风格）
参数:
U (sympy.Matrix): 正交变换矩阵
T (sympy.Matrix): 实数舒尔形式的上准三角矩阵
返回:
tuple: (U_csf, T_csf) 与MATLAB一致
异常:
ValueError: 如果输入矩阵无效
"""
try:
# 检查输入矩阵
expr = sp.sympify(input_str)
T, U = None, None
# 检查表达式是否为元组且长度为 2
if isinstance(expr, tuple) and len(expr) == 2:
# 如果是，则将元组中的元素分别赋值给 U 和 T
U, T = expr
else:
# 若不满足条件，返回输入参数错误的提示信息
return f"输入参数不对。"
# 将 sympy 矩阵 T 转换为 numpy 数组，数据类型为 float
if isinstance(T, list):
T_np = np.array(T, dtype=float)
# 将 sympy 矩阵 U 转换为 numpy 数组，数据类型为 float
if isinstance(U, list):
U_np = np.array(U, dtype=float)
# 调用 scipy.linalg 的 rsf2csf 函数，将实数舒尔形式转换为复数舒尔形式
T2, U2 = linalg.rsf2csf(T_np, U_np)
# 将转换后的 numpy 数组 T2 和 U2 转换为 sympy 矩阵并返回
return sp.Matrix(T2), sp.Matrix(U2)
except Exception as e:
raise ValueError(f"转换失败: {str(e)}")
# 修改后的示例代码
if __name__ == "__main__":
# 示例矩阵（调整参数顺序为U,T）
T_real = [
[4.8121, 1.1972, -2.2273, -1.0067],
[0, 1.9202, -3.0485, -1.8381],
[0, 0.7129, 1.9202, 0.2566],
[0, 0, 0, 1.3474]
]
U_real = [
[-0.4916, -0.49, -0.6331, -0.3428],
[-0.4980, 0.2403, -0.2325, 0.8001],
[-0.6751, 0.4288, 0.4230, -0.4260],
[-0.2337, -0.7200, 0.6052, 0.2466]
]
# 直接传入矩阵对象（不再使用字符串）
U_csf, T_csf = real_schur_complex_schur(f"{U_real},{T_real}")
print("复数舒尔形式 T_csf:")
print(T_csf)
# Matrix([[-0.491600000000000, -0.275620718655525 - 0.441127791460706*I, 0.213321990429959 + 0.56995511178321*I, -0.342800000000000],
# [-0.498000000000000, -0.101219107704011 + 0.216332669975526*I, -0.104614845510856 + 0.209310635744111*I, 0.800100000000000],
# [-0.675100000000000, 0.184153473371169 + 0.386031830568063*I, -0.186678509176258 - 0.380810317934446*I, -0.426000000000000],
# [-0.233700000000000, 0.263474425731044 - 0.648187775207568*I, 0.313452720631777 - 0.544837835493917*I, 0.246600000000000]])
print("\n酉矩阵 U_csf:")
print(U_csf)
# Matrix([[4.81210000000000, -0.969657284254384 + 1.07779222844236*I, -0.521202218250505 + 2.00515087738864*I, -1.00670000000000],
# [0, 1.9202 + 1.47420339505782*I, 2.33560000000000, 0.111711066825158 + 1.65476937445699*I],
# [0, 0, 1.9202 - 1.47420339505782*I, 0.800218674712874 + 0.231006920997586*I], [0, 0, 0, 1.34740000000000]])
黎曼和
rsums(f) 以交互方式通过 x 从 0 到 1 的中间黎曼和来近似 f(x) 的积分.
f — 被积函数,符号表达式,符号函数,符号数
a — 下限数字,符号数字.
b — 上限数字,符号数字.
示例1. 物理学：计算位移（速度积分）
物体以速度 v(t) = 5 m/s 运动，计算 [0, 10] 秒内的位移"
result1 = rsums(5, 0, 10)
print(f"近似位移: {result1} 米")
#近似位移: 50 米
示例2. 物理学：变加速运动的位移
加速度 a(t) = 2t m/s²，计算 [0, 5] 秒内的速度变化
result2 = rsums(2*x, 0, 5)
print(f"速度变化: {result2} m/s")
#速度变化: 25 m/s
示例3. 经济学：边际收益函数的总收益
边际收益 MR(q) = 100 - 2q，计算产量从 0 到 40 的总收益
result3 = rsums(100 - 2*x, 0, 40)
print(f"总收益: {result3}")
#总收益: 2400
示例4. 工程学：梁的弯矩计算
分布载荷 w(x) = x(4-x) kN/m，计算 [0, 4] 米上的总载荷
result4 = rsums(x*(4-x), 0, 4)
print(f"总载荷: {result4} kN")
#总载荷: 10.6672 kN
示例5. 概率论：概率密度函数积分
f(x) = x 在 [0, 1] 上的积分（均匀分布的期望）
result5 = rsums(x)
print(f"期望值: {result5}")
#期望值: 0.5
示例6. 生物学：种群增长模型
增长率 r(t) = 100/(1+t) 个体/年，计算 [0, 10] 年的总增长
result6 = rsums(100/(1+x), 0, 10)
print(f"总增长: {result6} 个体")
#总增长: 239.748277497110 个体
示例7. 化学：反应速率积分
反应速率 v(t) = e^(-0.1t) mol/s，计算 [0, 20] 秒内的总产物
result7 = rsums(exp(-0.1*x), 0, 20)
print(f"总产物: {result7} mol")
#总产物: 8.64650305852902 mol
示例8. 金融学：连续复利的终值
利率函数 r(t) = 0.05，计算 10 年后的增长因子
result8 = rsums(exp(0.05*x), 0, 10)
print(f"增长因子: {result8}")
#增长因子: 12.9744118989859
示例9. 环境科学：河流污染物总量
浓度分布 c(x) = 50*exp(-0.2*x) mg/L，计算 [0, 10] km 的总量
result9 = rsums(50*exp(-0.2*x), 0, 10)
print(f"污染物总量: {result9} mg")
#污染物总量: 216.162576463225 mg
示例10. 医学：药物代谢曲线下面积(AUC)
血药浓度 C(t) = 100*exp(-0.3*t) mg/L，计算 [0, 24] 小时的AUC
result10 = rsums(100*exp(-0.3*x), 0, 24)
print(f"AUC: {result10} mg·h/L")
#AUC: 333.012536028172 mg·h/L
示例11. 电学：交流电功率计算
瞬时功率 p(t) = sin²(t) W，计算 [0, π] 秒内的总能量
result11 = rsums(sin(x)**2, 0, @pi)
print(f"总能量: {result11} J")
#总能量: 1.57079632679490 J
示例12. 热力学：变温过程的热量
热容 C(T) = 2 + 0.1T J/°C，计算从 0°C 到 100°C 的热量
result12 = rsums(2 + 0.1*x, 0, 100)
print(f"总热量: {result12} J")
#总热量: 700.000000000000 J
示例13. 流体力学：变截面管道流量
流速 v(r) = 1 - r² m/s（抛物线分布），计算半径 [0, 1] 的流量
result13 = rsums(2*pi*x*(1-x**2), 0, 1)
print(f"总流量: {result13} m³/s")
#总流量: 1.57087486661124 m³/s
示例14. 声学：正弦声波的能量
声强 I(t) = sin²(2πt) W/m²，计算一个周期 [0, 1] 内的总能量
result14 = rsums(sin(2*@pi*x)**2, 0, 1)
print(f"周期能量: {result14} J/m²")
#周期能量: 0.5 J/m²
示例15. 光学：非均匀光源的光通量
光强分布 I(x) = cos(πx/2) W/sr，计算 [0, 1] 弧度内的总光通量
result15 = rsums(cos(@pi*x/2), 0, 1)
print(f"总光通量: {result15} W")
#总光通量: 0.636626317399378 W
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def riemann_sums_evaluation(input_str):
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
def evaluation_riemann_sum(f, a=0, b=1):
"""
计算函数f在区间[a, b]上的左端点黎曼和。
参数:
f: SymPy表达式，被积函数
a: 积分下限
b: 积分上限
返回:
黎曼和的数值结果，无法计算时返回None
"""
# 确定积分变量
method = 'midpoint'
n = 100
symbols = f.free_symbols
if len(symbols) == 0:
x = sp.Symbol('x') # 处理常数函数
elif len(symbols) == 1:
x = symbols.pop()
else:
return None
dx = (b - a) / n
riemann_sum = 0
for i in range(n):
if method == 'left':
x_i = a + i * dx
elif method == 'midpoint':
x_i = a + (i + 0.5) * dx # 中点法修正
elif method == 'right':
x_i = a + (i + 1) * dx
else:
raise ValueError("Invalid method")
riemann_sum += f.subs(x, x_i) * dx
return sp.N(riemann_sum)
# 处理输入为元组（函数, a, b）的情况
if isinstance(expr, tuple) and len(expr) == 3:
f_part, a_part, b_part = expr
# 检查函数变量数和a,b是否为数值
if (len(f_part.free_symbols) <= 1 and
a_part.is_number and
b_part.is_number):
result = evaluation_riemann_sum(f_part, a_part.evalf(), b_part.evalf())
error = result is None
else:
error = True
# 处理输入为单个表达式的情况（使用默认区间[0, 1]）
elif isinstance(expr, sp.Expr):
result = evaluation_riemann_sum(expr)
error = result is None
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1：x^2在[0,1]的黎曼和
print(riemann_sums_evaluation("x**2"))
#0.333325000000000
# 示例2：sin(x)在[0, pi]的积分
print(riemann_sums_evaluation("(sin(x), 0, pi)"))
# 2.00008224907099
# 示例3：常数函数3在[0,2]的积分
print(riemann_sums_evaluation("(3, 0, 2)"))
# 5.99999999999999
# 示例4：a > b的情况
print(riemann_sums_evaluation("(x, 2, 1)"))
# -1.50000000000000