pade 帕德近似
帕德近似
pade(f)在var＝0时返回表达式f的三阶Padé逼近(不知道三阶逼近是否是pade定义,后续查证,目前代码里使用泰勒六阶展开).
f —— 近似输入,符号函数|符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.interpolate import pade
import numpy as np
def coefficients_of_polynomial(input_str):
"""
提取多项式的系数列表
参数:
input_str - 输入字符串，支持格式：
1. 纯多项式表达式 (如 "x**2 + 2*x + 1")
2. 元组格式 (表达式, 变量) (如 "(x**2 + y, x)")
3. 元组格式 (表达式, [变量1, 变量2,...]) (如 "(x**2 + x*y, [x, y])")
locals_dict - 符号解析字典，默认为None
返回:
系数列表 (单变量) 或嵌套系数列表 (多变量)
错误时返回描述性错误信息
"""
try:
expr = sp.sympify(input_str)
result = []
if isinstance(expr, tuple) and len(expr) == 2:
poly_expr, variables = expr
# 处理多变量情况
if isinstance(variables, list):
for var in variables:
if not var.free_symbols:
raise ValueError(f"无效变量: {var}")
poly = sp.poly(poly_expr, var)
result.append(poly.coeffs())
# 处理单变量情况
else:
if not variables.free_symbols:
raise ValueError(f"无效变量: {variables}")
poly = sp.poly(poly_expr, variables)
result = poly.coeffs()
else:
poly = sp.poly(expr)
result = poly.coeffs()
return result
except Exception as e:
return f"错误: {e}"
def pade_approximant(input_str):
"""
对标 MATLAB 的 pade 函数，计算符号表达式的 Pade 近似式
参数格式：
"(表达式, 变量, 展开点, 分子阶数m, 分母阶数n)"
示例："(sin(x), x, 0, 3, 2)"
返回：
SymPy表达式 或 错误信息字符串
"""
try:
expr = sp.sympify(input_str, evaluate=False)
error = False
result = None
# 计算泰勒展开（确保阶数足够）
sym_var = None
x0 = 0
m = 2
n = 2
if isinstance(expr, tuple):
if len(expr) == 2:
f_expr, sym_var = expr
elif len(expr) == 3:
f_expr, sym_var, x0 = expr
else:
error = True
elif expr.free_symbols:
f_expr = expr
sym_var = list(expr.free_symbols)[0]
else:
error = True
taylor_expr = sp.series(f_expr, sym_var, x0).removeO()
x = sp.symbols('x')
# 提取泰勒级数的系数
coeffs = coefficients_of_polynomial(str(taylor_expr))
pade_approximant = pade(np.array(coeffs, dtype=float), m=m) # 分母多项式的阶数为 n-m
# 获取分子和分母多项式
numerator, denominator = pade_approximant
# 构建分子和分母多项式的表达式
numerator_expr = sum(coef * x ** i for i, coef in enumerate(numerator.c))
denominator_expr = sum(coef * x ** i for i, coef in enumerate(denominator.c))
# 构建 Pade 近似的符号表达式
pade_expr = numerator_expr / denominator_expr
# 打印 Pade 近似的符号表达式
result = sp.nsimplify(pade_expr, tolerance=0.01)
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 测试用例
if __name__ == "__main__":
# 示例1：sin(x) 的 [0] Pade近似
print(pade_approximant("sin(x),x"))
# 1/(100*(x**2 + 20*x + 280))
# 示例2：sin(x) 的 [0] Pade近似
print(pade_approximant("sin(x),x,1"))
# (x**8/100 - x**7/29 + x**6/100 - 3*x**5/53 - 14*x**4/61 + 16*x**3/63 + 3*x**2/41 + 18*x/43 + 76/61)/(x**2 + 32*x/37 + 21/94)
padecoef 时滞的帕德逼近
时滞的帕德逼近
[num,den] = padecoef(T,N) 以传递函数形式返回持续时滞 exp(-T*s) 的 N 阶 Pade 逼近。行向量 num 和 den 包含了按 s 的降幂排序的分子和分母系数。分子和分母均为 N 次多项式。
T — 时滞,标量
N — 逼近的阶,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.special import comb
def pade_coef(input_str):
"""
计算时滞函数 e^(-sT) 的 N 阶 Padé 逼近系数
参数:
T (float): 时滞
N (int): 逼近的阶
返回:
a (np.ndarray): 分子多项式的系数（升幂排列）
b (np.ndarray): 分母多项式的系数（升幂排列）
"""
try:
expr = sp.sympify(input_str)
error = False
s = sp.symbols('s')
numerator = 0
denominator = 0
result = None, None
if isinstance(expr, tuple) and len(expr) == 2:
t, n = expr
if t.is_number and n.is_integer:
T = float(t)
N = int(n)
for k in range(N + 1):
# 计算组合数因子
c = comb(2 * N - k, N) / comb(2 * N, N)
# 分子累加项：(-T s)^k
numerator += c * (-T * s) ** k
# 分母累加项：(T s)^k
denominator += c * (T * s) ** k
# 提取多项式系数（SymPy按降幂排列，需反转）
num_poly = sp.poly(numerator, s)
den_poly = sp.poly(denominator, s)
# 转换为升幂排列的NumPy数组
a = np.array(num_poly.all_coeffs()[::-1], dtype=float)
b = np.array(den_poly.all_coeffs()[::-1], dtype=float)
result = a, b
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"解析错误: {str(e)}"
# 示例验证（N=2, T=2时的标准Padé系数）
result = pade_coef("2,2")
print("分子多项式系数（升幂排列）:", result[0])
# [ 1. -1. 0.66666667]
print("分母多项式系数（升幂排列）:", result[1])
# [1. 1. 0.66666667]
paramci 概率分布参数的置信区间
概率分布参数的置信区间
ci = paramci(data)返回数组ci,其中包含概率分布data中每个参数的95%置信区间的下限和上限.
data —— 输入,数值向量
ci —— 置信区间, 两维列表. 以p-by-2数组的形式返回, 分别是mu和sigma参数的下限和上限.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.stats as stats
def param_confidence_intervals(input_str):
"""
计算输入数据的正态分布参数置信区间。
参数:
input_str: 输入数据的字符串表示，可以是 SymPy 矩阵或列表。
confidence_level: 置信水平，默认为 0.95。
返回:
置信区间数组或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def matrix_to_list(matrix):
"""
将 SymPy 向量矩阵转换为 Python 列表。
参数:
matrix: SymPy 矩阵对象。
返回:
如果是行向量或列向量，返回对应的列表；否则返回 None。
"""
try:
if not isinstance(matrix, sp.Matrix):
raise ValueError("输入必须是一个 SymPy 矩阵对象")
rows, cols = matrix.shape
# 处理行向量（1行n列）
if rows == 1:
return list(matrix.row(0))
# 处理列向量（n行1列）
elif cols == 1:
return [matrix[i, 0] for i in range(rows)]
else:
return None # 非向量矩阵
except Exception as e:
print(f"Error in matrix_to_list: {e}")
return None
def evaluation_paramci(data, confidence_level): # 接受confidence_level参数
n = len(data)
if n < 2:
raise ValueError("数据长度必须至少为2")
mu = np.mean(data)
s = np.std(data, ddof=1)
# 使用t分布计算均值置信区间
t_value = stats.t.ppf((1 + confidence_level) / 2, df=n - 1)
mu_se = s / np.sqrt(n)
mu_ci = (mu - t_value * mu_se, mu + t_value * mu_se)
# 标准差置信区间保持不变...
alpha = 1 - confidence_level
chi2_upper = stats.chi2.ppf(1 - alpha / 2, df=n - 1)
chi2_lower = stats.chi2.ppf(alpha / 2, df=n - 1)
sigma_ci_lower = np.sqrt((n - 1) * s ** 2 / chi2_upper)
sigma_ci_upper = np.sqrt((n - 1) * s ** 2 / chi2_lower)
sigma_ci = (sigma_ci_lower, sigma_ci_upper)
return np.array([mu_ci, sigma_ci])
if isinstance(expr, list):
matrix = sp.Matrix(expr)
data_list = matrix_to_list(matrix)
if data_list is None:
return "错误：输入必须为向量（行向量或列向量）"
data = np.array(data_list, dtype=float)
result = evaluation_paramci(data, confidence_level=0.95) # 传递confidence_level
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例验证
data = [1, 2, 3, 4, 5]
input_str = str(data) # 输入字符串应为"[1, 2, 3, 4, 5]"
ci = param_confidence_intervals(input_str)
print("均值置信区间:", ci[0])
# [1.03675684 4.96324316]
print("标准差置信区间:", ci[1])
# [0.94731267 4.54349039]
ParametricPlot 参数方程图
参数方程图
ParametricPlot(x(u), y(u))通过独立参数（如u）控制多个坐标变量，能描述复杂轨迹、动态过程和几何变换。解决隐式方程难以绘制的问题。
尤其擅长表达: 闭合曲线, 自相交图形, 时间演化路径.
Lissajous 曲线 - 电子信号分析
ParametricPlot([sin(7u),sin(8u)], u=[0,2@pi])
卫星天线 - 垂直压缩椭圆
ParametricPlot(2.4cos(u),1.8sin(u),u=[0, 2@pi])
螺旋进给 - 机械加工路径
ParametricPlot(u*cos(u),u*sin(u),u=[0,8@pi])
心电信号 - 医学波形生成
ParametricPlot(16(sin(u)^3), 13cos(u)-5cos(2u)-2cos(3u)-cos(4u),u=[0,2@pi])
混沌吸引子 - 气象预测模型
ParametricPlot(sin(3u)*cos(5u), sin(2u)*sin(3u), u=[0,10@pi])
齿轮啮合 - 机械设计
ParametricPlot(cos(u)+0.2*cos(5u),sin(u)-0.2*sin(5u),u=[0,2@pi])
ParametricPlot3D 参数方程三维图
参数方程三维图
ParametricPlot3D(f1(u,v), f2(u,v), f3(u,v))使用两个参数（通常用 u 和 v 表示）来定义三维空间中的点 (x, y, z)。描述很多用显式方程 z = f(x, y) 难以甚至无法描述的复杂曲面
螺旋线 (Helix) - 弹簧、DNA结构、螺纹的基础: 这是最基本的空间曲线之一。用于建模弹簧、螺丝螺纹、DNA 双螺旋结构、螺旋楼梯的中心线、粒子在磁场中的轨迹等。
ParametricPlot3D(r*cos(t),r*sin(t),2t,r=[0,2],t=[0,4@pi])
球面 (Sphere) - 行星、球体、气泡: 这是描述完美球体的标准方式。用于建模行星、球状容器、气泡、原子模型、3D 图形中的球体等。
ParametricPlot3D(2*sin(r)*cos(t),2*sin(r)*sin(t), 2*cos(r),r=[0,@pi],t=[0,2@pi])
环面 (Torus) - 甜甜圈、轮胎、环形线圈: 形状像甜甜圈或游泳圈。用于建模轮胎、环形磁铁、粒子加速器（如环形对撞机）的真空管、一些建筑结构、首饰（如戒指的主体）等。
ParametricPlot3D((3+cos(v))*cos(u), (3+cos(v))*sin(u), sin(v),u=[0,2@pi],v=[0,2@pi])
柱面 (Cylinder) - 管道、柱子、罐体: 描述一个无限长或有限高的圆柱侧面。用于建模管道、柱子、罐体、滚筒、轴等。
ParametricPlot3D(2*cos(u), 2*sin(u),v,u=[0,2@pi],v=[0,5])
锥面 (Cone) - 漏斗、锥形零件、山峰简化模型：描述一个圆锥面。用于建模漏斗、锥形齿轮、交通锥筒、简化版的山峰或火山、灯具罩等。
ParametricPlot3D((5-v)*(3/5)*cos(u),(5-v)*(3/5)*sin(u),v,u=[0,2@pi],v=[0,5])
莫比乌斯带 (Möbius Strip) - 单侧曲面、有趣的拓扑结构、传送带设计：这是一个著名的单侧曲面（只有一个面，一条边界）。它具有有趣的拓扑性质。虽然实际工程应用较少（因其应力分布不均），但在传送带设计、艺术装置、数学教学中常被提及。
ParametricPlot3D((1+0.3v*cos(u/2))*cos(u), (1+0.3v*cos(u/2))*sin(u),0.3v*sin(u/2),u=[0,2@pi],v=[-1,1])
旋转曲面 (Surface of Revolution) - 花瓶、碗、车削零件： 这是生成轴对称曲面的主要方法。用于建模花瓶、碗、酒杯、车床加工的零件（轴、销、盘）、压力容器封头、旋转抛物面天线等。
ParametricPlot3D(sqrt(z)*cos(u), sqrt(z)*sin(u),z,u=[0,2@pi],z=[0,4])
螺旋面 (Helicoid) - 螺旋楼梯、螺纹、某些蛋白质结构：像一个“螺旋形的坡道”。用于建模螺旋楼梯的曲面、某些类型的螺纹（如方牙螺纹的简化模型）、DNA 结构的某些简化表示、螺旋坡道（如停车场）等。
ParametricPlot3D(v*cos(u),v*sin(u),u/2,u=[0,4@pi],v=[0,2])
双曲抛物面（马鞍面）: 建筑结构（如悉尼歌剧院）、重力场模拟
ParametricPlot3D(u,v,u^2-v^2,v=[-2,2], u=[-2,2])
恩涅珀曲面（极小曲面）: 肥皂膜自然形成的极小化曲面
ParametricPlot3D(u-u^3/3+u*v^2, v-v^3/3+v*u^2, u^2-v^2, v=[-2,2], u=[-2,2])
机械齿轮曲面: 齿轮齿面建模、螺纹设计
ParametricPlot3D((3+sin(5v))*cos(u), (3+sin(5v))*sin(u), v, v=[-1,1], u=[0,2@pi])
心脏形曲面（医学成像）: 模拟心脏瓣膜开闭时的血流动力学，用于心血管疾病研究和人工心脏设计
ParametricPlot3D((1-cos(u))*cos(u)*cos(v),(1-cos(u))*sin(u)*cos(v),(1-cos(u))*sin(v),u=[0,2@pi],v=[0,@pi])
空气动力学翼型（航空航天）: NACA 0012翼型方程, 飞机机翼和涡轮叶片的流体动力学模拟，优化升力/阻力比
ParametricPlot3D(u,(0.6*(0.2969*sqrt(u)-0.1260*u-0.3516*u^2+0.2843*u^3-0.1036*u^4))*cos(v),(0.6*(0.2969*sqrt(u)-0.1260*u-0.3516*u^2+0.2843*u^3-0.1036*u^4))*sin(v),v=[0,2@pi], u=[0,1])
DNA纳米管（生物技术）: 分子自组装结构设计，药物输送载体建模
ParametricPlot3D((3+0.3*cos(6u))*cos(v),(3+0.3*cos(6u))*sin(v),u+0.5*sin(6u),u=[0,4@pi],v=[0,2@pi])
重力透镜曲面（天体物理）: 模拟大质量天体引起的时空弯曲，解释星系引力透镜观测数据
ParametricPlot3D(u,v,0.5*exp(-(u^2 + v^2)/4)*cos(sqrt(u^2+v^2)),u=[-4,4],v=[-4,4])
混凝土薄壳建筑（土木工程）: 现代建筑薄壳屋顶设计（如悉尼歌剧院），优化承重分布
ParametricPlot3D(u,v,(u^2-v^2)/10,u=[-5,5],v=[-5,5])
量子位能阱曲面（量子计算）: 量子计算机中囚禁离子能级建模，优化量子位控制
ParametricPlot3D(u,v,-10/(1+u^2+v^2)+0.1*(u^2+v^2),u=[-4,4],v=[-4,4])
珊瑚生长模型（海洋生物学）: 海洋生态系统模拟，研究酸化对珊瑚骨架结构的影响
ParametricPlot3D((2+sin(3v))*u*cos(v),(2+sin(3v))*u*sin(v),0.5*u^2+2*sin(2v),u=[0,3],v=[0,2@pi])
超导磁悬浮曲面（物理）: 磁悬浮列车系统设计，优化超导体与永磁体的相互作用力
ParametricPlot3D(u,v,2*exp(-(u^2+v^2))*cos(2@pi*sqrt(u^2+v^2)),u=[-2,2],v=[-2,2])
ParametricPlot3DLine 参数方程三维曲线图
参数方程三维曲线图
ParametricPlot3DLine(x(t), y(t), z(t))通过参数方程定义三维空间中的连续曲线。
参数t变化时，点 (x,y,z)的轨迹形成曲线。这种表示法能精确描述复杂几何形状的运动轨迹。
1. 螺旋弹簧（机械工程）
应用：弹簧力学仿真
特征：恒定半径圆柱螺旋，z线性增长
ParametricPlot3DLine(3cos(t),3sin(t),0.5t,t=[0,8@pi])
2. DNA 双螺旋（生物物理）
应用：分子结构可视化
特征：两条相位差π的螺旋链
ParametricPlot3DLine(0.8cos(t+@pi),0.8sin(t+@pi),0.2t, t=[0,10@pi])
3. 飓风轨迹（气象学）
应用：风暴路径预测
特征：半径递增的不规则螺旋，模拟气压变化
ParametricPlot3DLine((1+0.3t)*cos(t+sin(2t)),(1+0.3t)*sin(t+cos(2t)),0.1t^2,t=[0,6@pi])
4. 黑洞吸积盘（天体物理）
应用：相对论性粒子轨道
特征：径向衰减 + 垂直振荡，模拟时空扭曲
ParametricPlot3DLine(2t/(1+t^2)*cos(10t),2t/(1+t^2)*sin(10t),t/5*exp(-0.2t)*sin(20t),t=[0,15])
5. 圆锥螺旋线（如粒子在磁场中的运动）
半径随高度增加而增大
ParametricPlot3DLine(t*cos(t),t*sin(t),t,t=[0,10@pi])
6. 环面结
Torus Knot，一种空间结，用于描述分子结构或装饰
ParametricPlot3DLine((2+cos(3t))*cos(2),(2+cos(3t))*sin(2),sin(3t),t=[0,2@pi])
7. 球面上的玫瑰曲线（如天线辐射方向图）
4瓣玫瑰曲线在球面上
ParametricPlot3DLine(sin(4t)*cos(t),sin(4t)*sin(t),cos(4t),t=[0,2@pi])
8. 利萨如曲线（Lissajous curve）
应用: 用于电子学、信号处理
ParametricPlot3DLine(sin(2t),sin(3t+@pi/4),sin(4t+@pi/3),t=[0,2@pi])
9. 心脏线在三维中的旋转
应用: 如艺术设计
ParametricPlot3DLine(16sin(t)^3,13cos(t)-5cos(2t)-2cos(3t)-cos(4t),t,t=[0,2@pi])
10. 空间导弹拦截轨迹
应用: 简化模型，假设匀速
ParametricPlot3DLine(t,t^2,sin(t),t=[0,2@pi])
11. 涡旋环（Vortex Ring，如烟雾环）中单个粒子的轨迹
ParametricPlot3DLine(cos(t)*(1-cos(t)),sin(t)*(1-cos(t)),sin(t),t=[0,2@pi])
12. 莫比乌斯带（Mobius strip）的中心线
ParametricPlot3DLine((1+0.5cos(t/2))*cos(t),(1+0.5cos(t/2))*sin(t),0.5sin(t/2),t=[0,2@pi])
13. 粒子加速轨道（高能物理）
应用：大型强子对撞机粒子轨迹
特征：径向衰减 + 高频振荡
ParametricPlot3DLine(2t/(1+t^2)*cos(10t),2t/(1+t^2)*sin(10t),t/5*exp(-0.2t)*sin(20t),t=[0,15])
14. 血管网络（医学成像）
应用：脑血管三维重建
特征：主血管缠绕 + 分支波动
ParametricPlot3DLine(3cos(t)+0.5cos(7t),3sin(t)+0.5sin(7t),2sin(3t),t=[0,4@pi])
15. 三维玫瑰曲线
xy平面：投影为旋转的4瓣花瓣
z方向：呈波浪形起伏（8个波峰波谷）
整体效果：类似弹簧上的旋转花朵
ParametricPlot3DLine(2cos(p)*cos(4p),2sin(p)*cos(4p),cos(4p)^2+@pi,p=[0,2@pi])
16. 生物运动：水母游动
三叶投影 → 水母伞状体的收缩舒张
z轴波动 → 触手的规律摆动
ParametricPlot3DLine(3sin(t)+2sin(3t),cos(t)-2cos(3t),cos(5t),t=[0,2@pi]
partfrac 部分分数分解
部分分数分解
partfrac(expr,var)找到expr相对于var的部分分数分解.如果不指定var,则partfrac使用默认变量.
partfrac(expr,var,Name,Value）使用由一个或多个Name,Value对参数指定的附加选项来查找部分分数分解.
目前只支持Name=FactorMode,Value=full, 尚不支持Value=complex or Value=real
expr —— 有理表达式,符号表达,符号功能.
var —— 感兴趣的变量,符号变量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Poly
from sympy.abc import x
def partial_fraction_decomposition(input_str):
"""
执行有理函数的部分分式分解
参数:
input_str (str): 支持以下格式：
1. 单独有理式："(x+1)/(x**2-1)"
2. 带变量声明："( (x+1)/(x**2-y^2), y )"
3. 带分解模式："( (x+1)/(x**3-1), x, full=True )"
返回:
Expr/str: 分解后的表达式或错误信息
示例:
>>> partial_fraction_decomposition("(x+1)/(x**2-1)")
1/(x - 1)
>>> partial_fraction_decomposition("( (x+y)/(x**2 - y**2), y )")
1/(x - y)
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str, evaluate=False)
# 初始化参数
func = None
var = x # 默认使用x作为变量
full_decompose = False
# 参数类型处理
if isinstance(expr, tuple):
# 提取主表达式
func = expr[0]
# 处理变量声明
if len(expr) >= 2:
var = expr[1]
if not var.is_Symbol:
return "变量声明错误: 第二个参数必须是符号变量"
# 处理分解模式
if len(expr) >= 3:
if expr[2] == sp.S.true or expr[2] == True:
full_decompose = True
elif expr[2] == sp.S.false or expr[2] == False:
full_decompose = False
else:
return "模式参数错误: 第三个参数应为True/False"
elif isinstance(expr, sp.Expr):
func = expr
else:
return "输入格式错误: 需要有理表达式或参数元组"
# 验证表达式类型
if not isinstance(func, sp.Expr):
return "表达式类型错误: 输入必须是有理函数"
# 执行分解
result = sp.apart(func, var, full=full_decompose)
# 验证分解结果有效性
if result == func:
return "分解失败: 无法进行部分分式分解"
return result
except sp.SympifyError:
return "表达式解析失败: 请检查输入语法"
except TypeError as te:
return f"类型错误: {str(te)}"
except ValueError as ve:
return f"数值错误: {str(ve)}"
except Exception as e:
return f"未知错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 有效输入测试
y = sp.symbols('y')
test_cases = [
("简单分解", "(x+1)/(x**2-1)"),
# 1/(x - 1)
("多变量分解", "( (x+y)/(x**2 - y^2), y )"),
# 1/(x - y)
("高次分解", "( (x+1)/(x**3-1), x, True )"),
# RootSum(_w**2 + _w + 1, Lambda(_a, (_a**2/3 + _a/3)/(-_a + x))) + 2/(3*(x - 1))
("假分式处理", "(x^3 + 1)/(x^2 - 1)")
# x + 1/(x - 1)
]
for name, input_str in test_cases:
print(f"{name}测试:")
result = partial_fraction_decomposition(input_str)
if isinstance(result, sp.Expr):
print(f"输入: {input_str}")
print(f"结果: {result}")
else:
print(f"错误结果: {result}\n")
parzenwin parzen窗
parzen窗
w = parzenwin(L,sym=1) 返回一个长度为L个点的修正的parzen窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.signal as signal
def parzenwin_window(input_str):
"""
生成Parzen窗函数，对标MATLAB的parzenwin函数
参数:
input_str (str):
- 单独整数: 生成对称Parzen窗（默认行为，对应MATLAB parzenwin(N)）
- 元组格式: (N, sym_flag)
N: 窗口长度（正整数）
sym_flag: 对称标志（0: 非对称周期性窗口，1: 对称窗口，默认1）
返回:
list: 窗函数数值列表 或 错误信息字符串
示例:
>>> parzenwin_window("5")
[0.0, 0.25, 1.0, 0.25, 0.0] (对称窗)
>>> parzenwin_window("(5, 0)")
[0.0, 0.0625, 1.0, 0.0625, 0.0] (非对称窗)
"""
try:
expr = sp.sympify(input_str)
N = None
sym_flag = 1 # 默认对称窗，对齐MATLAB
# 参数解析分支
if isinstance(expr, tuple):
# 验证元组格式 (N, sym_flag)
if len(expr) != 2 or not expr[0].is_Integer or not expr[1].is_Integer:
return f"输入格式错误: 需要 (整数, 0/1)，实际输入 {input_str}"
N = int(expr[0])
sym_flag = int(expr[1])
# 验证参数范围
if N <= 0:
return f"窗口长度必须为正整数，当前输入 N={N}"
if sym_flag not in (0, 1):
return f"对称标志必须为0或1，当前输入 sym_flag={sym_flag}"
elif expr.is_Integer:
# 单独整数输入
N = int(expr)
if N <= 0:
return f"窗口长度必须为正整数，当前输入 N={N}"
else:
return f"输入类型错误: 需要整数或元组，实际输入 {type(expr)}"
# 生成窗函数（sym=True对应MATLAB对称行为）
window = signal.windows.parzen(N, sym=(sym_flag == 1))
# 转换为列表并保留6位小数（对齐MATLAB显示精度）
return [round(float(x), 6) for x in window]
except Exception as e:
return f"错误: {str(e)}"
if __name__ == "__main__":
# ----------------- 测试用例 -----------------
def print_test_case(input_str, comment):
print(f"\n{comment}: {input_str}")
result = parzenwin_window(input_str)
print(result if isinstance(result, str) else [round(x, 4) for x in result])
# 测试1: 对称窗（MATLAB默认行为）
print_test_case("5", "5点对称窗")
# [0.016, 0.424, 1.0, 0.424, 0.016]
# 测试2: 非对称窗
print_test_case("(5, 0)", "5点非对称窗")
# [0.0093, 0.25, 0.8611, 0.8611, 0.25]
pascal 帕斯卡矩阵
帕斯卡矩阵
P = pascal(n) 返回n阶帕斯卡矩阵. P是一个对称正定矩阵,其整数项来自帕斯卡三角形.P的逆矩阵具有整数项.
P = pascal(n,1) 返回帕斯卡矩阵的下三角乔列斯基因子(最高到列符号).P是对合矩阵,即该矩阵是它自身的逆矩阵.
P = pascal(n,2) 返回 pascal(n,1) 的转置和置换版本.在这种情况下,P是单位矩阵的立方根.
n - 是矩阵的阶数,非负整数标量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.linalg import pascal
def pascal_type_matrix(input_str):
"""
生成帕斯卡矩阵（Pascal Matrix），对标MATLAB的pascal函数
参数:
input_str (str):
- 单独整数: 生成对称帕斯卡矩阵 (n x n)
- 元组格式: (n, type)
type=1: 下三角Cholesky分解矩阵
type=2: 上三角转置矩阵
返回:
sp.Matrix: SymPy矩阵对象 或 错误信息字符串
示例:
>>> pascal_type_matrix("3")
Matrix([
[1, 1, 1],
[1, 2, 3],
[1, 3, 6]])
>>> pascal_type_matrix("(3, 1)")
Matrix([
[1, 0, 0],
[1, 1, 0],
[1, 2, 1]])
"""
try:
expr = sp.sympify(input_str)
result = None
# 情况1: 输入为元组 (n, type)
if isinstance(expr, tuple):
# 验证元组长度和元素类型
if len(expr) != 2 or not all(e.is_integer for e in expr):
return f"输入格式错误: 需要 (整数, 1/2)，实际输入 {input_str}"
n = int(expr[0]) # 转换为Python整数
kind_flag = int(expr[1])
# 根据类型标志生成矩阵
if kind_flag == 1:
arr = pascal(n, kind='lower')
elif kind_flag == 2:
arr = pascal(n, kind='upper')
else:
return f"类型参数错误: 仅支持1或2，实际输入 {kind_flag}"
# 情况2: 输入为单个整数
elif expr.is_integer:
n = int(expr)
arr = pascal(n) # 默认对称矩阵
else:
return f"输入类型错误: 需要整数或元组，实际输入 {type(expr)}"
# 将NumPy数组转换为SymPy矩阵
return sp.Matrix(arr.tolist())
except Exception as e:
return f"错误: {e}"
if __name__ == "__main__":
# ----------------- 测试用例 -----------------
def print_test_case(input_str, comment):
print(f"\n{comment}: {input_str}")
result = pascal_type_matrix(input_str)
print(result if isinstance(result, str) else result.applyfunc(lambda x: x))
# 测试1: 3x3对称矩阵（对应MATLAB pascal(3)）
print_test_case("3", "对称帕斯卡矩阵")
# Matrix([[1, 1, 1],
# [1, 2, 3],
# [1, 3, 6]])
# 测试2: 3阶下三角矩阵（对应MATLAB pascal(3,1)）
print_test_case("(3, 1)", "下三角分解矩阵")
# Matrix([[1, 0, 0],
# [1, 1, 0],
# [1, 2, 1]])
# 测试3: 3阶上三角矩阵（对应MATLAB pascal(3,2)）
print_test_case("(3, 2)", "上三角转置矩阵")
# Matrix([[1, 1, 1],
# [0, 1, 2],
# [0, 0, 1]])
pcg 求解线性系统 - 预条件共轭梯度法
求解线性系统 - 预条件共轭梯度法
x,flog = pcg(A,b) 尝试使用预条件共轭梯度法求解关于 x 的线性系统 A*x = b. 如果尝试成功, pcg 会显示一条消息来确认收敛. 如果 pcg 无法在达到最大迭代次数后收敛或出于任何原因暂停, 则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代序号的诊断消息.
x,flog = pcg(A,b,tol) 指定该方法的容差. 默认容差是 1e-6.
x,flog = pcg(A,b,tol,maxit) 指定要使用的最大迭代次数.如果 pcg 无法在 maxit 次迭代内收敛，将显示诊断消息.
x,flog = pcg(A,b,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解关于 y 的方程组 H−1AH−Ty=H−1b 来计算 x, 其中 y=HTx 且 H=M1/2=(M1M2)1/2. 该算法不显式形成 H.使用预条件子矩阵可以改善问题的数值属性和计算的效率.
A — 系数矩阵,矩阵
b — 线性方程的右侧,列向量
tol — 方法容差, [] 或 1e-6 (默认), 正标量
maxit — 最大迭代次数, [] 或 min(size(A,1),20) (默认) , 正整数标量
M — 预条件子矩阵（以单独参量指定）
x — 线性系统的解, 列向量
flag — 收敛标志, 标量
flag = 0 : 成功 - pcg 在 maxit 次迭代内收敛至所需容差 tol
flag = 1 : 失败 - pcg 执行了 maxit 次迭代，但未收敛
flag = 2 : 失败 - 预条件子矩阵 M 或 M = M1*M2 为病态
flag = 3 : 失败 - pcg 在经过两次相同的连续迭代后已停滞
flag = 4 : 失败 - 由 pcg 算法计算的标量数量之一变得太小或太大，无法继续计算
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.sparse.linalg import cg
import ast
def pcg_solver(A, b, tol=1e-6, maxiter=1000, M=None):
x, info = cg(A, b, tol=tol, maxiter=maxiter, M=M)
return x, info
def preconditioned_conjugate_gradient(input_str):
"""
解析输入字符串并执行PCG求解，输入格式对标MATLAB的pcg函数。
输入字符串格式:
"(A, b, tol, maxiter, M)" 其中：
- A: SymPy格式矩阵（例如'Matrix([[1,2],[3,4]])'）
- b: SymPy格式向量（例如'[5,6]'）
- tol, maxiter, M: 可选参数
返回:
(解向量, 状态码) 或错误信息
"""
try:
expr = ast.literal_eval(input_str)
if not isinstance(expr, tuple) or len(expr) < 2:
return f"输入错误: 需要至少(A, b)，当前输入为 {input_str}"
# 转换矩阵 A
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
if A_sym is None:
return f"矩阵A解析错误: {expr[0]}"
A_np = np.array(A_sym.tolist(), dtype=float)
# 转换向量 b
b_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
if b_sym is None or b_sym.shape[1] != 1: # 检查是否为列向量
return f"向量b解析错误: {expr[1]}"
b_np = np.array(b_sym.tolist(), dtype=float).flatten()
# 检查维度一致性
if A_np.shape[0] != b_np.shape[0]:
return f"维度不匹配: A的行数({A_np.shape[0]}) ≠ b的长度({b_np.shape[0]})"
# 解析可选参数
tol, maxiter, M_np = 1e-6, 1000, None
if len(expr) >= 3:
tol = float(expr[2])
if len(expr) >= 4:
tol = float(expr[2])
maxiter = int(expr[3])
if len(expr) >= 5:
tol = float(expr[2])
maxiter = int(expr[3])
M_sym = sp.Matrix(expr[4]) if isinstance(expr[4], list) else None
if M_sym is None or M_sym.shape != A_sym.shape:
return f"预条件矩阵M解析错误: {expr[4]}"
M_np = np.array(M_sym.tolist(), dtype=float)
# 调用求解器
x, info = pcg_solver(A=A_np, b=b_np, tol=tol, maxiter=maxiter, M=M_np)
return (x.tolist(), info)
except Exception as e:
return f"错误: {str(e)}"
# ------------------------- 示例测试 -------------------------
if __name__ == "__main__":
# 示例1: 基本用法（2x2系统）
input_str1 = "([[4,1],[1,3]]), [1, 2]"
result1 = preconditioned_conjugate_gradient(input_str1)
print("示例1 解:", result1[0], "状态码:", result1[1])
# 解: [0.09090909090909091, 0.6363636363636364] 状态码: 0
# 示例2: 使用预条件矩阵
input_str2 = "([[25,15,-5],[15,18,0],[-5,0,11]], [1, 2, 3], 1e-8, 1000, [[5,0,0],[0,6,0],[0,0,3]])"
result2 = preconditioned_conjugate_gradient(input_str2)
print("\n示例2 解:", result2[0], "状态码:", result2[1])
# 解: [0.06814814814814821, 0.054320987654320946, 0.30370370370370375] 状态码: 0
pchip 分段三次Hermite插值多项式
分段三次Hermite插值多项式
p = pchip(x,y,xq) 返回与 xq 中的查询点对应的插值向量 p。p 的值由 x 和 y 的保形分段三次插值确定。
pp = pchip(x,y) 返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp。
x — 样本点, 向量
y — 样本点处的函数值, 向量 | 矩阵
xq — 查询点, 标量 | 向量 | 矩阵
p — 查询点位置的插值, 标量 | 向量 | 矩阵
pp — 分段多项式, 结构体
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.interpolate import PchipInterpolator
from sympy import sympify, lambdify, pprint
def pchip_interpolation(input_str):
"""
实现MATLAB风格的pchip插值，支持数值计算和符号表达式
参数格式：
- 模式1（返回插值函数）："(x, y)"
- 模式2（直接计算插值）："(x, y, xx)"
返回：
- 模式1：可调用函数
- 模式2：NumPy数组
- 错误时返回字符串描述
示例：
>>> pchip_interpolation("([1,2,3,4], [2,1,4,3])")
>>> pchip_interpolation("([1,2,3,4], [2,1,4,3], [1.5, 2.5])")
array([1.6875, 2.375 ])
"""
try:
# 解析输入
expr = sympify(input_str)
# 参数数量验证
if not isinstance(expr, tuple) or len(expr) not in (2, 3):
return "输入错误：需要(x, y) 或 (x, y, xx)格式"
# 提取数据点
x = np.array(expr[0], dtype=float)
y = np.array(expr[1], dtype=float)
# 数据校验
if x.shape != y.shape:
return "错误：x和y维度不匹配"
if len(x) < 2:
return "错误：至少需要两个数据点"
if not np.all(np.diff(x) > 0):
return "错误：x必须单调递增"
# 创建插值器
pchip = PchipInterpolator(x, y)
# 模式1：返回可调用函数
if len(expr) == 2:
return pchip_symbolic(x, y)
# 模式2：直接计算结果
if len(expr) == 3:
xx = np.array(expr[2], dtype=float)
return pchip(xx)
except Exception as e:
return f"错误：{str(e)}"
# 增强版符号表达式生成（可选功能）
def pchip_symbolic(x_nodes, y_nodes):
"""
生成符号表达式形式的PCHIP插值函数
参数：
x_nodes - 节点x值（列表或数组）
y_nodes - 节点y值（列表或数组）
返回：
SymPy Piecewise表达式
"""
x = sp.symbols('x')
pchip = PchipInterpolator(x_nodes, y_nodes)
# 获取分段多项式系数
breaks = pchip.x
coeffs = pchip.c
pieces = []
for i in range(len(breaks) - 1):
a, b = breaks[i], breaks[i + 1]
c = coeffs[:, i]
poly = sum(c[k] * (x - a) ** (3 - k) for k in range(4))
pieces.append((poly, (x >= a) & (x <= b)))
pieces.append((sp.nan, True)) # 超出范围返回nan
return sp.Piecewise(*pieces)
# 测试案例
if __name__ == "__main__":
# 测试模式1：返回插值函数
print("测试1：创建插值函数")
interp_func = pchip_interpolation("([1, 2, 3, 4], [2, 1, 4, 3], 2.5)")
print("类型检查:", type(interp_func).__name__) # 应显示 'function'
print("f(2.5) =", interp_func)
# 2.5
# 测试模式2：直接计算
print("\n测试2：直接插值计算")
result = pchip_interpolation("([1, 2, 3, 4], [2, 1, 4, 3], [1.5, 2.5, 3.5])")
print("插值结果:", result)
# [1.125 2.5 3.875]
# 测试模式3：返回结构体
print("\n测试3：返回结构体")
result = pchip_interpolation("([1, 2, 3, 4], [2, 1, 4, 3])")
print(result)
# Piecewise((-3.0*x - 1.0*(x - 1.0)**3 + 3.0*(x - 1.0)**2 + 5.0, (x >= 1.0) & (x <= 2.0)),
# (-48.0*(0.5*x - 1)**3 + 36.0*(0.5*x - 1)**2 + 1.0, (x >= 2.0) & (x <= 3.0)),
# (4.0 - 27.0*(0.333333333333333*x - 1)**3, (x >= 3.0) & (x <= 4.0)), (nan, True))
pDiff 复合函数偏导数
复合函数偏导数
p = pDiff(f(x),[vars],[expressions],var)
f(x) -- 符号表达式,符号函数.
vars -- 依赖变量列表.
expressions -- 依赖方程列表.
var -- 求偏导数的独立变量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def partial_derivatives(input_str):
"""
计算复合函数的偏导数，对标 mathstudio 的 pDiff 函数。
参数:
input_str (str): 输入字符串，格式为四元组 '(f, (原变量), (新变量), 目标变量)'
其中 f 是原函数，原变量和新变量为符号变量或表达式，目标变量为求导变量。
返回:
sp.Expr 或 str: 偏导数表达式或错误信息。
示例:
>>> partial_derivatives('(x**2, (x,), (u,), u)')
2*u
>>> partial_derivatives('(x*y, (x,y), (u+v, u-v), u)')
2*u
数值计算示例:
expr = partial_derivatives('(x*y, (x,y), (u+v, u-v), u)')
expr.subs({'u': 2, 'v': 3}) # 得到数值结果 4
"""
try:
# 将输入字符串转换为 Sympy 表达式（自动解析元组结构）
expr = sp.sympify(input_str)
# 检查输入是否为四元组（函数, 原变量列表, 新变量列表, 目标变量）
if not (isinstance(expr, tuple) and len(expr) == 4):
return f"输入错误: 输入应为四元组，实际输入为 {input_str}"
f, independent_vars, new_vars, target_var = expr
# 检查原变量和新变量数量是否一致
if len(independent_vars) != len(new_vars):
return f"变量数量不匹配: 原变量{len(independent_vars)}个，新变量{len(new_vars)}个"
# 创建变量替换字典（例如 {x: u+v, y: u-v}）
substitution = {old: new for old, new in zip(independent_vars, new_vars)}
f_substituted = f.subs(substitution)
# 计算替换后的函数对目标变量的偏导数
derivative = sp.diff(f_substituted, target_var)
return derivative
except Exception as e:
return f"错误: {e}"
# 示例1：基本使用（符号计算）
input_str1 = '(x**2, (x,), (u,), u)'
result1 = partial_derivatives(input_str1)
print("符号导数:", result1)
# 2*u
# 示例2：复合函数求导
input_str2 = '(x*y, (x,y), (u+v, u-v), u)'
result2 = partial_derivatives(input_str2)
print("复合函数导数:", result2)
# 2*u
# 示例3：数值计算
u_value = 2
v_value = 3
numerical_result = result2.subs({'u': u_value, 'v': v_value})
print(f"当 u={u_value}, v={v_value} 时，导数值为:", numerical_result)
# 4
pdist 成对观测值之间的两两距离
成对观测值之间的两两距离
D = pdist(X,Distance,DistParameter) 返回 X 中成对观测值之间的距离.
X是数值矩阵
Distance是距离类型,默认是'euclidean'(欧几里德)
其他距离类型包括euclidean,canberra,cityblock,chebyshev,minkowski,hamming,cosine,correlation,jaccard
DistParameter,正整数标量,当Distance类型是minkowski距离时候有效,默认为2
D是两两距离,数值向量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.spatial import distance
def distance_paired_first(input_str):
"""
计算成对观测值之间的两两距离，对标MATLAB的pdist函数。
参数:
input_str: 输入字符串，格式为"(data, metric, ...)"，其中:
- data: 二维数据矩阵，每行代表一个观测值
- metric: 距离度量类型（如'euclidean', 'minkowski'等）
- 后续参数: 如minkowski距离的p值
返回:
距离列表或错误信息。
"""
try:
expr = sp.sympify(input_str)
points = None
distance_type = 'euclidean' # 默认距离类型
p = 2 # Minkowski距离的默认p值
# 解析输入参数
if isinstance(expr, tuple):
# 检查数据部分
if len(expr) < 1:
return "错误: 输入必须包含数据矩阵"
data_expr = expr[0]
# 转换数据为矩阵
X = sp.Matrix(data_expr) if isinstance(data_expr, list) else None
if X is None:
return "错误: 无法解析数据矩阵"
points = np.array(X.tolist(), dtype=float)
# 检查数据维度并处理一维数据
if points.ndim == 1:
points = points.reshape(-1, 1)
elif points.ndim != 2:
return "错误: 数据必须是二维数组"
# 解析距离类型
if len(expr) >= 2:
distance_type = str(expr[1]).lower()
# 检查是否支持的距离类型
supported_metrics = ["euclidean", "canberra", "cityblock", "chebyshev",
"minkowski", "hamming", "cosine", "correlation", "jaccard"]
if distance_type not in supported_metrics:
return f"错误: 不支持的距离类型 '{distance_type}'"
# 处理Minkowski的p参数
if distance_type == 'minkowski' and len(expr) >= 3:
p = float(expr[2]) # 允许浮点数
else:
distance_type = 'euclidean'
else:
# 输入不是元组，尝试解析为数据矩阵
X = sp.Matrix(expr) if isinstance(expr, list) else None
if X is None:
return "错误: 输入格式不正确，应为 (data, metric, ...)"
points = np.array(X.tolist(), dtype=float)
# 检查数据维度并处理一维数据
if points.ndim == 1:
points = points.reshape(-1, 1)
elif points.ndim != 2:
return "错误: 数据必须是二维数组"
# 计算距离
if distance_type == 'minkowski':
distances = distance.pdist(points, metric=distance_type, p=p)
else:
distances = distance.pdist(points, metric=distance_type)
return distances.tolist()
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例1：二维数据的欧氏距离
example1 = '([[1, 2], [3, 4]], "euclidean")'
print("示例1结果:", distance_paired_first(example1))
# [2.8284271247461903]
# 示例2：一维数据转换为二维后的曼哈顿距离
example2 = '([1, 2, 3, 4], "cityblock")'
print("示例2结果:", distance_paired_first(example2))
# [1.0, 2.0, 3.0, 1.0, 2.0, 1.0]
# 示例3：带p参数的Minkowski距离
example3 = '([[1, 2], [3, 4]], "minkowski", 3)'
print("示例3结果:", distance_paired_first(example3))
# [2.5198420997897464]
pdist2 两组观测值之间的两两距离
两组观测值之间的两两距离
D = pdist(X,Y,Distance,DistParameter) 返回返回 X 和 Y 中每对观测值之间的距离.
X,Y是数值矩阵
Distance是距离类型,默认是'euclidean'(欧几里德),
其他距离类型包括euclidean,canberra,cityblock,chebyshev,minkowski,hamming,cosine,correlation,jaccard
DistParameter,正整数标量,当Distance类型是minkowski距离时候有效,默认为2
D是两两距离,数值矩阵. 是一个压缩形式的距离矩阵(通过squareform还原完整矩阵).
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.spatial import distance
def distance_paired_second(input_str):
"""
计算两组观测值之间的两两距离，对标MATLAB的pdist2函数。
参数:
input_str (str): 输入字符串，格式为'(list1, list2, metric, p)'，其中：
list1: 第一组观测值，可以是二维列表或一维列表。
list2: 第二组观测值，可以是二维列表或一维列表。
metric (str, 可选): 距离类型，支持的类型包括：
'euclidean', 'canberra', 'cityblock', 'chebyshev', 'minkowski',
'hamming', 'cosine', 'correlation', 'jaccard'。默认为 'euclidean'。
p (int, 可选): 仅当metric为'minkowski'时有效，指定闵可夫斯基距离的p值，默认为2。
返回:
sympy.Matrix: 距离矩阵，元素(i,j)表示list1的第i个观测值与list2的第j个观测值的距离。
如果输入有误，返回错误信息字符串。
示例:
>>> distance_paired_second("([[1, 2], [3, 4]], [[5, 6], [7, 8]], 'euclidean')")
Matrix([
[5.830951894845301, 8.246211251235321],
[2.8284271247461903, 5.830951894845301]])
"""
try:
expr = sp.sympify(input_str)
if isinstance(expr, tuple) and len(expr) >= 2 and isinstance(expr[0], (list, sp.Matrix)) and isinstance(expr[1],
(list,
sp.Matrix)):
# 转换为NumPy数组
list1 = np.array(expr[0], dtype=float)
list2 = np.array(expr[1], dtype=float)
# 确保处理为二维数组
if list1.ndim == 1:
list1 = list1.reshape(1, -1)
if list2.ndim == 1:
list2 = list2.reshape(1, -1)
# 检查观测值的维度是否一致
if list1.shape[1] != list2.shape[1]:
return f"输入错误: 观测值的维度不一致 {input_str}"
# 处理距离类型
distance_type = 'euclidean'
valid_metrics = ["euclidean", "canberra", "cityblock", "chebyshev", "minkowski",
"hamming", "cosine", "correlation", "jaccard"]
if len(expr) >= 3:
distance_type = str(expr[2]).lower()
if distance_type not in valid_metrics:
return f"输入错误: 不支持的距离类型 {distance_type}"
# 处理Minkowski距离的p参数
p_param = 2
if distance_type == 'minkowski':
if len(expr) >= 4:
try:
p_param = int(expr[3])
except:
return f"输入错误: Minkowski距离的p参数必须为整数 {expr[3]}"
# 计算距离矩阵
if distance_type == 'minkowski':
dist_matrix = distance.cdist(list1, list2, distance_type, p=p_param)
else:
dist_matrix = distance.cdist(list1, list2, distance_type)
return sp.Matrix(dist_matrix)
else:
return f"输入错误: 输入格式不正确 {input_str}"
except Exception as e:
return f"错误: {e}"
input_str = "([[1, 2], [3, 4]], [[5, 6], [7, 8]], 'euclidean')"
result = distance_paired_second(input_str)
print(result)
# Matrix([[5.65685424949238, 8.48528137423857],
# [2.82842712474619, 5.65685424949238]])
input_str = "([1, 2], [3, 4], 'cityblock')"
result = distance_paired_second(input_str)
print(result)
# Matrix([[4.00000000000000]])
input_str = "([[1, 2], [3, 4]], [[5, 6]], 'minkowski', 1)"
result = distance_paired_second(input_str)
print(result)
# Matrix([[8.00000000000000],
# [4.00000000000000]])
perms 所有可能的排列
所有可能的排列
P = perms(v) 返回的矩阵包含了向量v中元素的所有排列.P的每一行包含v中n个元素的一个不同排列.矩阵P具有与v相同的数据类型,包含n!行和n列.
v — 项目集合, 数值,logical或char值的向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import itertools
def possible_all_arrangements(input_str):
"""
对标 MATLAB 的 perms 函数，生成所有可能的排列（反向字典序）
参数:
input_str: 输入数据的字符串表示，可以是列表或矩阵形式，例如 "[1,2,3]" 或 "Matrix([1,2,3])"
返回:
SymPy 矩阵对象，每行代表一种排列，排列顺序与 MATLAB 一致
错误时返回字符串描述
"""
try:
# 若失败则尝试解析为 SymPy 矩阵
expr = sp.sympify(input_str)
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is None:
return f"输入格式错误: {input_str}"
# 检查矩阵维度
rows, cols = matrix.shape
if rows == 1: # 行向量
elements = list(matrix.row(0))
elif cols == 1: # 列向量
elements = [matrix[i, 0] for i in range(rows)]
else:
return "错误：输入必须是行向量或列向量"
# 生成所有排列并反转顺序以匹配 MATLAB
permutations = list(itertools.permutations(elements))[::-1] # 关键步骤：反转顺序
# 转换为 SymPy 矩阵（每行一个排列）
return sp.Matrix(permutations)
except Exception as e:
return f"错误:{e}"
# 示例验证 ==============================================
if __name__ == "__main__":
# 测试用例 1：标准输入
input_str1 = "[2, 4, 6]"
result1 = possible_all_arrangements(input_str1)
print("测试1 输入:", input_str1)
print("输出排列矩阵:\n", result1)
# Matrix([[6, 4, 2],
# [6, 2, 4],
# [4, 6, 2],
# [4, 2, 6],
# [2, 6, 4],
# [2, 4, 6]])
# 测试用例 2：SymPy 矩阵输入
input_str2 = "[[5], [4], [9]]" # 列向量
result2 = possible_all_arrangements(input_str2)
print("\n测试2 输入:", input_str2)
print("输出排列矩阵:\n", result2)
# Matrix([[9, 4, 5],
# [9, 5, 4],
# [4, 9, 5],
# [4, 5, 9],
# [5, 9, 4],
# [5, 4, 9]])
permute 置换数组维度
置换数组维度
B = permute(A,dimorder) 按照向量 dimorder 指定的顺序重新排列数组的维度。例如，permute(A,[2 1]) 交换矩阵 A 的行和列维度。通常，输出数组的第 i 个维度是输入数组的维度 dimorder(i)。
A — 输入数组, 向量 | 矩阵
dimorder — 维度顺序, 行向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def permutation_order_vector(input_str):
"""
实现类似 MATLAB 的 permute 函数，处理向量和二维矩阵的维度排列。
参数格式要求:
input_str 应为有效表示矩阵和维度的字符串，例如："([[1,2,3]], [2])" 或 "([[1,2],[3,4]], [2,1])"
返回:
重新排列后的矩阵 或 错误信息字符串
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def numpy_permute(A, dimorder):
"""
MATLAB permute函数的NumPy实现
参数:
A: 输入数组 (NumPy ndarray)
dimorder: 维度顺序列表 (基于1的索引，如MATLAB)
返回:
维度置换后的数组
"""
# 将MATLAB的1-based索引转换为NumPy的0-based索引
axes = [d - 1 for d in dimorder]
# 使用NumPy的transpose函数
return np.transpose(A, axes=axes)
def sympy_permute(A, dimorder):
"""
MATLAB permute函数的SymPy实现
参数:
A: 输入数组 (SymPy Array 或可转换为Array的对象)
dimorder: 维度顺序列表 (基于1的索引，如MATLAB)
返回:
维度置换后的SymPy Array
"""
# 验证维度数量
ndim = A.rank()
if len(dimorder) != ndim:
raise ValueError(f"维度顺序列表长度({len(dimorder)})必须等于数组维度数({ndim})")
# 将MATLAB的1-based索引转换为0-based索引
axes = [d - 1 for d in dimorder]
# 执行维度置换
return sp.permutedims(A, axes)
if isinstance(expr, tuple) and len(expr) == 2:
sp_array = sp.Array(expr[0]) if isinstance(expr[0], list) else None
if sp_array is None:
raise ValueError(f"输入错误:{input_str}")
if not isinstance(expr[1], list):
raise ValueError(f"输入错误:{input_str}")
if sp_array.free_symbols:
result = sympy_permute(sp_array, expr[1])
else:
z_arr = np.asarray(sp_array)
result = numpy_permute(z_arr, expr[1])
else:
error = True
return result if not error else f"输入错误:{input_str}"
except Exception as e:
return f"运行时错误：{str(e)}"
# 示范数值数组
print(permutation_order_vector("[1,2,3],[1]"))
# [1 2 3]
# 示范符号数组
print(permutation_order_vector("[[[x,y],[z,x+y]],[[y,z],[x,x*y*z]]],[3,1,2]"))
# [[[x, z],
# [y, x]],
# [[y, x + y],
# [z, x*y*z]]]
phaseplot 相频图
相频图
phaseplot(sys) 显示系统输出信号相对于输入信号的相位偏移（以度为单位）随频率变化的特性.
低通滤波器（一阶系统）: RC电路噪声滤除。
确认系统是低通特性（相位滞后随频率增加），并帮助估计截止频率附近（约10 rad/s）的相位滞后.
phaseplot(1,0.1s+1)
直流电机速度控制: 电机调速系统
有助于识别系统类型（含积分器，低频相位为-90°），并评估稳定性风险.
phaseplot(5,0.5*s^2+s)
质量-弹簧-阻尼系统: 机械振动分析
分析机械振动中的共振现象和相位反转.
phaseplot(1,s^2+0.6s+1)
超前补偿器（PD控制器）: 提高系统稳定裕度
体现了PD控制器提升稳定裕度的能力.
phaseplot(2s+10,0.02s+1)
带通滤波器: 通信信号选频
识别带通特性（中心频率约10 rad/s）和群延迟特性。
phaseplot(0.1s,(s+1)(0.01s+1))
RC低通滤波器
ADC采样前的抗混叠滤波
频率范围 相位偏移
0.1 rad/s -5.7°
10 rad/s -45°
100 rad/s -84.3°
phaseplot([1],[0.1,1])
伺服电机位置控制
工业机器人关节定位精度保障
频率范围 相位偏移
1 rad/s -11.3°
7.07 rad/s -90° (共振点)
50 rad/s -172.9°
phaseplot([50],[1,5,50])
心电图工频陷波器
医疗设备中消除50 Hz电源干扰
频率范围 相位偏移
40 Hz +87°
50 Hz 0° (相位跳变)
60 Hz -87°
phaseplot([1,0,10000],[1,100,10000])
飞机俯仰控制系统
民航客机自动驾驶模式稳定性验证
频率范围 相位偏移
0.5 rad/s -28°
5 rad/s -112° (增益穿越频率)
50 rad/s -256°
phaseplot([0.5,1],[0.01,0.1,1,1])
5G毫米波相位补偿
Massive MIMO基站波束成形校准
频率范围 相位偏移
1e9 rad/s +4.3°
5e10 rad/s +35.5° (工作频段)
1e11 rad/s +52.1°
phaseplot([0.02,1],[0.002,1])
电动汽车电池管理系统
锂电池健康状态(SOH)在线诊断
频率范围 相位偏移
10 rad/s -21.8°
31.6 rad/s -90° (特征频率)
100 rad/s -158.2°
phaseplot([1000],[1,40,1000])
叠加相频图
phaseplot([sys1],[sys2]...[sysN]) 叠加相频图是将多个系统的相位响应曲线绘制在同一坐标系中的可视化工具，主要用于分析不同系统在频域中的相位特性.
控制系统稳定性分析（航空航天）:
LQR控制器在关键频段(1-10 rad/s)相位滞后减少15°，显著提升飞行稳定性
phaseplot([1,0.2s+1], # 基础PID控制器
[1,0.05*s^2+0.1s+1]) # 改进型LQR控制器
滤波器设计对比（音频处理）
高保真音频选贝塞尔，抗干扰选切比雪夫
phaseplot([1,s^2+1.4s+1], # 巴特沃斯滤波器
[1,s^2+s+1], # 切比雪夫滤波器
[1,s^2+0.5s+1]) # 贝塞尔滤波器
通信系统多径效应分析（5G网络）:
识别2.4GHz频段多径引起的相位突变，优化MIMO天线配置
phaseplot([1,0.01s+1], # 视距信道
[1,0.1*s^2+0.3s+1], # 城市多径信道
[1,0.5*s^2+0.1s+1]) # 室内密集多径
机械系统谐振分析（汽车悬架）:
低阻尼系统在谐振频率(≈2Hz)相位突变180°
高阻尼系统相位变化平缓，乘坐舒适性更好
phaseplot([1,0.3*s^2+0.8*s+1], # 标准阻尼
[1,0.3*s^2+1.2*s+1], # 高阻尼
[1,0.3*s^2+0.4s+1]) # 低阻尼
pi_digits 圆周率的小数精度
圆周率的小数精度
pi_digits(n)返回圆周率小数点后面的精度.
n -- 输入, 标量, 正整数.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def pi_digits_calculate(input_str):
"""
计算π的数值到指定有效数字位数（对标MATLAB的vpa(pi, n)功能）
参数:
input_str (str): 表示所需有效数字位数的字符串，范围1-500之间的整数
返回:
sp.Float: π的近似值，保留指定有效数字
str: 错误信息（当输入无效时）
示例:
>>> pi_digits_calculate("10")
3.1415926536
>>> pi_digits_calculate("1000")
'错误：输入值需为1-500之间的整数'
"""
try:
# 将输入字符串解析为SymPy表达式
expr = sp.sympify(input_str)
# 检查是否为元组（非法输入情况）
if isinstance(expr, tuple):
return f"输入错误: 不支持元组类型输入 - {input_str}"
# 验证输入是否为整数
if not expr.is_Integer:
return f"输入错误: 需输入整数 - {input_str}"
n = int(expr)
# 验证数值范围
if not (1 <= n <= 500):
return f"错误：输入值需为1-500之间的整数，当前输入：{n}"
# 计算并返回π值（保留n位有效数字）
return sp.N(sp.pi, n)
except sp.SympifyError:
return f"解析错误: 无法解析输入字符串 - {input_str}"
except Exception as e:
return f"运行时错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 正常案例
print("正常案例：")
print(pi_digits_calculate("10"))
# 3.141592654
print(pi_digits_calculate("5"))
# 3.1416
PiecewisePlot 分段函数图
分段函数图
PiecewisePlot(Piecewise(f(x)))处理分段函数，会区分不同区间或点集进行渲染（例如，用曲线表示连续部分，用点表示离散部分）。
x -- 输入, 符号变量.
阶梯电价模型
Piecewise(
(0.5 * x, And(x >= 0, x <= 100)), # 第一档电价
(50 + 0.8*(x-100), And(x > 100, x <= 200)), # 第二档
(130 + 1.2*(x-200), x > 200), # 第三档
(0, True) # 默认值（实际不会触发）
)
PiecewisePlot(Piecewise((0.5x, And(x>= 0, x<=100)),(50+0.8*(x-100), And(x >100,x<=200)),(130+1.2*(x-200),x> 200),(0,True)))
汽车制动距离模型
Piecewise(
(0.01*v**2, v <= 30), # 低速制动
(9 + 0.015*(v-30)**2, v > 30), # 高速制动
(0, True) # 默认值
)
PiecewisePlot(Piecewise((0.01*v^2, v<=30),(9+0.015*(v-30)^2, v>30),(0, True)))
医疗输液控制
Piecewise(
(20, And(t >= 0, t < 10)), # 初始快速滴注
(5, And(t >= 10, t < 60)), # 维持速率
(0, t >= 60), # 结束输液
(0, True) # 时间无效
)
PiecewisePlot(Piecewise((20, And(t>=0, t<10)),(5, And(t>=10, t<60)),(0, t>=60),(0, True)))
温度阈值控制系统
Piecewise(
(0, T < 10), # 低温关闭
(0.5*(T-10), And(T >= 10, T < 30)), # 线性调节
(10, T >= 30) # 高温限幅输出
)
PiecewisePlot(Piecewise((0, T<10),(0.5*(T-10), And(T>=10, T<30)),(10, T>= 30)))
绝对值函数
Piecewise(
(-x, x < 0), # 负区间
(x, x >= 0) # 非负区间
)
PiecewisePlot(Piecewise((-x, x<0),(x, x>=0)))
pinv 摩尔-彭罗斯伪逆
摩尔-彭罗斯伪逆
B = pinv(A) 返回矩阵 A 的摩尔-彭罗斯伪逆.
B是摩尔-彭罗斯伪逆,一种矩阵,可在不存在逆矩阵的情况下作为逆矩阵的部分替代.此矩阵常被用于求解没有唯一解或有许多解的线性方程组.
A — 输入矩阵,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def pseudo_inverse_matrix(input_str):
"""
计算给定矩阵的Moore-Penrose伪逆（对标MATLAB的pinv函数）。
参数:
input_str (str): 字符串形式的矩阵，支持以下格式：
- SymPy矩阵字符串，如"Matrix([[1, 0], [0, 0]])"
- 嵌套列表字符串，如"[[1, 0], [0, 0]]"
返回:
SymPy矩阵对象: 输入矩阵的伪逆（若输入有效）。
错误信息字符串: 若输入无效或计算失败。
示例:
>>> pseudo_inverse_matrix("[[1, 0], [0, 0]]")
Matrix([[1, 0], [0, 0]])
>>> pseudo_inverse_matrix("[[1, 2], [3, 4]]")
Matrix([
[ -2, 1],
[1.5, -0.5]])
>>> pseudo_inverse_matrix("5")
'输入错误: 5'
>>> pseudo_inverse_matrix("[[1, 2], [3, 4, 5]]")
'错误: ...' # 具体错误信息取决于异常
"""
try:
# 将输入字符串解析为SymPy表达式
expr = sp.sympify(input_str)
error_flag = False
result = None
# 检查并转换矩阵
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is not None:
# 计算伪逆
result = matrix.pinv()
else:
error_flag = True
return result if not error_flag else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 有效输入示例
print(pseudo_inverse_matrix("[[1, 0], [0, 0]]"))
# Matrix([[1, 0],
# [0, 0]])
print(pseudo_inverse_matrix("[[1, 2], [3, 4]]"))
# Matrix([[-2, 1],
# [3/2, -1/2]])
planerot 吉文斯平面旋转
吉文斯平面旋转
[G,y] = planerot(x),其中x是一个包含2个分量的列向量,返回一个2×2正交矩阵G,使得y = G*x有y(2) = 0.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def plane_givens_rotation(input_str):
"""
计算二维向量的吉文斯平面旋转矩阵（对标MATLAB的planerot函数）
参数:
input_str (str): 字符串形式的二维向量，支持格式：
- SymPy矩阵字符串: "Matrix([[3], [4]])"
- 嵌套列表字符串: "[[3], [4]]" 或 "[3, 4]"
返回:
tuple: (G, y) 元组，包含：
- G (Matrix): 2x2 吉文斯旋转矩阵
- y (Matrix): 旋转后向量，第二个元素为0
若输入有效则返回该元组
错误信息字符串: 输入无效时返回错误描述
示例:
>>> plane_givens_rotation("[3, 4]")
(Matrix([
[3/5, 4/5],
[-4/5, 3/5]]), Matrix([
[5],
[0]]))
>>> plane_givens_rotation("[[0], [0]]")
(Matrix([
[1, 0],
[0, 1]]), Matrix([
[0],
[0]]))
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def compute_givens(x):
"""计算吉文斯旋转的核心逻辑"""
a, b = x[0], x[1]
# 处理零向量特殊情况
if a == 0 and b == 0:
return sp.eye(2), sp.Integer(0)
# 计算旋转参数
r = sp.sqrt(a ** 2 + b ** 2).evalf()
c = a / r
s = b / r
# 构建吉文斯矩阵
G = sp.Matrix([
[c, s],
[-s, c]
]).evalf()
return G, r
# 矩阵有效性检查
matrix_expr = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix_expr is not None:
# 验证是否为二维向量
if matrix_expr.rows * matrix_expr.cols == 2:
x = matrix_expr.reshape(2, 1) # 强制转换为列向量
G, r = compute_givens(x)
y = G @ x # 矩阵乘法计算旋转后向量
result = (G, y)
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__":
# 正常案例测试
test_cases = [
"[3, 4]",
# (Matrix([[ 0.6, 0.8],
# [-0.8, 0.6]]),
# Matrix([[5.0],
# [ 0]]))
"[[1], [0]]",
# (Matrix([[1.0, 0],
# [ 0, 1.0]]),
# Matrix([[1.0],
# [ 0]]))
"[5, 12]",
# (Matrix([[ 0.384615384615385, 0.923076923076923],
# [-0.923076923076923, 0.384615384615385]]),
# Matrix([[13.0],
# [ 0]]))
]
for case in test_cases:
print(f"输入: {case}")
print("结果:", plane_givens_rotation(case))
print("-" * 40)
Plot 一元函数图
一元函数图 Plot(f(x))
描绘的是自变量 x 和因变量 y = f(x) 之间的静态关系. 直观展示函数在整个定义域或某个区间上的整体形态、变化趋势和关键特征。
x -- 变量
核心概念与基础性质:
多项式函数 (展示根、极值点、拐点)： 清晰展示三次多项式的行为：一个局部极大值、一个局部极小值、一个拐点、三个实根（或一个实根两个复根，取决于具体值）。观察导数为零处（极值点）和二阶导数为零处（拐点）。
Plot(x^3-6x^2+4x+12,x=[-2,8])
三角函数组合 (展示叠加、拍频): 展示不同频率三角函数的叠加效果。可以观察到低频包络线（由频率差决定）和高频振荡（由频率和决定），即“拍频”现象。
Plot(sin(x)+cos(3x),x=[-2@pi,2@pi])
指数衰减 (展示渐近行为): 经典的自然衰减模型（如放射性衰变、RC电路放电）。清晰展示函数从1开始单调递减，并无限趋近于x轴（y=0）但永不触及（渐近线）
Plot(exp(-x),x=[0,5])
对数函数 (展示增长缓慢、定义域限制): 展示对数函数的特性：增长极其缓慢（对比指数函数），是偶函数（因为 x²+1），始终在x轴上方（值域 [0, ∞)），在x=0处取得最小值0
Plot(ln(x^2+1),x=[-5,5])
应用场景与模型:
阻尼振荡 (物理 - 弹簧振子、RLC电路)： 模拟受阻尼的简谐振动。指数部分 e^{-0.2x} 代表振幅的衰减（阻尼），sin(x) 代表振荡。图像显示振荡幅度随时间逐渐减小。
Plot(exp(-0.2x)*sin(x),x=[0,20])
正态分布 (概率统计): 展示著名的钟形曲线。对称于y轴，在 x=0 处取得最大值 1/√(2π) ≈ 0.4，在 |x| > 2 后迅速衰减接近0。这是统计学中最核心的分布之一。
Plot(exp((-x^2)/2)/sqrt(2@pi),x=[-3,3])
需求曲线 (经济学): 展示典型的需求定律：价格(x)上涨，需求量(f(x))下降。函数单调递减，有下界（这里假设最小需求为0.5）。曲线形状反映需求弹性
Plot(1/(x+1)+0.5,x=[0,10])
Sinc 函数 (信号处理): 在信号处理（特别是采样定理）中极其重要。函数关于y轴对称，在 x=0 处值为1，在 x = ±kπ (k=1,2,3...) 处值为0。振幅随着 |x| 增大而衰减。主瓣和旁瓣清晰可见。
Plot(sin(x)/x,x=[-20,20])
有趣性质与特殊函数:
绝对值函数与振荡 (展示尖点与振荡): 函数在 x=0 处有尖点（不可导）。|x| 使得振幅随 |x| 增大而增大，sin(x) 提供振荡。图像像一条振幅不断增大的波浪线，且在原点“尖锐”地穿过x轴。
Plot(abs(x)*sin(x),x=[-10,10])
指数爆炸与振荡: 与阻尼振荡相反。指数部分 e^{0.1x} 代表振幅的指数级增长，sin(x) 提供振荡。图像显示振荡幅度随时间迅速增大。可类比某些不稳定系统或正反馈。
Plot(exp(0.1x)*sin(x),x=[0,30])
函数极限的经典例子: 当 x 趋近于0时，函数值在 [-1, 1] 之间无限次、无限快速地振荡。图像在0附近非常密集，无法趋近于任何一个特定的值，直观展示极限不存在。这是分析学中的一个经典例子。
Plot(sin(1/x),x=[-1,1])
多曲线函数图 Plot(f1(x), f2(x),...)
多个相关函数置于同一坐标系下，极大地增强了单图的信息量和分析能力，直观地揭示了仅靠单个函数图或分开的多个单函数图难以捕捉的对比关系和差异特征。
比较奇偶性/对称性: 直观展示 sin(x) 是奇函数 (关于原点对称)，sin(-x) 与之重合或镜像
Plot(sin(x),sin(-x),x=[-2@pi,2@pi])
展示参数影响 (阻尼振动): 清晰对比不同阻尼系数 (0.1, 0.5, 2.0) 对简谐振子振幅衰减速度的影响
Plot(exp(-0.1x)*sin(x),exp(-0.5x)*sin(x),exp(-2x)*sin(x),x=[0,20])
比较不同阶多项式的增长：清晰展示线性、二次、三次函数在 x > 1 区域增长速度的显著差异 (x^3 > x^2 > x)，直观理解函数阶数对增长趋势的影响。
Plot(x,x^2,x^3,x=[0,5])
展示基础函数的振荡差异：比较不同频率（ω=1 vs ω=2）正弦波的振荡周期和密度。直观显示频率加倍如何使波形在相同 x 范围内振荡次数翻倍。
Plot(sin(x),sin(2x), x=[0,4@pi])
比较衰减速率：对比不同衰减常数（λ=1 vs λ=0.5）的指数衰减。清晰显示衰减常数越小（0.5 < 1），函数值衰减得越慢。
Plot(exp(-x), exp(-0.5x),x=[0,10])
理解函数界限： 直观验证正弦函数的取值范围 [-1, 1] 以及它在哪些 x 点达到边界。
Plot(sin(x),1,-1,x=[0,2@pi])
比较不同优化器的收敛 (简化)：用简单函数模拟优化过程中损失随迭代次数 x 的下降。1/x 代表较快下降（如 Momentum），1/√x 代表较慢下降（如基础 SGD）。对比收敛速度差异。
Plot(1/x,1/(x^0.5), x=[1,100])
比较概率分布密度形状：对比两种重要概率密度函数：正态分布（钟形、轻尾）和贝塔分布（尖峰、重尾）。突出它们在峰值陡峭度和尾部衰减速度上的显著差异。
Plot([normpdf(x,2,1), betapdf(x,2,1)], x=[-5, 5])
时域函数图 Plot(f(x,T))
在绘制函数 f(x+T) 时，我们不是用当前时间 x 的值，而是用 T 时间单位之后（或之前）的时间点 (x + T) 的值。 观察和分析一个函数（信号）在时间轴上发生平移后的效果。
x -- 变量
T -- 动画参数
基础波形动画:
标准正弦波平移: 完整正弦波沿x轴向左平滑移动
Plot(sin(x + T))
压缩正弦波: 频率加倍的正弦波快速向左移动
Plot(sin(2x + T))
衰减波包: 振幅指数衰减的波包向右传播
Plot(e^(-0.1x) * sin(x + T))
调制与复合波形:
调幅波(AM): 载波频率5Hz，振幅以0.25Hz低频振荡
Plot((1 + 0.5*sin(0.5T)) * sin(5x + T))
调频波(FM): 瞬时频率随T周期性变化
Plot(sin(3x + 2*sin(0.3x+T)*x))
物理现象模拟:
驻波: 波节位置固定，波腹振幅周期性变化
Plot(sin(x-T) + sin(x+T))
行波反射: 入射波与衰减反射波叠加
Plot(sin(x-T) + 0.7*sin(-x-T))
特殊波形:
混沌振荡: 双频率耦合产生准周期运动
Plot(sin(x+T) + 0.5*sin(πx + 2T))
多曲线时域动画图 Plot(f1(x,T), f2(x,T), ...)
将多条曲线、时间演化和动画结合在一起，用于展示随时间变化的多组数据或函数关系。 揭示动态系统中多个相关量随时间和另一个变量共同演化的复杂关系。
x -- 变量
T -- 动画参数
基础波形组:
两条正弦波交替呈现"峰-谷"互补
Plot(
sin(x + T), # 基准波
sin(x + T + π/2) # 超前90°的波
)
频率阶梯: 不同频率波形同步平移，展示谐波关系
Plot(
sin(x + T), # 基波
sin(2x + T), # 二倍频波
sin(3x + T) # 三倍频波
)
调制对比: 振幅变化 vs 相位抖动的直观对比
Plot(
(1 + 0.5*sin(T)) * sin(x), # 调幅波
sin(x + 2*sin(T)) # 调相波
)
物理现象组:
行波与驻波: 观察行波如何叠加形成固定波节的驻波
Plot(
sin(x - T), # 向右传播的行波
sin(x - T) + sin(x + T) # 形成的驻波
)
多普勒频移模拟: 运动波源的波形压缩/拉伸效应
Plot(
sin(10x), # 静止波源
sin(10*(x - 0.5*sin(T))) # 运动波源的波
)
偏振合成: 线偏振→圆偏振的动态转变
Plot(
sin(T)*sin(x), # X方向分量
cos(T)*cos(x), # Y方向分量
sin(T)*sin(x) + cos(T)*cos(x) # 合成光波
)
Plot3D 三维函数图
单曲面三维函数图
Plot3D(f(x,y)) 直观呈现两个自变量如何共同影响因变量（非简单叠加）, 将抽象数学描述（如偏微分方程解）转化为可理解的几何形态.
波动干涉模型（物理学）: 模拟两个点波源（位置 (2,0) 和 (-2,0)）的干涉波纹。三维曲面清晰展示波峰叠加（相长干涉）和波谷抵消（相消干涉）的区域，直观解释干涉条纹的形成原理。
Plot3D(sin(@pi*sqrt((x-2)^2+y^2))+sin(@pi*sqrt((x+2)^2+y^2)),x=[-5,5],y=[-5,5])
优化问题的目标函数曲面（数学/运筹学）: 绘制 Beale函数（经典优化测试函数）。曲面呈现多个局部极小值和狭窄的全局最优谷底，突出优化算法（如梯度下降）易陷入局部最优的挑战。
Plot3D((1.5-x+x*y)^2+(2.25-x+x*y^2)^2+(2.625-x+x*y^3)^2,x=[-4,4],y=[-4,4])
热传导稳态分布（工程学）: 模拟矩形金属板在 x=0 边恒温加热时的稳态温度分布。展示温度沿 x 方向正弦振荡、沿 y 方向指数衰减的物理规律，高温区和低温区一目了然。
Plot3D(50*sin(@pi*x/2)*exp(-@pi*y/2),x=[0,2],y=[0,3])
经济学中的效用曲面（微观经济学）: Cobb-Douglas效用函数，描述消费者对两种商品 (x,y) 的满意度。曲面呈现边际效用递减（随 x 或 y 增加，效用增长变缓）和无差异曲线投影（等高线）。
Plot3D(x^0.4*y^0.6,x=[0,5],y=[0,5])
量子力学概率密度（物理学）: 氢原子 2pₓ 轨道的概率密度分布。环形山峰反映电子出现的概率在原点处为零，在特定半径处达到峰值，直观展示量子态的非均匀特性。
Plot3D((x^2+y^2)*exp(-(x^2+y^2)),x=[-3,3],y=[-3,3])
地形高程模型（地理学）: Matlab的Peaks函数，模拟多峰山脉地形。同时显示主峰、次峰、山谷和鞍部，用于测试地形渲染算法或路径规划。
Plot3D(3*(1-x)^2*exp(-x^2-(y+1)^2)-10*(x/5-x^3-y^5)*exp(-x^2-y^2)-(1/3)*exp(-(x+1)^2-y^2),x=[-3,3], y=[-3,3])
流体速度场势函数（流体力学）: 不可压缩流体的速度势。曲面等高线对应流线，梯度方向即流速方向，直观呈现涡旋和流动模式。
Plot3D((cos(2*@pi*x)/2+sin(2*@pi*y)/2),x=[0,1],y=[0,1])
多曲面三维函数图
Plot3D(f1(x,y),f2(x,y)...) 同时可视化多个三维函数（两个自变量）能揭示单曲面无法呈现的空间关系、对称性、包络效应和系统行为，在科学和工程中具有关键作用.
电磁场正交分量: 展示电磁波传播中电场与磁场的空间正交性和能量传递相位差
Plot3D(sin(x)*cos(y),cos(x)*sin(y),x=[-@pi,@pi],y=[-@pi,@pi])
量子叠加态概率密度: 对比量子粒子在不同能级的概率分布差异，交叉区域显示量子隧穿效应
Plot3D(exp(-(x-1)^2-y^2),exp(-(x+1)^2-y^2),x=[-3,3],y=[-3,3])
流体速度剖面: 对比管道流体中层流与湍流的速度分布差异，波纹曲面显示涡旋结构
Plot3D(1-(x^2+y^2),1-sqrt(x^2+y^2)+0.2*sin(10*sqrt(x^2 + y^2)),x=[-1,1],y=[-1,1])
地形与水位线: 显示洪水淹没区域（曲面交线以下），用于灾害模拟
Plot3D(3*exp(-x^2-(y-1)^2),0.5,x=[-3,3],y=[-3,3])
神经网络损失曲面对比: 对比不同优化问题的复杂性，非凸曲面的局部极小值陷阱清晰可见
Plot3D(x^2+y^2,x^2-y^2+0.5*sin(3@pi*x)*cos(2@pi*y),x=[-2,2],y=[-2,2])
声波干涉模式: 叠加曲面显示相长/相消干涉条纹，峰值重合区为最大振幅点
Plot3D(cos(@pi*sqrt((x-2)^2+y^2)),cos(@pi*sqrt((x+2)^2+y^2)),x=[-5,5],y=[-5,5])
经济生产可能性边界: 显示不同生产技术的产出边界，交线为最优资源配置路径
Plot3D(10-x^2-y^2, 8-0.5*(x-1)^2-0.5*(y-1)^2,x=[0,3],y=[0,3])
pochhammer 阶乘幂
阶乘幂
pochhammer(x,n)返回 pochhammer符号(x)n
x —— 输入,数字,向量,矩阵,符号数字,符号变量,符号向量,符号矩阵,符号表达式
n —— 输入,数字,向量,矩阵,符号数字,符号变量,符号向量,符号矩阵,符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import poch
def pochhammer_rising_factorial(input_str):
"""
计算Pochhammer升阶乘符号（对标MATLAB的pochhammer函数）
参数:
input_str (str): 字符串形式的输入，应为包含两个元素的元组：
- 格式示例: "(x, n)" 或 "(Matrix([[1,2]]), 3)"
- 元素可以是标量或矩阵
返回:
Matrix/Expr: 计算结果，有以下情况：
- 标量输入返回符号表达式
- 矩阵输入返回同型矩阵
- 错误时返回描述字符串
示例:
>>> pochhammer_rising_factorial("(5, 3)")
210
>>> pochhammer_rising_factorial("(Matrix([[1,2],[3,4]]), 2)")
Matrix([
[ 2, 6],
[12, 20]])
>>> pochhammer_rising_factorial("(a, 3)")
a*(a + 1)*(a + 2)
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 验证输入为二元组
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = poch(*params) # x和n为数值类型
elif any(e.free_symbols for e in expr):
result = sp.rf(*expr).rewrite(sp.gamma)
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__":
# 正常输入测试
tests = [
("(5, 3)", "标量计算"),
# 210.0
("(a, 3)", "符号计算"),
# a*(a + 1)*(a + 2)
]
for case, desc in tests:
print(f"测试案例: {desc}")
print("输入:", case)
print("输出:", pochhammer_rising_factorial(case))
print("=" * 60)
poisscdf 泊松累积分布函数
泊松累积分布函数
y = poisscdf(x,lambda) 使用 lambda 中的速率参数计算在 x 中每个值处的泊松累积分布函数值.
x 和 lambda 可以是大小都相同的标量,向量或矩阵. 如果只有一个参量是标量, poisscdf会将其扩展为其维数与另一个参量相同的常量数组.
x — 计算泊松 cdf 所基于的值,标量值,标量值组成的数组
lambda — 速率参数,正值,正值组成的数组
y — 泊松 cdf 值,标量值,标量值组成的数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Poisson, cdf
from scipy.stats import poisson
def poisson_distribution_cdf(input_str):
"""
计算泊松累积分布函数（对标MATLAB的poisscdf函数）
参数:
input_str (str): 字符串形式的输入，应为包含两个元素的元组：
- 格式示例: "(x, lambda)" 或 "(Matrix([[0,1,2]]), 2)"
- x: 事件数（支持标量/矩阵）
- lambda: 泊松分布参数（正数，支持标量/矩阵）
返回:
Matrix/Expr: 计算结果，有以下情况：
- 标量输入返回数值
- 矩阵输入返回同型矩阵
- 错误时返回描述字符串
示例:
>>> poisson_distribution_cdf("(2, 1)")
0.919698602928606
>>> poisson_distribution_cdf("(Matrix([[0,1,2]]), 2)")
Matrix([[0.135335283236613, 0.406005849709838, 0.676676416183063]])
>>> poisson_distribution_cdf("(-1, 2)")
'错误: Lambda参数必须为正数'
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 验证输入格式为二元组
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = poisson.cdf(*params)
elif any(e.free_symbols for e in expr):
x, lam = expr
# 参数验证
# 创建泊松随机变量
X = Poisson("poisson", lam)
result = cdf(X)(x)
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f'错误:{str(e)}'
# 测试用例
if __name__ == "__main__":
# 正常案例测试
tests = [
("(2, x)", "符号计算"),
# (-(-2*(x**2/2 + x + 1)*exp(-x) + 2)*exp(x)/2 + exp(x))*exp(-x)
("(2.5, 3)", "标量计算 "),
# 0.42319008112684364
]
for case, desc in tests:
print(f"测试案例: {desc}")
print("输入:", case)
print("输出:", poisson_distribution_cdf(case))
print("=" * 60)
poissfit 泊松参数估计
泊松参数估计
[lambdahat,lambdaci] = poissfit(data) 还在 lambdaci 中给出 95% 的置信区间.
[lambdahat,lambdaci] = poissfit(data,alpha) 给出 100(1 - alpha)% 的置信区间。例如，alpha = 0.001 得出 99.9% 的置信区间.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import numpy as np
from scipy import stats
def poisson_parameter_estimation(input_str):
"""
泊松分布参数估计（对标MATLAB的poissfit函数）
参数:
input_str (str): 字符串形式的输入，支持两种格式：
1. 纯数据列表: "[1,2,3,4,5]"
2. 带置信水平的元组: "([0,1,2,3], 0.1)"
返回:
tuple: (lambda_hat, (ci_lower, ci_upper)) 包含：
- lambda_hat (float): λ的估计值
- ci_lower (float): 置信区间下限
- ci_upper (float): 置信区间上限
若输入有效则返回该元组
错误信息字符串: 输入无效时返回错误描述
示例:
>>> poisson_parameter_estimation("[1,1,2,3,5]")
(2.4, (1.3859038246098037, 3.702993444456167))
>>> poisson_parameter_estimation("([2,2,2], 0.01)")
(2.0, (1.0793132286998841, 3.687433386962528))
"""
try:
# 解析输入字符串
parsed = ast.literal_eval(input_str)
# 初始化参数
data = None
alpha = 0.05 # 默认置信水平
# 解析输入格式
if isinstance(parsed, tuple):
if len(parsed) != 2:
raise ValueError("元组输入需要包含数据和alpha两个元素")
data, alpha = parsed
elif isinstance(parsed, list):
data = parsed
else:
raise TypeError("输入格式不正确")
# 转换数据类型
data_array = np.array(data, dtype=float)
# 参数验证
if data_array.size == 0:
raise ValueError("输入数据不能为空")
if np.any(data_array < 0):
raise ValueError("泊松数据不能包含负值")
if not (0 < alpha < 1):
raise ValueError("置信水平alpha必须在0到1之间")
# 计算参数估计
n = len(data_array)
total_counts = np.sum(data_array)
lambda_hat = total_counts / n
# 计算置信区间 (基于卡方分布)
chi2_lower = stats.chi2.ppf(alpha / 2, 2 * total_counts)
chi2_upper = stats.chi2.ppf(1 - alpha / 2, 2 * (total_counts + 1))
ci_lower = chi2_lower / (2 * n)
ci_upper = chi2_upper / (2 * n)
return (lambda_hat, (ci_lower, ci_upper))
except Exception as e:
return f"错误: {str(e)}"
# 测试用例
if __name__ == "__main__":
# 正常案例测试
test_cases = [
("[1,1,2,3,5]", "标准输入"),
# (2.4, (1.2401150217444434, 4.192317009635392))
("([0,0,0,0], 0.1)", "全零数据"),
# (0.0, (nan, 0.7489330683884974))
("([5], 0.01)", "单数据点"),
# (5.0, (1.0779282406523194, 14.149759411023012))
("([2,2,2,2], 0.2)", "均匀数据")
# (2.0, (1.1640295442245008, 3.248677885329651))
]
for case, desc in test_cases:
print(f"测试案例: {desc}")
print("输入:", case)
print("输出:", poisson_parameter_estimation(case))
print("=" * 60)
poissinv 泊松逆累积分布函数
泊松逆累积分布函数
X = poissinv(P,lambda)返回最小值X,使得在X处计算的泊松cdf等于或超过P,使用lambda中的平均参数. P和lambda可以是大小相同的向量或矩阵. 标量输入被扩展为与其他输入具有相同维度的常数数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy import stats
from sympy.stats import Poisson, quantile
def poisson_inverse_cumulative(input_str):
"""
泊松逆累积分布函数（对标MATLAB的poissinv函数）
参数:
input_str (str): 字符串形式的输入，应为包含两个元素的元组：
- 格式示例: "(0.9, 5)" 或 "(Matrix([[0.9, 0.95]]), 5)"
- p: 累积概率值（支持标量/矩阵），范围[0,1)
- lambda: 泊松分布参数（正数，支持标量/矩阵）
返回:
Matrix/ndarray: 计算结果，有以下情况：
- 标量输入返回整数
- 矩阵输入返回同型矩阵
- 错误时返回描述字符串
示例:
>>> poisson_inverse_cumulative("(0.95, 5)")
9
>>> poisson_inverse_cumulative("([[0.5, 0.9], [0.95, 0.99]], 10)")
[[6 13]
[14 18]]
>>> poisson_inverse_cumulative("(1.1, 5)")
'错误: 概率值必须在[0, 1)区间'
"""
try:
# 安全解析输入字符串
expr = sp.sympify(input_str)
error = False
result = None
def eval_poisson_inverse_cdf(p_val, lambda_val):
return stats.poisson.ppf(float(p_val), float(lambda_val)).astype(int)
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 = stats.poisson.ppf(*params)
elif any(e.free_symbols for e in expr):
p, lambda_sym = expr
X = Poisson("X", lambda_sym)
result = quantile(X)(p)
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 测试用例
if __name__ == "__main__":
# 正常案例测试
tests = [
("(0.8, 2)", "标量计算"),
# 3.0
("(0, x)", "边界值测试"),
# Piecewise((_x, (exp(x) - sinh(x))*exp(-x) >= 0))
]
for case, desc in tests:
print(f"测试案例: {desc}")
print("输入:", case)
print("输出:", poisson_inverse_cumulative(case))
print("=" * 60)
poisspdf 泊松概率密度函数
泊松概率质量函数
y = poisspdf(x,lambda) 使用 lambda 中的速率参数计算 x 中每个值处的泊松概率密度函数值.
x 和 lambda 可以是大小都相同的标量,向量或矩阵. 如果只有一个参量是标量, poisspdf 会将其扩展为其维数与另一个参量相同的常量数组.
x — 计算泊松 pdf 时所基于的值,标量值,标量值组成的数组
lambda — 速率参数,正值,正值组成的数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Poisson, density
from scipy.stats import poisson as sci_poisson
def poisson_distribution_pdf(input_str):
"""
计算泊松概率质量函数（对标MATLAB的poisspdf函数）
参数:
input_str (str): 字符串形式的输入，应为包含两个元素的元组：
- 格式示例: "(x, lambda)" 或 "(Matrix([[0,1,2]]), 5)"
- x: 事件数（非负整数，支持标量/矩阵）
- lambda: 泊松参数（正数，支持标量/矩阵）
返回:
Matrix/Expr: 计算结果，有以下情况：
- 标量输入返回符号表达式
- 矩阵输入返回同型矩阵
- 错误时返回描述字符串
示例:
>>> poisson_distribution_pdf("(2, 1)")
0.183939720585721
>>> poisson_distribution_pdf("(Matrix([[0,1,2]]), 2)")
Matrix([[0.135335283236613, 0.270670566473225, 0.270670566473225]])
>>> poisson_distribution_pdf("(-1, 2)")
'错误: 事件数必须为非负整数'
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = sci_poisson.pmf(*params)
elif any(e.free_symbols for e in expr):
x, lambda_val = expr
X = Poisson("poisson", lambda_val)
result = density(X)(x)
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f'错误: {str(e)}'
# 测试用例
if __name__ == "__main__":
# 正常输入测试
tests = [
("(2, 1)", "标量计算"),
# 0.18393972058572114
("(k, 3)", "符号计算"),
# 3**k*exp(-3)/factorial(k)
]
for case, desc in tests:
print(f"测试案例: {desc}")
print("输入:", case)
print("输出:", poisson_distribution_pdf(case))
print("=" * 60)
poissrnd 来自泊松分布的随机数
来自泊松分布的随机数
r = poissrnd(lambda) 从速率参数 lambda 所指定的泊松分布中生成随机数.
lambda 可以是标量,向量或矩阵.
r = poissrnd(lambda,sz1,sz2) 从具有标量速率参数 lambda 的泊松分布中生成随机数数组,其中 sz1,sz2 表示矩阵维度的大小.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.stats as stats
import numpy as np
def eval_poisson_rnd(lambda_param, size=None):
"""
生成泊松分布的随机数。
参数:
lambda_param: 泊松分布的λ参数。
size: 输出的维度，默认为None（生成单个随机数）。
返回:
根据size参数返回单个随机数、列表或嵌套列表。
"""
lambda_val = float(lambda_param)
if lambda_val < 0:
raise ValueError("泊松分布的λ参数必须非负。")
poisson_dist = stats.poisson(lambda_val)
random_numbers = poisson_dist.rvs(size=size)
# 将结果转换为Python原生类型
if size is None:
return int(random_numbers)
elif isinstance(random_numbers, np.ndarray):
return random_numbers.tolist()
elif isinstance(random_numbers, (int, np.integer)):
return int(random_numbers)
else:
return list(random_numbers)
def poisson_random_numbers(input_str):
"""
根据输入字符串生成泊松分布的随机数或矩阵。
参数:
input_str: 输入字符串，可以是λ参数、维度或矩阵。
返回:
生成的随机数、列表或SymPy矩阵。若输入错误则返回错误信息。
示例:
>>> poisson_random_numbers("3") # 单个随机数
>>> poisson_random_numbers("3,5") # 包含5个随机数的列表
>>> poisson_random_numbers("3,2,3") # 2x3的随机矩阵
>>> poisson_random_numbers("[[2,3],[4,5]]") # 2x2矩阵，每个元素对应λ生成的随机数
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
# 处理元组输入，如(λ, m, n)或(λ, size)
if len(expr) == 3 and expr[1].is_number and expr[2].is_number:
lambda_val = expr[0]
rows = int(expr[1])
cols = int(expr[2])
rand_data = eval_poisson_rnd(lambda_val, size=(rows, cols))
result = sp.Matrix(rand_data)
elif len(expr) == 2 and expr[1].is_number:
lambda_val = expr[0]
size = int(expr[1])
rand_data = eval_poisson_rnd(lambda_val, size=size)
result = rand_data # 直接返回列表
else:
error = True
elif isinstance(expr, list):
matrix = sp.Matrix(expr)
# 输入是矩阵，生成对应形状的随机矩阵
rows, cols = matrix.shape
rand_matrix = []
for i in range(rows):
row = []
for j in range(cols):
lambda_ij = matrix[i, j]
row.append(eval_poisson_rnd(lambda_ij, size=None))
rand_matrix.append(row)
result = sp.Matrix(rand_matrix)
elif expr.is_number:
# 单个λ参数生成单个随机数
result = eval_poisson_rnd(expr, size=None)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：生成λ=3的单个随机数
print("示例1:", poisson_random_numbers("3"))
# 0
# 示例2：生成λ=3的5个随机数
print("示例2:", poisson_random_numbers("3,5"))
# [4, 2, 1, 2, 2]
# 示例3：生成λ=3的2x3随机矩阵
print("示例3:", poisson_random_numbers("3,2,3"))
# Matrix([[4, 2, 4],
# [4, 3, 4]])
# 示例4：生成与输入矩阵对应的随机矩阵
print("示例4:", poisson_random_numbers("[[2,3],[4,5]]"))
# Matrix([[3, 4],
# [3, 4]])
poisstat 泊松均值和方差
泊松均值和方差
M = poisstat(lambda) 使用 lambda 中的均值参数返回泊松分布的均值. M的大小就是lambda的大小.
[M,V] = poisstat(lambda) 还返回泊松分布的方差 V.
对于带参数λ的泊松分布，均值和方差都等于λ.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Poisson, E, variance
def evaluation_poisson_stat(x_val):
"""
计算泊松分布的期望和方差。
当x_val为符号时，返回符号表达式。
参数:
x_val: 泊松分布的λ参数，可以是数值或符号。
返回:
(期望, 方差) 的元组，类型为SymPy表达式或数值。
"""
if x_val.is_number and float(x_val) < 0:
raise ValueError("泊松分布的λ参数必须非负。")
# 定义泊松随机变量
X = Poisson('X', x_val)
mean_val = E(X)
var_val = variance(X)
return mean_val, var_val
def poisson_mean_variance(input_str):
"""
根据输入计算泊松分布的均值和方差。
支持数值、符号、矩阵输入，对标MATLAB的poisstat函数。
参数:
input_str: 输入字符串，可以是数值、符号表达式或矩阵。
返回:
均值矩阵/值和方差矩阵/值的元组，若输入错误则返回错误信息。
示例:
>>> poisson_mean_variance("3") # 返回 (3, 3)
>>> poisson_mean_variance("lambda") # 返回 (lambda, lambda)
>>> poisson_mean_variance("[[2,3],[4,5]]") # 返回 (2x2矩阵, 2x2矩阵)
"""
try:
expr = sp.sympify(input_str)
error = False
M, V = None, None
# 处理元组输入（不支持）
if isinstance(expr, tuple):
error = True
raise ValueError("不支持元组输入，请输入数值、符号或矩阵。")
# 处理数值/符号输入
elif expr.is_number or expr.free_symbols:
# 检查数值是否为负
if expr.is_number and float(expr) < 0:
raise ValueError("λ参数必须非负。")
M, V = evaluation_poisson_stat(expr)
else:
error = True
return (M, V) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：数值输入
print("示例1:", poisson_mean_variance("3"))
# (3, 3)
# 示例2：符号输入
lambda_sym = sp.symbols('x')
print("示例2:", poisson_mean_variance(lambda_sym))
# (x, x)
PolarPlot 极坐标图
极坐标图
PolarPlot(r=f(theta)) 在极坐标系中绘制的图形，其中每个点由两个参数定义：
半径 r: 表示点到原点的距离。
角度 thetaθ: 表示点相对于正x轴的角度（弧度制，通常从0到2π 或更广范围）。
这种函数图适合表示周期性、旋转对称或圆形分布的图案，常见于物理、工程、艺术和自然界（如花瓣、声波、天线辐射模式）。
圆: 表示理想化的圆形路径或物体，如车轮轨迹、圆形池塘的边界。
在工程中，用于简化模型（如轴承滚珠的路径）。
PolarPlot(2)
心形线: 常见于声学设备，如心形指向性麦克风的拾音模式（对前方声音敏感，后方抑制）
也用于天线设计，优化信号方向性。
PolarPlot(2(1+cos(@theta)))
玫瑰曲线: 模拟周期性自然现象
如花瓣排列（植物学）、声波干涉图样（物理实验），或艺术设计（如教堂彩窗图案）。
PolarPlot(3cos(4@theta))
阿基米德螺线: 等间距螺旋线
机械工程：卷起的绳子、钟表发条。
天文学：简化模型表示星系旋臂（如银河系）。
日常用品：DVD 光轨、风扇叶片的设计基础。
PolarPlot(0.5@theta)
对数螺线: 螺距逐渐增大的螺旋线
生物学：鹦鹉螺壳、向日葵种子排列（黄金角优化）。
气象学：飓风云图的结构。
艺术：用于建筑和雕塑（如古希腊柱头设计）。
PolarPlot(exp(0.1@theta))
多曲线极坐标图
PolarPlot(f1(theta),f2(theta)) 通过叠加多条r(θ)函数，直观对比不同物理过程
天线辐射模式对比: 对比信号强度分布，优化无线网络覆盖。
PolarPlot(cos(@theta)^2,0.5*sin(@theta)^2)
机械振动分析: 检测旋转机械（如涡轮机）的异常振动模式。
PolarPlot(2+sin(5@theta),2+1.5*sin(5@theta))
双星系统轨道模拟: 模拟引力相互作用下的双星运动轨迹。
PolarPlot(3/(1+0.9*cos(@theta)),3/(1+0.9*cos(@theta+@pi)))
声波干涉分析: 可视化声波叠加产生的干涉瓣方向。
PolarPlot(2*abs(cos(3@theta)),2*abs(sin(3@theta)))
polezeroplot 零极点图
零极点图
polezeroplot(sys)是复平面（s平面或z平面）上标出系统传递函数零点（Zeros）和极点（Poles）位置的图形，用于分析系统特性（稳定性、频率响应等）
系统在频率 ω=1 rad/s 处增益最小（零点抑制）
陷波滤波器（Notch Filter），用于消除特定频率干扰（如 50Hz 工频噪声）。
polezeroplot(s^2+1,1)
所有极点在左半平面 → 稳定系统。
低通滤波器，极点聚集在实轴负半轴，衰减高频信号。
polezeroplot(1,s^4+4s^3+6s^2+5s+2)
带阻滤波器
零点在虚轴 ±2j 处 → 抑制 ω=2 rad/s 的信号（如消除机械振动噪声）。
polezeroplot(s^2+4,s^2+s+4)
不稳定系统
极点在右半平面（s=1）→ 系统不稳定（如正反馈电路）。
polezeroplot(s,s^2-1)
一阶低通滤波器（RC电路）
极点在左半平面 → 系统稳定。
polezeroplot([1],[1,1])
二阶带通滤波器（RLC电路）
允许特定频率通过（如收音机选频），零点在原点抑制直流分量，虚轴极点位置控制中心频率。
polezeroplot([1],[1,1,1])
全通滤波器（相位校正）
幅度响应平坦，但调整相位（如延迟均衡）。零极点关于虚轴对称 → 全通特性。
polezeroplot([1,-1],[1,1])
不稳定系统（振荡器）
极点位于右半平面 → 系统发散（如未阻尼振荡），实际需避免。
polezeroplot([1],[1,0,-1])
数字滤波器（音频均衡器）
增强特定音频频率（如人声），极点在单位圆内保证稳定，零点位置控制频响凹陷。
polezeroplot([1,-0.8],[1,-1.2,0.5])
陷波滤波器（消除工频干扰）
消除特定频率（如 50 Hz 电源噪声），零点强制该频率增益为 0，极点位置控制带宽。
polezeroplot([1,0,100],[1,2,101])
叠加零极点图
polezeroplot([sys1],[sys2]...[sysN])同一复平面上绘制多个系统传递函数的零极点.
滤波器组性能对比: 设计音频处理系统中的两个滤波器
系统1零点在 \pm j±j → 在 ω=1 rad/s 处强烈抑制信号
系统2极点在 -0.5 \pm j1.94−0.5±j1.94 → 在 ω≈2 rad/s 处有谐振峰
polezeroplot([s^2+1,s^2+4],[1,s^2+s+4])
控制器与被控对象交互: 电机控制系统
控制器零点吸引被控对象极点 → 可加速系统响应
无右半平面极点 → 闭环系统稳定
polezeroplot([1,s^2+2s+5],[0.5s^2+2s+3,s])
通信系统多径干扰抑制: 接收端两个抑制干扰的子系统
陷波器零点与极点接近虚轴 → 对 3 rad/s 干扰有深度抑制
抑制器极点远离虚轴 → 不影响陷波器性能
polezeroplot([s^2+9,s^2+10],[s+2,s^2+4s+20])
机械振动控制系统: 汽车悬架系统模型
主动控制器在 s=0s=0 处的极点 → 需检查积分器稳定性
主动系统零点 \pm 2j±2j 与被动系统谐振频率（7j）无冲突 → 可有效抑制 2 rad/s 路面振动
polezeroplot([s+1,s^2+2s+10],[5*(s^2+4),s^3+10s^2+50s])
pol2cart 将极坐标或柱坐标转换为笛卡尔坐标
将极坐标或柱坐标转换为笛卡尔坐标
[x,y] = pol2cart(theta,rho) 将极坐标数组theta和rho的对应元素变换为二维笛卡尔坐标或xy坐标.
[x,y,z] = pol2cart(theta,rho,z) 将柱坐标数组theta,rho和z的对应元素变换为三维笛卡尔坐标或xyz坐标.
theta — 角坐标,标量,向量,矩阵
rho — 径向坐标,标量,向量,矩阵
z — 仰角坐标,标量,向量,矩阵
x, y, z — 笛卡尔坐标,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def polar_to_cartesian(input_str):
"""
将极坐标或柱坐标转换为笛卡尔坐标，对标MATLAB的pol2cart函数。
参数:
input_str: 输入字符串，可以是以下形式：
- 二维坐标: "(theta, rho)"
- 三维坐标: "(theta, rho, z)"
其中theta和rho可以是数值、符号或矩阵
返回:
- 二维情况: (x, y) 的元组（数值/矩阵）
- 三维情况: (x, y, z) 的元组
若输入错误则返回错误信息。
示例:
>>> polar_to_cartesian("(0.785, 5)") # 二维数值输入
>>> polar_to_cartesian("([0, 1.57], [2, 3])") # 二维矩阵输入
>>> polar_to_cartesian("(0.785, 5, 10)") # 三维数值输入
"""
try:
expr = sp.sympify(input_str)
error = False
M, V = None, None
def pol2cart_sym(theta, rho):
# 使用SymPy的三角函数以支持符号计算
x = rho * sp.cos(theta)
y = rho * sp.sin(theta)
return x.evalf(), y.evalf()
def pol2cart_sci(theta, rho=None):
"""
将极坐标(或柱坐标)转换为笛卡尔坐标，对标MATLAB的pol2cart函数
参数:
theta: 角度(弧度)，可以是标量、向量或多维数组
rho: 径向距离(可选)，默认为单位圆
z: 高度(柱坐标的z坐标，可选)
返回:
x, y: 笛卡尔坐标x和y
(如果提供z，则返回x, y, z)
"""
# 处理不同输入情况
# 计算x和y坐标
x = rho * np.cos(theta)
y = rho * np.sin(theta)
return x, y
# 检查 expr 是否为元组类型
if isinstance(expr, tuple):
if all(e.is_number for e in expr):
params = tuple(float(e) for e in expr)
theta, rho = pol2cart_sci(*params[:2])
elif any(e.free_symbols for e in expr):
theta, rho = pol2cart_sym(*expr)
# 如果 expr 不是元组类型，说明输入不符合要求，将错误标志 error 设为 True
else:
error = True
if len(expr) >= 3:
return theta, rho, expr[2]
else:
return theta, rho
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：二维标量输入
print("示例1:", polar_to_cartesian("(0.785, 5)"))
# (3.5369413458359986, 3.53412590552683)
# 示例2：三维输入
print("示例3:",
polar_to_cartesian("(0.785, 5, 10)")
)
# (3.5369413458359986, 3.53412590552683, 10)
# 示例3：符号计算
theta_sym, rho_sym = sp.symbols('theta rho')
print("示例4:",
polar_to_cartesian(f"{theta_sym}, {rho_sym}")
)
# (rho*cos(theta), rho*sin(theta))
poly 具有指定根的多项式或特征多项式
具有指定根的多项式或特征多项式
p = poly(r)（其中r是向量）返回多项式的系数,其中多项式的根是r的元素.
p = poly(A)（其中A是n×n矩阵）返回矩阵det(λI – A) 的特征多项式的n+1个系数.
r — 多项式根,向量
A — 输入矩阵,矩阵
p — 多项式系数,行向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def polynomial_eig_matrix(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
def is_vector(matrix):
"""判断是否是向量（单行或单列）"""
return matrix.rows == 1 or matrix.cols == 1
# 尝试判断表达式是否为矩阵
if isinstance(expr, list):
matrix = sp.Matrix(expr)
if is_vector(matrix):
# 若矩阵为向量，将其转换为 NumPy 数组并展平
A = np.array(matrix.tolist(), dtype=float).ravel()
elif matrix.is_square:
# 若矩阵为方阵，将其转换为 NumPy 数组
A = np.array(matrix.tolist(), dtype=float)
else:
# 若矩阵既不是向量也不是方阵，设置错误标志
error = True
# 计算矩阵或向量的特征多项式的系数
result = np.poly(A)
# 设置 NumPy 输出格式，抑制科学计数法，精确到小数点后 8 位
np.set_printoptions(suppress=True, precision=8)
else:
# 若表达式不是矩阵，设置错误标志
error = True
return list(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：数值向量输入
print("示例1:", polynomial_eig_matrix("[1, 2]"))
# [1.0, -3.0, 2.0]
# 示例2：数值矩阵输入
print("示例2:", polynomial_eig_matrix("[[2,1],[1,2]]"))
# [1.0, -4.0, 3.0]
polyder 多项式微分
多项式微分
k = polyder(p) 返回 p 中的系数表示的多项式的导数.
p,q = polyder(a,b) 返回 p = 多项式a和b的乘积的导数, q= 多项式a和b的商的导数.
p是多项式系数,指定为向量.
a,b是多项式系数(作为单独参量),为行向量.
k是微分多项式系数,行向量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import sympy as sp
def polynomial_partial_derivative(input_str, mode='product'):
"""
对标MATLAB的polyder函数，计算多项式的导数。
参数:
input_str (str): 输入的多项式或元组的字符串表示，如 "[1,2,3]" 或 "([1,2], [3,4])"
mode (str): 如果是两个多项式，指定是 'product'（乘积导数，默认）或 'quotient'（商导数）。
返回:
单个多项式导数: 系数列表，如 [2, 2] 对应 2x + 2。
两个多项式乘积导数: 系数列表。
两个多项式商导数: 返回元组（分子系数，分母系数）。
错误时返回错误信息字符串。
"""
try:
expr = ast.literal_eval(input_str)
# 处理单个多项式
if isinstance(expr, list):
poly = sp.Poly(expr, sp.symbols('x'))
deriv = poly.diff()
return deriv.all_coeffs()
# 处理两个多项式的情况
elif isinstance(expr, tuple) and len(expr) == 2:
a, b = expr
if isinstance(a, list) and isinstance(b, list):
a_poly = sp.Poly(a, sp.symbols('x'))
b_poly = sp.Poly(b, sp.symbols('x'))
if mode == 'product':
# 计算乘积的导数
product_deriv = (a_poly * b_poly).diff()
return product_deriv.all_coeffs()
elif mode == 'quotient':
# 计算商导数的分子和分母
a_deriv = a_poly.diff()
b_deriv = b_poly.diff()
numerator = a_deriv * b_poly - a_poly * b_deriv
denominator = b_poly ** 2
return (numerator.all_coeffs(), denominator.all_coeffs())
else:
return f"错误: 未知模式 '{mode}'，应为 'product' 或 'quotient'"
else:
return f"输入错误: 元组元素必须为列表"
else:
return f"输入格式错误: 应为列表或包含两个列表的元组"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 单个多项式导数
p = "[3, 2, 1]"
print("单个多项式导数:", polynomial_partial_derivative(p))
# [6, 2]
# 示例2: 乘积导数（默认）
p_q = "([1, 2], [3, 4])"
print("\n乘积导数（默认）:", polynomial_partial_derivative(p_q))
# [6, 10]
# 示例3: 商导数（显式指定）
print("\n商导数:", polynomial_partial_derivative(p_q, mode='quotient'))
# ([-2], [9, 24, 16])
polyeig 多项式特征值问题
多项式特征值问题
[X,e] = polyeig(A0,A1,...,Ap) 返回大小为 n×n*p 的矩阵X,其列是特征向量,和p次多项式特征值问题的特征值e.
polyeig用于求解任意多项式特征值问题.比如在转子动力学的情况下,必须解决一个二次特征值问题.
eig用于求解广义特征值问题.
A0,A1,...,Ap — 系数方阵(作为单独参量),矩阵.
e — 特征值,向量.
X — 特征向量,矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy import linalg
import ast
import numpy as np
from sympy import Matrix
from scipy import linalg
def polynomial_eigen_value(input_str):
"""
求解多项式特征值问题：(A0 + λ A1 + ... + λ^p Ap)x = 0
参数:
input_str (str): 表示矩阵元组的字符串，例如 "([[1,0], [0,1]], [[2,1], [1,2]])"
返回:
tuple: (特征向量矩阵, 特征值数组) 或错误信息字符串
示例:
>>> polynomial_eigen_value("([[1,0], [0,1]], [[0,1], [1,0]])")
(Matrix([
[ 1.0, -1.0],
[ 1.0, 1.0]]), Matrix([[-1.0], [1.0]]))
"""
try:
# 将输入字符串转换为Python元组
matrices_tuple = ast.literal_eval(input_str)
# 验证输入是否为非空元组
if not isinstance(matrices_tuple, tuple) or len(matrices_tuple) == 0:
return "输入错误: 需要矩阵元组，例如 ([[1,0],[0,1]], [[0,1],[1,0]])"
# 转换为NumPy数组并验证矩阵
np_matrices = []
n = None
for mat in matrices_tuple:
# 转换为SymPy矩阵进行格式验证
sym_mat = Matrix(mat)
# 验证是否为方阵
if sym_mat.rows != sym_mat.cols:
return f"矩阵错误: 矩阵{sym_mat}不是方阵"
# 记录第一个矩阵的维度
if n is None:
n = sym_mat.rows
elif sym_mat.rows != n:
return f"维度不一致: 所有矩阵应为{n}x{n}"
# 转换为NumPy数组用于计算
np_matrices.append(np.array(sym_mat.tolist(), dtype=np.float64))
# 多项式阶数验证
if len(np_matrices) < 1:
return "需要至少一个矩阵"
def solve_polyeig(*A):
"""
核心求解函数，将多项式特征值问题转化为广义特征值问题
算法参考MATLAB polyeig实现
"""
# 矩阵维度一致性验证
for Ai in A:
if Ai.shape != A[0].shape:
raise ValueError("所有矩阵必须具有相同维度")
n = A[0].shape[0]
p = len(A) - 1 # 多项式阶数
# 构造广义特征值问题矩阵
if p == 0: # 标准特征值问题
return linalg.eig(A[0])
else:
# 构建块矩阵
C = np.zeros((n * p, n * p))
np.fill_diagonal(C[:n * (p - 1), n:], 1) # 填充次对角线
C[-n:, :] = -np.hstack(A[:p]) # 最后一行填充前p个矩阵
D = np.eye(n * p)
D[-n:, -n:] = A[p] # 右下角块为最后一个矩阵
# 求解广义特征值问题
eigvals, eigvecs = linalg.eig(C, D)
# 提取有效特征向量
valid_indices = ~np.isinf(eigvals)
eigvals = eigvals[valid_indices]
eigvecs = eigvecs[:n, valid_indices] # 取前n行
# 排序：按实部升序排列
sort_idx = np.argsort(eigvals.real)
return eigvecs[:, sort_idx], eigvals[sort_idx]
# 执行求解
eigenvectors, eigenvalues = solve_polyeig(*np_matrices)
# 转换为SymPy矩阵返回
return (Matrix(eigenvectors.astype(np.float64)),
Matrix(eigenvalues.astype(np.float64)))
except (SyntaxError, ValueError) as e:
return f"输入解析错误: {str(e)}"
except np.linalg.LinAlgError as e:
return f"数值计算错误: {str(e)}"
except Exception as e:
return f"未知错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 有效输入测试
test_cases = [
("二阶系统",
"([[1,0],[0,1]], [[0,1],[1,0]], [[2,0],[0,2]])")
]
for name, input_str in test_cases:
print(f"{name}测试:")
print(polynomial_eigen_value(input_str))
print()
# (Matrix([[0.577350269189626, 0.577350269189626, -0.577350269189625, -0.577350269189625],
# [0.577350269189626, 0.577350269189626, 0.577350269189626, 0.577350269189626]]),
#
# Matrix([[-0.25],
# [-0.25],
# [ 0.25],
# [ 0.25]]))
polyfit 多项式曲线拟合
多项式曲线拟合
p = polyfit(x,y,n) 返回次数为n的多项式 p(x) 的系数，该阶数是y中数据的最佳拟合（基于最小二乘指标）.p中的系数按降幂排列,p的长度为 n+1
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是多项式拟合的次数, 正整数标量.
Y = p1*x+p2
Y = p1*x^2+p2*x+p3
Y = p1*x^3+p2*x^2+...+p4
Y = p1*x^9+p2*x^8+...+p10
以此类推，最大到poly9
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import numpy as np
def polynomial_fit(input_str):
"""
对标MATLAB的polyfit函数，实现多项式曲线拟合，返回降序排列的系数列表。
参数:
input_str (str): 输入字符串，格式应为 "([x1,x2,...], [y1,y2,...], n)"，
其中n为多项式阶数。
返回:
list: 拟合后的多项式系数列表，按降幂排列（如 [an, an-1, ..., a0]）。
str: 如果输入错误或拟合失败，返回错误信息。
"""
try:
# 解析输入字符串为Python对象（元组）
expr = ast.literal_eval(input_str)
# 验证输入格式是否为 (x_list, y_list, n)
if not (isinstance(expr, tuple) and len(expr) == 3):
return f"输入格式错误: 应为三元组 (x列表, y列表, 阶数)，例如 '([1,2,3], [4,5,6], 2)'"
x_list, y_list, n = expr
# 进一步验证子元素类型
if not (isinstance(x_list, list) and isinstance(y_list, list) and isinstance(n, int)):
return "输入错误: 前两个元素应为列表，第三个为整数"
if len(x_list) != len(y_list):
return f"输入错误: x与y数据长度不一致（{len(x_list)} vs {len(y_list)}）"
if n < 0 or n >= len(x_list):
return f"阶数错误: n应为非负整数且小于数据点数（当前n={n}, 数据点数={len(x_list)}）"
# 转换为NumPy数组并确保为一维
x_array = np.array(x_list, dtype=float).ravel()
y_array = np.array(y_list, dtype=float).ravel()
# 执行多项式拟合（系数按降幂排列）
coefficients = np.polyfit(x_array, y_array, deg=n)
# 转换为列表并保留小数精度（可选）
return coefficients.round(6).tolist()
except SyntaxError:
return "语法错误: 输入字符串格式不正确，请检查括号和逗号"
except ValueError as ve:
return f"数值错误: {str(ve)}"
except Exception as e:
return f"未知错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例1: 二次拟合
input1 = "([1, 2, 3, 4], [1, 4, 9, 16], 2)"
print("示例1 拟合结果:", polynomial_fit(input1))
# [1.0, -0.0, 0.0]
# 示例2: 线性拟合
input2 = "([0, 1, 2], [1, 3, 5], 1)"
print("示例2 拟合结果:", polynomial_fit(input2))
# [2.0, 1.0]
polyijfit 多项式ij曲线拟合
多项式ij曲线拟合
对于多项式曲面，模型名称为'polyij',其中i是x的次数,j是y的次数.i和j的最大值均为5.多项式的次数是i和j的最大值.
每项中x的次数将小于等于i,每项中y的次数将小于等于j.
p = polyijfit(x,y,n) 返回次数为n的多项式 p(x) 的系数，该阶数是y中数据的最佳拟合（基于最小二乘指标）.p中的系数按降幂排列,p的长度为 n+1
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是多项式拟合的次数, 正整数标量.
poly21: Z = p00 + p10*x + p01*y + p20*x^2 + p11*x*y
poly13: Z = p00 + p10*x + p01*y + p11*x*y + p02*y^2 + p12*x*y^2 + p03*y^3
poly55: Z = p00 + p10*x + p01*y +...+ p14*x*y^4 + p05*y^5
# 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 polynomial_ij_fit(input_str):
"""
二维多项式曲线拟合函数，支持多种预设多项式模型
参数:
input_str (str): 输入字符串，格式应为 "(x_list, y_list, z_list, model_type)"
- x_list, y_list: 自变量数据列表
- z_list: 因变量数据列表
- model_type: 预设模型类型，可选 21, 13, 55
返回:
拟合后的多项式表达式字符串
或错误信息字符串
"""
try:
expr = sp.sympify(input_str)
error = False
maxfev = 10000
# 定义多项式拟合函数
def poly21(xy, p00, p10, p01, p20, p11):
x, y = xy
return p00 + p10 * x + p01 * y + p20 * x ** 2 + p11 * x * y
def poly13(xy, p00, p10, p01, p11, p02, p12, p03):
x, y = xy
return p00 + p10 * x + p01 * y + p11 * x * y + p02 * y ** 2 + p12 * x * y ** 2 + p03 * y ** 3
def poly55(xy, p00, p10, p01, p11, p20, p02, p12, p21, p03, p13, p04, p22, p14, p05):
x, y = xy
return (p00 + p10 * x + p01 * y + p11 * x * y + p20 * x ** 2 + p02 * y ** 2 +
p12 * x * y ** 2 + p21 * x ** 2 * y + p03 * y ** 3 + p13 * x * y ** 3 +
p04 * y ** 4 + p22 * x ** 2 * y ** 2 + p14 * x * y ** 4 + p05 * y ** 5)
if isinstance(expr, tuple):
if len(expr) > 2:
x, y, n = expr[0], expr[1], int(expr[2])
else:
x, y = expr[0], expr[1]
n = 1
if isinstance(x, list):
x_data = np.array(x, dtype=float).ravel()
if isinstance(y, list):
y_data = np.array(y, dtype=float).ravel()
else:
error = True
x_data, y_data = np.meshgrid(x_data, y_data)
z_data = 1 + 2 * x_data + 3 * y_data + 4 * x_data ** 2 + 5 * x_data * y_data
z_data += 0.1 * np.random.normal(size=z_data.shape) # 添加一些噪声
z_data = z_data.flatten()
x_data = x_data.flatten()
y_data = y_data.flatten()
xy_data = np.vstack((x_data, y_data))
if n == 21:
# 扁平化输入数据
# 设置初始猜测值
initial_guess = np.ones(poly21.__code__.co_argcount - 1) # 所有参数初始值设置为1
# 进行曲线拟合
params, params_covariance = curve_fit(poly21, xy_data, z_data, p0=initial_guess, maxfev=maxfev)
expression_with_vars = " + ".join(
[f"{params[0]:.4f}"] +
[f"{params[1]:.4f} * x"] +
[f"{params[2]:.4f} * y"] +
[f"{params[3]:.4f} * x^2"] +
[f"{params[4]:.4f} * x * y"]
)
elif n == 13:
initial_guess = np.ones(poly13.__code__.co_argcount - 1) # 所有参数初始值设置为1
# 进行曲线拟合
params, params_covariance = curve_fit(poly13, xy_data, z_data, p0=initial_guess, maxfev=maxfev)
expression_with_vars = " + ".join(
[f"{params[0]:.4f}"] +
[f"{params[1]:.4f} * x"] +
[f"{params[2]:.4f} * y"] +
[f"{params[3]:.4f} * x * y"] +
[f"{params[4]:.4f} * y^2"] +
[f"{params[5]:.4f} * x * y^2"] +
[f"{params[6]:.4f} * y^3"]
)
elif n == 55:
initial_guess = np.ones(poly55.__code__.co_argcount - 1) # 所有参数初始值设置为1
# 进行曲线拟合
params, params_covariance = curve_fit(poly55, xy_data, z_data, p0=initial_guess, maxfev=maxfev)
expression_with_vars = " + ".join(
[f"{params[0]:.4f}"] +
[f"{params[1]:.4f} * x"] +
[f"{params[2]:.4f} * y"] +
[f"{params[3]:.4f} * x * y"] +
[f"{params[4]:.4f} * x^2"] +
[f"{params[5]:.4f} * y^2"] +
[f"{params[6]:.4f} * x * y^2"] +
[f"{params[7]:.4f} * x^2 * y"] +
[f"{params[8]:.4f} * y^3"] +
[f"{params[9]:.4f} * x * y^3"] +
[f"{params[10]:.4f} * y^4"] +
[f"{params[11]:.4f} * x^2 * y^2"] +
[f"{params[12]:.4f} * x * y^4"] +
[f"{params[13]:.4f} * y^5"]
)
else:
error = True
return sp.simplify(expression_with_vars) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例用法
if __name__ == "__main__":
# 生成测试数据 (z = 1 + 2x + 3y + 4x² + 5xy)
x = np.linspace(0, 1, 10)
y = np.linspace(0, 1, 10)
X, Y = np.meshgrid(x, y)
# 转换为列表格式
input_data = (
X.ravel().tolist(),
Y.ravel().tolist(),
21
)
# 执行拟合
print(polynomial_ij_fit(str(input_data)))
# 4.0158*x**2 + 5.0097*x*y + 1.9801*x + 2.9941*y + 1.0065
polyint 多项式积分
多项式积分
q = polyint(p,k) 使用积分常量 k 返回 p 中系数所表示的多项式积分.
q = polyint(p) 假定积分常量 k = 0.
p是多项式系数,指定为向量.
k是积分常量,指定为数值标量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import sympy as sp
def polynomial_indefinite_integral(input_str):
"""
对标MATLAB的polyint函数，计算多项式的不定积分，支持指定积分常数。
参数:
input_str (str): 输入的多项式或元组的字符串表示，如 "[3,2,1]" 或 "([3,2,1],5)"。
- 单个列表：多项式系数，积分常数默认为0。
- 元组（系数列表, k）：积分常数为k。
返回:
list: 积分后多项式的系数列表，包含指定的积分常数。
str: 如果输入错误或发生异常，返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查表达式是否为元组且长度为 2
if isinstance(expr, tuple) and len(expr) == 2:
# 解包元组，得到多项式系数列表和积分常数
coeffs, k = expr
# 检查系数是否为列表且积分常数是否为数字
if isinstance(coeffs, list) and k.is_number:
# 创建一个以 'x' 为变量的 SymPy 多项式对象
poly1 = sp.Poly(coeffs, sp.symbols('x'))
# 对多项式进行积分，并加上指定的积分常数
integrated_poly = poly1.integrate() + float(k)
# 获取积分后多项式的所有系数
result = integrated_poly.all_coeffs()
else:
# 如果系数不是列表或积分常数不是数字，标记输入错误
error = True
# 检查表达式是否为列表
elif isinstance(expr, list):
# 创建一个以 'x' 为变量的 SymPy 多项式对象
poly_expr = sp.Poly(expr, sp.symbols('x'))
# 对多项式进行积分，积分常数默认为 0
integrated_poly = poly_expr.integrate()
# 获取积分后多项式的所有系数
result = integrated_poly.all_coeffs()
else:
# 如果表达式既不是指定格式的元组也不是列表，标记输入错误
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 单个多项式，默认k=0
p1 = "[3,0,-4,10,-25]"
print("示例1:", polynomial_indefinite_integral(p1))
# [3/5, 0, -4/3, 5, -25, 0]
# 示例2: 指定积分常数k=5
p2 = "([3, 2, 1], 5)"
print("示例2:", polynomial_indefinite_integral(p2))
# [1.00000000000000, 1.00000000000000, 1.00000000000000, 5.00000000000000]
# 示例3: 低次多项式测试
p3 = "[1, 2]"
print("示例3:", polynomial_indefinite_integral(p3))
# [1/2, 2, 0]
polylog 多重对数函数
多重对数函数
Li＝polylog(n,x)返回阶数n和自变量x的多对数.
n —— 多对数的阶数, 数字,数组,符号数,符号变量,符号函数,符号表达式,符号数组.
x —— 多对数的自变量, 数字,数组,符号数,符号变量,符号函数,符号表达式,符号数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def polynomial_logarithm(input_str):
"""
计算多对数函数 Polylog(n, x)（对标MATLAB的polylog函数）
参数:
input_str: 输入字符串，格式应为"(n, x)"，其中：
n: 整数标量 或 数值矩阵
x: 数值标量 或 数值矩阵
示例："(2, 0.5)"
"(Matrix([[1,2],[3,4]]), 0.5)"
"(2, Matrix([[0.1,0.2],[0.3,0.4]]))"
返回:
SymPy 矩阵/浮点数: 计算结果
str: 错误信息
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 计算结果
def elementwise_polylog(n_val, x_val):
"""元素级polylog计算"""
try:
return sp.polylog(n_val, x_val).evalf()
except:
return sp.nan # 无效计算返回nan
if isinstance(expr, tuple) and len(expr) == 2:
n, x = expr
result = elementwise_polylog(n, x)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 测试示例
if __name__ == "__main__":
# 正常案例
print("标量计算:")
print(polynomial_logarithm("(2, 0.5)"))
# 0.582240526465013
print(polynomial_logarithm("(a, 0.5)"))
# polylog(a, 0.5)
polymul 多项式乘积
多项式乘积
k＝polymul(f,g,...)返回多个多项式相乘,并且可以有多个变量.
f, g -- 输入, 符号表达式.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def poly_multivariate_product(input_str):
"""
多变量多项式的polymvp函数，实现多项式乘法，返回乘积多项式的系数列表（降幂排列）。
参数:
input_str (str): 输入字符串，表示多个多项式的元组，例如 "(x^2+1,x+2,x-2)"。
返回:
str: 乘积多项式表达式。
str: 如果输入错误，返回错误信息。
"""
try:
# 使用 sympy 的 sympify 函数将输入字符串转换为 SymPy 表达式，evaluate=False 表示不进行计算
expr = sp.sympify(input_str, evaluate=False)
# 初始化错误标志为 False，表示目前没有错误
error = False
# 初始化结果为 1，用于后续多项式相乘
result = 1
# 检查转换后的表达式是否为元组类型
if isinstance(expr, tuple):
# 遍历元组中的每个元素
for item in expr:
# 将当前元素转换为 SymPy 的多项式对象，并与结果相乘
result *= sp.Poly(item)
else:
# 如果表达式不是元组类型，将错误标志设为 True
error = True
return result.as_expr() if (not error or result != 1) else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 测试示例
if __name__ == "__main__":
# 正常案例
print("单变量多项式:")
print(poly_multivariate_product("(x**2+1,x+2,x-2)"))
# x**4 - 3*x**2 - 4
print("\n多变量多项式:")
print(poly_multivariate_product("(x**2+3*x+4,x+y+2)"))
# x**3 + x**2*y + 5*x**2 + 3*x*y + 10*x + 4*y + 8
polynomialDegree 多项式次数
多项式次数
polynomialDegree(p)通过vars返回多项式p相对于p中所有变量的次数
polynomialDegree(p，vars)返回p相对于vars中变量的阶数.
p - 多项式,符号表达式.
vars - 多项式变量,符号变量向量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Poly, Symbol, sympify
from ast import literal_eval
def polynomial_of_degree(input_str):
"""
对标 Wolfram 的 PolynomialDegree 函数，计算多项式的最高次数。
参数:
input_str (str): 输入字符串，支持两种格式：
1. 单个多项式表达式 (计算所有变量的总次数)
示例: "3*x**2 + 2*y**3 + 1"
2. 元组 (多项式表达式, 变量列表)
示例: "(3*x**2 + 2*y**3 + 1, ['x', 'y'])"
示例: "(x**2*y + x*y**2, [x, y])"
返回:
int: 多项式的最高次数
str: 错误信息（若输入无效）
"""
try:
# 解析输入字符串
expr = sp.sympify(input_str)
error = False
result = None
# 处理元组输入（多项式 + 变量列表）
if isinstance(expr, tuple) and len(expr) == 2:
# 处理元组输入（多项式 + 变量列表）
poly_expr, variables = expr
if isinstance(variables, list):
sym_vars = variables
elif variables.free_symbols:
sym_vars = [variables]
else:
error = True
poly = sp.poly(poly_expr, *sym_vars)
result = poly.total_degree()
# 处理单个多项式输入（自动检测所有变量）
elif expr.free_symbols:
expr_poly = sp.poly(expr)
result = expr_poly.total_degree()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except SyntaxError:
return "语法错误: 输入字符串格式不正确"
except ValueError as ve:
return f"值错误: {str(ve)}"
except TypeError as te:
return f"类型错误: {str(te)}"
except Exception as e:
return f"未知错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例1: 计算多变量多项式的总次数
input1 = "3*x**2 + 2*y**3 + 1"
print(f"示例1 输入: {input1}")
print("输出:", polynomial_of_degree(input1))
# 3
# 示例2: 指定变量子集计算次数
input2 = "(x**2*y + x*y**2, [x, y])"
print(f"\n示例2 输入: {input2}")
print("输出:", polynomial_of_degree(input2))
# 3
polyval 多项式计算
多项式计算
y = polyval(p,x) 计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）
p是多项式系数，指定为向量.
x是查询点,指定为向量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def polynomial_value(input_str):
"""
计算多项式在给定点的值，对标MATLAB的polyval函数。
参数：
input_str (str): 输入字符串，表示多项式系数和计算点。格式应为SymPy可解析的元组，
例如："([3, 0, 1], 2)" 或 "([1, -2], [5, 6])"。
返回：
float/list/str: 计算结果。若输入错误或计算异常，返回错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查表达式是否为元组且长度为2
if isinstance(expr, tuple) and len(expr) == 2:
# 从元组中提取多项式系数部分和计算点部分
c_part, x_part = expr[0], expr[1]
# 使用numpy的poly1d函数根据多项式系数创建多项式对象
p = np.poly1d(c_part)
# 检查计算点部分是否为列表
if isinstance(x_part, list):
# 如果是列表，对列表中的每个元素计算多项式的值
result = [np.polyval(p, item) for item in x_part]
# 检查计算点部分是否为一个数字
elif x_part.is_number:
# 如果是数字，计算该点处多项式的值
result = np.polyval(p, x_part)
else:
# 如果计算点既不是列表也不是数字，标记为输入错误
error = True
else:
# 如果表达式不是长度为2的元组，标记为输入错误
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：计算多项式 3x² + 1 在x=2处的值
input_str1 = "([3, 0, 1], 2)"
print(f"输入: {input_str1} -> 结果: {polynomial_value(input_str1)}")
# 13
# 示例2：计算多项式 x - 2 在x=5和x=6处的值
input_str2 = "([1, -2], [5, 6])"
print(f"输入: {input_str2} -> 结果: {polynomial_value(input_str2)}")
# [3, 4]
# 示例3：符号输入（符号系数）
input_str3 = "([a, 1], 2)"
print(f"输入: {input_str3} -> 结果: {polynomial_value(input_str3)}")
# 2*a + 1
polyvalm 矩阵多项式计算
矩阵多项式计算
Y = polyvalm(p,X) 以矩阵方式返回多项式p的计算值.此计算方式等同于使用多项式p替换矩阵X
p是多项式系数，指定为向量.
x是方阵,指定为矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def polynomial_val_matrix(input_str):
"""
计算矩阵多项式的值，对标MATLAB的polyvalm函数。
参数:
input_str (str): 输入字符串，格式应为 "([c0, c1, ..., cn], [[m11, m12,...], ...])"，
其中系数列表表示多项式，矩阵为方阵。
返回:
SymPy矩阵/str: 计算结果。若输入错误或计算异常，返回错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def eval_poly_matrix(P, A):
"""
核心计算函数：计算矩阵多项式的值 P(A)
参数:
P (sp.Poly): 多项式对象
A (sp.Matrix): 方阵
返回:
sp.Matrix: 计算结果矩阵
"""
# 初始化结果矩阵为同型零矩阵
res = sp.zeros(*A.shape)
# 获取系数列表（从最高次到常数项）并反转，以便从x^0开始处理
coefficients = P.all_coeffs()[::-1]
for power, coeff in enumerate(coefficients):
res += coeff * (A ** power) # 累加每一项：coeff * A^power
return res
# 检查输入是否为包含两个元素的SymPy元组
if isinstance(expr, tuple) and len(expr) == 2:
coeff_part, matrix_part = expr[0], expr[1]
# 处理多项式系数部分
if isinstance(coeff_part, list):
try:
# 将SymPy列表转换为Python浮点数列表
coefficients = [float(c) for c in coeff_part]
except TypeError:
error = True # 包含非数值型系数
else:
error = True # 系数部分不是列表
# 处理矩阵部分
A = sp.Matrix(matrix_part) if isinstance(matrix_part, list) else None
if A is None:
error = True # 矩阵无效或非方阵
elif A.is_square == False:
error = True
if not error:
try:
# 创建多项式对象（例如 [3,0,1] -> 3x² + 1）
poly = sp.Poly(coefficients, sp.symbols('x'))
result = eval_poly_matrix(poly, A)
except Exception as e:
return f"计算错误: {e}"
else:
error = True
else:
error = True # 输入结构不符合要求
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：计算 3A² + I，其中A = [[1,2],[3,4]]
input_str1 = "([3,0,1],[[1,2],[3,4]])"
A = sp.Matrix([[1, 2], [3, 4]])
print(f"输入: {input_str1}")
print("函数计算结果:\n", polynomial_val_matrix(input_str1))
# Matrix([[22.0000000000000, 30.0000000000000],
# [45.0000000000000, 67.0000000000000]])
PolygonAngle 多边形顶点角度
多边形顶点角度
Angles = PolygonAngle([(x,y),(x1,y1),(x2,y2)...])计算多边形各顶点内角.
x, y -- 数值，多边形顶点的平面坐标。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import math
import ast
def polygon_area(points):
"""计算多边形有向面积以判断顶点方向（正数为逆时针，负数为顺时针）"""
area = 0.0
n = len(points)
for i in range(n):
x1, y1 = points[i]
x2, y2 = points[(i + 1) % n]
area += (x2 - x1) * (y2 + y1) # Shoelace公式简化版
return area / 2.0
def vector_angle(v1, v2, is_ccw):
"""计算两个向量间的内角（考虑多边形方向）"""
# 计算向量长度
len1 = math.hypot(*v1)
len2 = math.hypot(*v2)
# 处理零向量
if len1 == 0 or len2 == 0:
return 0.0
# 计算点积和叉积
dot = (v1[0] * v2[0] + v1[1] * v2[1]) / (len1 * len2)
cross = (v1[0] * v2[1] - v1[1] * v2[0]) / (len1 * len2)
# 基础夹角（0~π）
base_angle = math.acos(max(min(dot, 1), -1)) # 防止浮点误差
# 根据叉积符号判断转向
if is_ccw:
return base_angle if cross > 0 else 2 * math.pi - base_angle
else:
return base_angle if cross < 0 else 2 * math.pi - base_angle
def shapely_polygon_angles(input_str):
"""
计算多边形各顶点内角（单位：度）
输入格式示例：[[x1,y1], [x2,y2], ...]
"""
try:
# 输入解析与校验
polygon = ast.literal_eval(input_str)
if not isinstance(polygon, (list, tuple)) or len(polygon) < 3:
return "错误：需要至少3个顶点的列表"
# 转换坐标并拷贝首节点处理环形结构
points = [tuple(map(float, p)) for p in polygon]
points.append(points[0]) # 闭合多边形
# 计算多边形方向
area = polygon_area(points[:-1])
is_ccw = area > 0
# 计算各顶点角度
angles = []
for i in range(len(points) - 1):
# 获取三点坐标
a, b, c = points[i], points[i + 1], points[(i + 2) % (len(points) - 1)]
# 计算向量
v_prev = (a[0] - b[0], a[1] - b[1]) # BA向量
v_next = (c[0] - b[0], c[1] - b[1]) # BC向量
# 计算带方向的内角
angle = vector_angle(v_prev, v_next, is_ccw)
angles.append(math.degrees(angle))
return angles
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 案例1：正六边形（逆时针）
hexagon = [[1, 0], [0.5, 0.866], [-0.5, 0.866],
[-1, 0], [-0.5, -0.866], [0.5, -0.866]]
print("正六边形角度:", shapely_polygon_angles(str(hexagon)))
# [120.00072778082736, 120.00072778082736, 119.99854443834525, 120.00072778082736, 120.00072778082736, 119.99854443834525]
# 案例2：凹多边形（星形）
star = [(0, 0), (2, 3), (5, 0), (2, -3), (-1, 0)]
print("\n凹多边形角度:", shapely_polygon_angles(str(star)))
# [78.69006752597979, 90.0, 90.0, 45.00000000000001, 236.30993247402026]
PolygonArea 多边形的面积
多边形的面积
PolygonArea([[x1,y1], [x2,y2], ..., [xn,yn]])返回多边形的面积。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import shapely.geometry as sh
def shapely_polygon_area(input_str):
"""
计算多边形的面积，对标Wolfram的PolygonArea函数。
参数：
input_str (str): 表示多边形顶点列表的字符串，格式应为[[x1,y1], [x2,y2], ..., [xn,yn]]或类似结构。
返回：
float/str: 计算成功返回面积（浮点数），否则返回错误信息字符串。
示例：
>>> shapely_polygon_area("[[0,0], [0,1], [1,1]]")
0.5
>>> shapely_polygon_area("[(0,0), (1,0), (1,1), (0,1)]")
1.0
"""
try:
# 解析输入字符串为Python对象（列表或元组）
vertices = ast.literal_eval(input_str)
# 检查输入是否为顶点集合
if not isinstance(vertices, (list, tuple)):
return f"输入错误: 需要顶点列表，例如 [[0,0], [0,1], [1,1]]，当前输入类型为 {type(vertices).__name__}"
# 验证每个顶点的格式
for i, point in enumerate(vertices):
# 检查是否为二维坐标
if not isinstance(point, (list, tuple)) or len(point) != 2:
return f"顶点格式错误: 第{i + 1}个顶点{point}应为包含两个数值的列表/元组"
# 检查坐标值是否为数字
try:
x, y = float(point[0]), float(point[1])
except (TypeError, ValueError):
return f"坐标值错误: 第{i + 1}个顶点{point}包含非数值类型"
# 创建多边形并计算面积
polygon = sh.Polygon(vertices)
return abs(polygon.area) # 返回绝对值以保持与Wolfram一致
except (SyntaxError, ValueError) as e:
return f"解析错误: 输入格式无效 - {str(e)}"
except Exception as e:
return f"计算错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
test_cases = [
("[[0,0], [0,1], [1,1]]"),
# 0.5
("[(0,0), (1,0), (1,1), (0,1)]"),
# 1.0
("[[0,0], [0,2], [2,2], [2,0]]"),
# 4.0
]
for input_str in test_cases:
result = shapely_polygon_area(input_str)
print(f"输入: {input_str}")
print(f"实际: {result}\n")
PolygonDiffArea 两个几何体的差异面积
两个几何体的差异面积
PolygonDiffArea(A, B, type)返回两个多边形AB的交集面积。
type为不同的交集面积:
0 -- A-B的差集面积
1 -- A∩B的交集面积
2 -- AΔB的对称差面积
3 -- A∪B的并集面积
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import shapely.geometry as sh
from shapely.errors import TopologicalError
def shapely_difference_area(input_str):
"""
计算两个几何图形的指定操作面积（支持点/线/多边形）
参数：
input_str (str): 格式为 "([[x1,y1],...], [[x1,y1],...], 操作码)" 的字符串，其中操作码：
0 -> A-B的差集面积
1 -> A∩B的交集面积
2 -> AΔB的对称差面积
3 -> A∪B的并集面积
返回：
float/str: 成功返回面积值，失败返回错误描述
示例：
>>> shapely_difference_area('([[0,0],[0,1],[1,1]], [[0.5,0.5],[0.5,1.5],[1.5,1.5]], 0)')
0.25
"""
try:
# 解析输入字符串
data = ast.literal_eval(input_str)
# 验证输入结构
if not (isinstance(data, tuple) and len(data) == 3):
return "输入格式错误：需要三元组 (几何体A, 几何体B, 操作码)"
a_points, b_points, operation = data
operation = int(operation) # 确保操作码是整数
# 定义几何体构建函数
def build_geometry(points):
"""将顶点列表转换为Shapely几何对象"""
if not isinstance(points, (list, tuple)):
raise ValueError("顶点数据必须是列表/元组")
# 验证单个几何体的顶点数据
if len(points) == 0:
raise ValueError("空几何体")
for i, p in enumerate(points):
if not isinstance(p, (list, tuple)) or len(p) != 2:
raise ValueError(f"顶点{i}格式错误：{p}")
try:
float(p[0]), float(p[1])
except ValueError:
raise ValueError(f"顶点{i}包含非数值坐标：{p}")
# 根据顶点数量创建几何对象
if len(points) == 1:
return sh.Point(points[0])
elif len(points) == 2:
return sh.LineString(points)
else:
return sh.Polygon(points)
# 创建几何对象
try:
a = build_geometry(a_points)
b = build_geometry(b_points)
except ValueError as e:
return f"几何体错误：{str(e)}"
# 执行几何操作
try:
if operation == 0:
result = a.difference(b)
elif operation == 1:
result = a.intersection(b)
elif operation == 2:
result = a.symmetric_difference(b)
elif operation == 3:
result = a.union(b)
else:
return f"操作码错误：支持0-3，当前是{operation}"
except TopologicalError:
return "几何操作失败：请检查几何体有效性"
return abs(result.area) # 返回绝对值保持结果稳定
except (SyntaxError, ValueError) as e:
return f"解析错误：{str(e)}"
except Exception as e:
return f"未知错误：{str(e)}"
# 测试案例
if __name__ == "__main__":
tests = [
# 有效输入测试
('([[0,0],[0,1],[1,1]], [[0.5,0.5],[0.5,1.5],[1.5,1.5]], 0)'),
# 0.375
('([(0,0),(1,1)], [[0,0],[0,1],[1,0]], 1)'),
# 0.0
('([[0,0],[0,2],[2,2],[2,0]], [[1,1],[1,3],[3,3]], 2)'),
# 5.0
('([[0,0],[0,1],[1,1]], [[0.5,0],[0.5,1],[1.5,1]], 3)'),
# 0.875
]
for input_str in tests:
print(f"输入: {input_str}")
result = shapely_difference_area(input_str)
print(f"实际: {result}\n")
PolygonLength 多边形的周长
多边形的周长
PolygonLength(A)返回多边形A的周长
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import ast
import shapely.geometry as sh
from shapely.errors import TopologicalError
def shapely_polygon_length(input_str):
"""
计算多边形的周长。
参数:
input_str (str): 表示多边形顶点的字符串，格式应为列表，例如 "[(0,0), (0,1), (1,1)]"。
返回:
float/str: 成功时返回多边形周长，失败时返回错误信息字符串。
示例:
>>> shapely_polygon_length("[(0,0), (0,1), (1,1)]")
3.414213562373095
"""
try:
# 将输入字符串转换为Python列表
coords = ast.literal_eval(input_str)
# 验证输入格式
if not (isinstance(coords, list) and
all(isinstance(p, (tuple, list)) and len(p) == 2 for p in coords)):
return f"输入格式错误: 需要二维坐标点列表，例如 [(x1,y1), (x2,y2), ...]"
# 创建多边形对象（Shapely会自动闭合多边形）
polygon = sh.Polygon(coords)
# 返回计算后的周长
return polygon.length
except (SyntaxError, ValueError) as e:
# 处理格式解析错误
return f"输入解析错误: {str(e)}"
except TopologicalError as e:
# 处理无效几何图形（如自相交）
return f"几何错误: {str(e)}"
except Exception as e:
# 处理其他未知错误
return f"计算错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 有效输入测试
test_cases = [
("三角形", "[(0,0), (0,1), (1,1)]"),
# 3.414213562373095
("正方形", "[(0,0), (0,1), (1,1), (1,0)]"),
# 4.0
("五边形", "[(0,0), (1,0), (2,1), (1,2), (0,1)]")
# 6.242640687119285
]
for name, input_str in test_cases:
result = shapely_polygon_length(input_str)
print(f"输入: {input_str}")
print(f"输出: {result}\n")
polydiv 多项式的商和余数
多项式的商和余数
q, r = polydiv(f,g)返回分子f和分母g多项式相除的商和余数.分子和分母都应该是多项式函数.返回的结果是一个列表,其中第一项为商,第二项为余数.
f, g -- 输入, 符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def polynomial_division(input_str):
"""
执行多项式长除法运算
参数:
input_str (str): 表示多项式对的字符串，格式为"(被除式, 除式)"，例如："(x**2 + 1, x - 1)"
返回:
tuple/str: 成功时返回(商, 余式)元组，失败时返回错误信息字符串
示例:
>>> polynomial_division("(x**3 - 1, x - 1)")
(x**2 + x + 1, 0)
>>> polynomial_division("(2*x**4 + 3*x, x**2 + 1)")
(2*x**2 - 2, 3*x + 2)
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str, evaluate=False)
# 验证输入格式
if not (isinstance(expr, tuple) and len(expr) == 2):
return f"输入格式错误: 需要包含两个多项式的元组，例如 (被除式, 除式)"
# 转换为多项式对象
dividend = sp.Poly(expr[0])
divisor = sp.Poly(expr[1])
# 检查除数有效性
if divisor.degree() < 0:
return "错误: 除式不能为零多项式"
# 执行多项式除法
quotient, remainder = dividend.div(divisor)
# 转换为标准表达式形式返回
return (quotient.as_expr(), remainder.as_expr())
except sp.SympifyError:
return f"表达式解析失败: 请检查输入格式是否正确"
except sp.PolynomialError as pe:
return f"多项式错误: {str(pe)}"
except Exception as e:
return f"未知错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 有效输入测试
x = sp.symbols('x')
test_cases = [
("一次式相除", "(x**2 - 1, x - 1)"),
# (x + 1, 0)
("高次式相除", "(x**5 + 1, x**2 + 1)"),
# (x**3 - x, x + 1)
("带余数情况", "(2*x**3 + 3*x + 1, x**2 + 1)"),
# (2*x, x + 1)
("低次除高次", "(x + 1, x**2 + 1)")
# (0, x + 1)
]
for name, input_str in test_cases:
print(f"{name}测试:")
result = polynomial_division(input_str)
if isinstance(result, tuple):
print(f"结果: {result}\n")
else:
print(f"错误结果: {result}\n")
polygcd 多项式的最大公约数
多项式的最大公约数
p = polygcd(a,b)返回a,b两个表达式之间的最大公约数
a, b -- 输入, 符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def polynomial_gcd(input_str):
"""
计算两个多项式的最大公约数（GCD）（对标MATLAB的gcd函数）
参数:
input_str (str): 表示两个多项式的字符串，格式应为"(poly1, poly2)"
示例："(x**2 - 1, x + 1)"
返回:
sp.Expr: 多项式GCD的符号表达式
str: 错误信息（当输入无效时）
示例:
>>> polynomial_gcd("(x**2 - 1, x + 1)")
x + 1
>>> polynomial_gcd("(x**3 - 1, x**2 - 1)")
x - 1
"""
try:
# 尝试解析输入字符串为SymPy表达式
expr = sp.sympify(input_str, evaluate=False)
# 检查输入是否为包含两个元素的元组
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入错误：需要两个多项式组成的元组，当前输入：{input_str}"
# 尝试将两个表达式转换为多项式
poly1 = sp.Poly(expr[0])
poly2 = sp.Poly(expr[1])
# 检查多项式是否有效
if poly1.is_number or poly2.is_number:
return "错误：输入必须为符号多项式，不能仅为数值"
# 计算最大公约数
gcd_poly = sp.gcd(poly1, poly2)
# 返回化简后的表达式形式
return gcd_poly.as_expr()
except sp.SympifyError:
return f"解析错误：无法解析输入字符串 - {input_str}"
except ValueError as ve:
return f"多项式转换错误：{str(ve)}"
except Exception as e:
return f"运行时错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 正常案例
test_cases = [
("(x**2 - 1, x + 1)"),
# x + 1
("(x**3 - 1, x**2 - 1)"),
# x - 1
("(2*x**2 + 4*x, x**2 - 4)")
# x + 2
]
print("正常案例测试：")
for input_str in test_cases:
result = polynomial_gcd(input_str)
print(f"输入：{input_str}")
print(f"实际：{result}\n")
polylcm 多项式的最小公倍数
多项式的最小公倍数
p = polylcm(a,b)返回a,b两个表达式之间的最小公倍数
a, b -- 输入, 符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def polynomial_lcm(input_str):
"""
计算两个多项式的最小公倍数（LCM）（对标MATLAB的lcm函数）
参数:
input_str (str): 表示两个多项式的字符串，格式应为"(poly1, poly2)"
示例："(x**2 - 1, x + 1)"
返回:
sp.Expr: 多项式LCM的符号表达式
str: 错误信息（当输入无效时）
示例:
>>> polynomial_lcm("(x**2 - 1, x + 1)")
x**2 - 1
>>> polynomial_lcm("(x**3 - 1, x**2 - 1)")
x**4 - x**3 - x + 1
"""
try:
# 解析输入字符串为SymPy表达式（禁用自动化简保持原始形式）
expr = sp.sympify(input_str, evaluate=False)
# 检查输入是否为包含两个元素的元组
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入错误：需要两个多项式组成的元组，当前输入：{input_str}"
# 尝试将两个表达式转换为多项式
try:
poly1 = sp.Poly(expr[0])
poly2 = sp.Poly(expr[1])
except (sp.PolynomialError, ValueError) as e:
return f"多项式转换错误：{str(e)}"
# 检查是否为有效多项式
if poly1.is_number and poly2.is_number:
return sp.lcm(poly1.as_expr(), poly2.as_expr())
# 计算最小公倍数
lcm_poly = sp.lcm(poly1, poly2)
# 返回标准化的表达式形式
return sp.simplify(lcm_poly.as_expr())
except sp.SympifyError:
return f"解析错误：无法解析输入字符串 - {input_str}"
except Exception as e:
return f"运行时错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 正常案例
test_cases = [
("(x**2 - 1, x + 1)"),
# x**2 - 1
("(x**3 - 1, x**2 - 1)"),
# x**4 + x**3 - x - 1
("(2*x**2 + 4*x, x**2 - 4)")
# 2*x*(x**2 - 4)
]
print("正常案例测试：")
for input_str in test_cases:
result = polynomial_lcm(input_str)
print(f"输入：{input_str}")
print(f"实际：{result}\n")
potential 向量场的势
向量场的势
potential(V,X) 计算矢量场 V 相对于笛卡尔坐标系中的矢量 X 的势. 矢量场 V 必须是梯度场.
potential(V,X,Y) 使用 Y 作为积分基点来计算矢量场 V 相对于 X 的势.
V — 矢量场,符号表达式或函数的三维符号向量(默认)
X — 输入,三个符号变量的向量
Y — 输入,数值向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.vector import CoordSys3D, scalar_potential
def potential_vector_field(input_str):
"""
计算向量场的标量势（对标MATLAB的potential函数）
参数格式要求：
输入应为字符串，格式为：
- 二维场： (([F1, F2], [x, y]) 或 (([F1, F2], [x, y], [x0, y0])
- 三维场： (([F1, F2, F3], [x, y, z]) 或 (([F1, F2, F3], [x, y, z], [x0, y0, z0]))
返回：
- 成功：标量势表达式（带积分常数）或带基点的数值结果
- 失败：错误信息字符串
"""
try:
expr = sp.sympify(input_str, evaluate=False)
error = False
result = None
# 输入格式验证
if not isinstance(expr, tuple) or len(expr) not in (2, 3):
return f"输入错误：需要格式为 (分量列表, 变量列表) 或 (分量列表, 变量列表, 基点)，当前输入：{input_str}"
# 解析输入参数
V_components = expr[0] # 向量场分量
X_vars = expr[1] # 变量列表
base_points = expr[2] if len(expr) == 3 else None # 基点（可选）
# 维度验证
dim = len(X_vars)
if dim not in (2, 3) or len(V_components) != dim:
return f"维度错误：支持2D/3D向量场，当前输入维度：分量数={len(V_components)}, 变量数={dim}"
# 创建坐标系
N = CoordSys3D('N') # 使用3D坐标系处理2D/3D场
# 构建向量场表达式
subs_dict = {X_vars[0]: N.x, X_vars[1]: N.y}
vec_components = [
V_components[0].subs(subs_dict),
V_components[1].subs(subs_dict)
]
if dim == 3:
subs_dict[X_vars[2]] = N.z
vec_components.append(V_components[2].subs(subs_dict))
# 创建向量场对象
if dim == 2:
V = vec_components[0] * N.i + vec_components[1] * N.j
else:
V = vec_components[0] * N.i + vec_components[1] * N.j + vec_components[2] * N.k
# 计算标量势
try:
phi = scalar_potential(V, N)
except ValueError as e:
return sp.nan
# 还原原始变量符号
reverse_subs = {N.x: X_vars[0], N.y: X_vars[1]}
if dim == 3:
reverse_subs[N.z] = X_vars[2]
phi = phi.subs(reverse_subs)
# 处理基点（修正部分）
if base_points:
if len(base_points) != dim:
return f"基点维度错误：需要{dim}个基点值，当前收到{len(base_points)}"
# 计算基点处的势值
phi_base = phi.subs(dict(zip(X_vars, base_points)))
# 消除积分常数，使基点处的势为0
phi = phi - phi_base
return phi
except sp.SympifyError:
return f"解析错误：无效输入格式 - {input_str}"
except Exception as e:
return f"运行时错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 有效案例
# 案例1：二维保守场
print("案例1：二维保守场")
input_2d = "([x,y,z*exp(z)],[x,y,z])"
print(potential_vector_field(input_2d))
# x**2/2 + y**2/2 + (z - 1)*exp(z)
# 案例2：三维带基点
print("\n案例2：三维带基点计算")
input_3d_base = "([x,y,z*exp(z)],[x,y,z],[0,0,0])"
print(potential_vector_field(input_3d_base))
# x**2/2 + y**2/2 + (z - 1)*exp(z) + 1
# 案例3：非保守场
print("\n案例3：非保守场")
input_nonconservative = "(([y, x + y], [x, y]))"
print(potential_vector_field(input_nonconservative))
# x*y + y**2/2
powTwo 浮点数的以2为底的幂运算和缩放
浮点数的以2为底的幂运算和缩放
Y = pow2(E) 计算2的E次方
Y = pow2(X,E) 计算X乘以2的E次方
E是标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import exp2
def pow_2_e(input_str):
"""
对标 MATLAB 的 pow2 函数，实现两种计算模式：
模式1：pow2(X) = 2.^X
模式2：pow2(F,E) = F.*2.^E
参数:
input_str: 输入字符串，支持以下格式：
- 单个数值/矩阵: "2" 或 "Matrix([[1,2],[3,4]])"
- 元组模式: "(F, E)"，其中F和E可为标量或矩阵
返回:
SymPy 矩阵/数值 或 错误信息字符串
"""
try:
# 解析输入字符串为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
f, e = params
result = f * exp2(e)
elif any(e.free_symbols for e in expr):
"""元素级计算 f * 2^e"""
f, e = expr
try:
result = f * sp.Pow(2, e).evalf()
except:
return sp.nan
else:
error = True
elif expr.is_number:
z = float(expr)
result = exp2(z)
elif expr.free_symbols:
result = sp.Pow(2, expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except sp.SympifyError:
return f"解析错误：无效输入格式 - {input_str}"
except Exception as e:
return f"运行时错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 模式1测试
print("模式1：2^X")
print(pow_2_e("3"))
# 8.0
# 模式2测试
print("\n模式2：F*2^E")
print(pow_2_e("(2, 3)"))
# 16.0
print(pow_2_e("(x, y)"))
# 2.0**y*x
PowerExpand 幂扩展
幂扩展
PowerExpand(x)是power_exp提示的包装扩展
x -- 输入,符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def power_expand(input_str):
"""
将输入的数学表达式字符串展开幂表达式，对标MathStudio的powerExpand函数。
参数:
input_str (str): 要展开的数学表达式字符串。
返回:
sp.Expr 或 str: 展开后的表达式，如果输入错误则返回错误信息字符串。
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查表达式是否包含符号或是纯数值
if expr.free_symbols or expr.is_number:
# 第一步：展开指数中的加法和减法，例如 exp(x + y) -> exp(x)*exp(y)
result = sp.expand_power_exp(expr)
# 第二步：展开幂的基底，例如 (a^b)^c -> a^(b*c)，force=True强制展开不考虑假设条件
result = sp.expand_power_base(result, force=True)
else:
error = True
# 返回处理后的结果或错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 捕获异常并返回错误信息
return f"Error: {e}"
# 示例代码
if __name__ == "__main__":
# 测试案例列表，每个案例为元组（输入表达式， 期望输出）
test_cases = [
("exp(x + y)"),
# exp(x)*exp(y)
("(x**2)**(1/2)"),
# sqrt(x**2)
("(3*a)**n"),
# 3**n*a**n
]
# 遍历测试案例并验证结果
for input_expr in test_cases:
print(f"输入: {input_expr}")
output = power_expand(input_expr)
print(f"实际输出: {output}\n")
powerfit 幂拟合
多项式ij曲线拟合
p = powerfit(x,y,n) 返回次数为n的多项式 p(x) 的系数.
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是多幂拟合的次数, 正整数标量.
power1: Y = a*x^b
power2: Y = a*x^b+c
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.optimize import curve_fit
import numpy as np
def power_fit_nonlinear(input_str):
"""
使用幂函数对数据进行非线性拟合，对标MATLAB的powerfit函数。
参数:
input_str: 输入字符串，格式为"(x_data, y_data, n)"，其中n为模型类型（1或2）。
返回:
拟合的SymPy表达式或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
expression = None
maxfev = 10000 # 最大函数评估次数
# 定义幂函数模型
def power_law1(x, a, b):
"""基础幂函数: y = a*x^b"""
return a * x ** b
def power_law2(x, a, b, c):
"""带常数项的幂函数: y = a*x^b + c"""
return a * x ** b + c
if isinstance(expr, tuple):
# 解析输入参数
if len(expr) > 2:
x, y, n = expr[0], expr[1], int(expr[2])
else:
x, y = expr[0], expr[1]
n = 1 # 默认为基础幂函数模型
# 转换为矩阵并获取数据
X = sp.Matrix(x) if isinstance(x, list) else None
Y = sp.Matrix(y) if isinstance(y, list) else None
if X is None or Y is None:
return "错误: 无法转换输入数据为矩阵"
x_data = np.array(X, dtype=float).flatten()
y_data = np.array(Y, dtype=float).flatten()
# 根据模型类型进行拟合
if n == 1:
initial_guess = [1.0, 1.0]
params, _ = curve_fit(power_law1, x_data, y_data, p0=initial_guess, maxfev=maxfev)
a, b = params
expression = f"{a:.4f} * x**{b:.4f}"
elif n == 2:
initial_guess = [1.0, 1.0, 0.0]
params, _ = curve_fit(power_law2, x_data, y_data, p0=initial_guess, maxfev=maxfev)
a, b, c = params
expression = f"{a:.4f} * x**{b:.4f} + {c:.4f}"
else:
return f"错误: 不支持的模型类型 n={n}，仅支持1或2"
return sp.sympify(expression) # 直接转换为SymPy表达式
else:
return "错误: 输入格式应为 (x_data, y_data, n)"
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例1：y = 2x² 拟合
example1 = "([1, 2, 3, 4], [2, 8, 18, 32], 1)"
print("拟合结果1:", power_fit_nonlinear(example1))
# 2.0*x**2.0
# 示例2：y = 3x^1.5 + 4 拟合
example2 = "([1, 2, 3, 4], [7.0, 8.2426, 11.588, 16.0], 2)"
print("拟合结果2:", power_fit_nonlinear(example2))
# 0.3904*x**2.3044 + 6.5197
ppval 计算分段多项式
计算分段多项式
v = ppval(pp,xq) 在查询点 xq 处计算分段多项式 pp
pp — 分段多项式, 结构体
xq — 查询点, 向量
v — 查询点处的分段多项式值, 向量 | 矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.interpolate import PPoly
def polynomials_segmented(input_str):
"""
实现MATLAB ppval函数，计算分段多项式在给定点的值
参数格式：
"(breaks, coefs, pieces, order, xi)" -> 返回插值结果
或
"(breaks, coefs, pieces, order)" -> 返回符号表达式
输入要求：
- breaks: 分段区间断点，单调递增
- coefs: 多项式系数矩阵，形状(pieces, order)
- pieces: 分段数
- order: 多项式阶数
- xi: 查询点
返回：
符号表达式/数值结果 或 错误信息
"""
try:
# 解析输入字符串
expr = sp.sympify(input_str)
# 提取参数
params = []
if isinstance(expr, tuple):
params = list(expr)
else:
params = [expr]
# 参数验证
if len(params) < 4:
return "错误：至少需要breaks, coefs, pieces, order四个参数"
# 转换数据类型
breaks = np.array(params[0], dtype=np.float64)
coefs = np.array(params[1], dtype=np.float64)
pieces = int(params[2])
order = int(params[3])
xi = params[4] if len(params) > 4 else None
# 输入验证
if len(breaks) != pieces + 1:
return f"错误：breaks长度应为pieces+1 ({pieces + 1}), 实际为{len(breaks)}"
if not np.all(np.diff(breaks) > 0):
return "错误：breaks必须单调递增"
if coefs.shape != (pieces, order):
return f"错误：coefs形状应为({pieces}, {order}), 实际为{coefs.shape}"
# 构造分段多项式对象（兼容MATLAB的系数排列顺序）
pp = PPoly(coefs.T, breaks) # Scipy需要转置系数矩阵
# 符号表达式构造
t = sp.symbols('t')
piecewise_expr = []
for i in range(pieces):
# 提取当前区间的系数（注意MATLAB的系数顺序）
coeff = coefs[i, :]
# 构造多项式表达式：coeff[0]*(t-break_i)^{order-1} + ... + coeff[-1]
poly_terms = [
coeff[j] * (t - breaks[i]) ** (order - 1 - j)
for j in range(order)
]
poly = sum(poly_terms)
# 定义区间条件
if i == pieces - 1:
cond = (t >= breaks[i]) & (t <= breaks[i + 1])
else:
cond = (t >= breaks[i]) & (t < breaks[i + 1])
piecewise_expr.append((poly, cond))
# 构建完整表达式
final_expr = sp.Piecewise(*piecewise_expr)
# 根据输入类型返回结果
if xi is not None:
# 数值计算
xi = np.asarray(xi, dtype=np.float64)
return pp(xi)
else:
return final_expr
except Exception as e:
return f"错误：{str(e)}"
if __name__ == "__main__":
# ================== 测试用例 ==================
# 测试1：符号表达式
input1 = "([0, 1, 2], [[1, 2], [3, 4]], 2, 2)"
expr = polynomials_segmented(input1)
print("测试1 符号表达式:\n", expr)
# Piecewise((1.0*t + 2.0, (t < 1.0) & (t >= 0)),
# (3.0*t + 1.0, (t >= 1.0) & (t <= 2.0)))
# 测试2：数值计算
input2 = "([0, 1, 2, 3], [[1,2], [3,4], [5,6]], 3, 2, [0.5,1.5,2.5])"
print("\n测试2 插值结果:", polynomials_segmented(input2))
# [2.5 5.5 8.5]
# 测试3：复杂分段多项式
input3 = (
"([0, 1, 3, 4], "
"[[2, 3, 1, 5], [4, 1, 6, 2], [7, 3, 2, 1]], "
"3, 4, [0.5, 2.0, 3.5])"
)
print("\n测试3 插值结果:", polynomials_segmented(input3))
# [ 6.5 13. 3.625]
prctile 数据集的百分比
数据集的百分比
P = prctile(A,p)返回区间[0,100]中百分比P的输入数据A中元素的百分位数.
如果A是一个向量,那么p是一个标量.
如果A是矩阵,那么p是行向量或矩阵,其中P的行数等于长度(P).
P = prctile(A, p, 0)返回包含每列百分比的行向量.
P = median(A, p, 1)返回包含每行百分比的行向量.
A - 输入数组, 向量, 矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def percentiles_cal_array(input_str):
"""
对标 MATLAB 的 prctile 函数，计算数据的百分位数。
参数:
input_str: 输入字符串，格式为 SymPy 可解析的元组或矩阵表达式。
示例:
- 矩阵和百分位: "(Matrix([[1,2],[3,4]]), 50)"
- 指定轴: "(Matrix([[1,2],[3,4]]), 50, 1)" (MATLAB的dim=1对应axis=0)
返回:
计算结果的列表或数值，出错时返回错误信息。
"""
try:
# 修正输入表达式
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
result = None
# 检查表达式是否为元组
if isinstance(expr, tuple):
if len(expr) == 3:
M = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
# 检查是否为矩阵，百分位数是否在 0 到 100 之间，轴是否为 0 或 1
if M is not None and 0 <= int(expr[1]) <= 100 and int(expr[2]) in [0, 1]:
a = np.array(M)
p = int(expr[1])
axis = int(expr[2])
# 计算百分位数
result = np.percentile(a, p, axis=axis, method="hazen")
elif len(expr) == 2:
M = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
# 检查是否为矩阵，百分位数是否在 0 到 100 之间
if M is not None and 0 <= int(expr[1]) <= 100:
a = np.array(M)
p = int(expr[1])
# 计算百分位数
result = np.percentile(a, p, method="hazen")
return result if result is not None else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例测试
if __name__ == "__main__":
# 示例1: 计算矩阵各列的中位数 (对标MATLAB prctile(A, 50, 1))
input_str = "([0.5377,1.8339,-2.2588,0.8622,0.3188,-1.3077,-0.4336],42)"
print(percentiles_cal_array(input_str))
# -0.102544000000000
prdisc 采用银行贴现率计算收益率的折价证券价格
采用银行贴现率计算收益率的折价证券价格
Price = prdisc(Settle,Maturity,Face,Discount)返回收益率以银行贴现率报价的证券的价格（例如，美国国库券）。
Settle - 结算日期
Maturity - 到期日
Discount - 证券银行贴现率
Face - 赎回（面值）值，指定为数值。
Basis - 用于时间因素计算的天数基础
Price - 折扣证券的价格，以数值形式返回。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
from datetime import datetime
def price_discounted_security(Settle,Maturity,Face,Discount,Basis=0):
"""计算贴现证券价格，对标MATLAB的prdisc函数"""
try:
# 解析输入参数
def date_to_datetime(date_str):
"""将日期字符串转换为datetime对象"""
try:
return datetime.strptime(date_str, '%Y-%m-%d')
except ValueError:
try:
return datetime.strptime(date_str, '%Y%m%d')
except ValueError:
raise ValueError(f"日期格式错误: {date_str}，支持YYYY-MM-DD或YYYYMMDD格式")
def days_between_dates(settle, maturity, basis=0):
"""计算两个日期之间的天数，根据日计数基准调整"""
d1 = date_to_datetime(settle)
d2 = date_to_datetime(maturity)
if d2 <= d1:
raise ValueError("到期日必须晚于结算日")
# 计算实际天数差
days = (d2 - d1).days
# 根据不同的日计数基准调整（简化版）
if basis == 0: # 美国(30/360)
# 简化版30/360计算
y1, m1, d1 = d1.year, d1.month, d1.day
y2, m2, d2 = d2.year, d2.month, d2.day
# 调整31日为30日
d1 = min(d1, 30)
d2 = min(d2, 30)
days = 360 * (y2 - y1) + 30 * (m2 - m1) + (d2 - d1)
elif basis == 1: # 实际/实际
pass # 已使用实际天数
elif basis == 2: # 实际/360
pass # 已使用实际天数（分母将在后续处理）
elif basis == 3: # 实际/365
pass # 已使用实际天数（分母将在后续处理）
elif basis == 4: # 欧洲(30/360)
# 简化版欧洲30/360计算
y1, m1, d1 = d1.year, d1.month, d1.day
y2, m2, d2 = d2.year, d2.month, d2.day
# 调整31日为30日
if d1 == 31:
d1 = 30
if d2 == 31:
d2 = 30
days = 360 * (y2 - y1) + 30 * (m2 - m1) + (d2 - d1)
else:
raise ValueError(f"不支持的日计数基准: {basis}")
return days
# 提取必要参数
settle = Settle
maturity = Maturity
face = float(Face) # 默认面值100
discount = float(Discount)
basis = int(Basis) # 默认日计数基准0
# 计算天数差
days = days_between_dates(settle, maturity, basis)
# 计算价格（核心公式）
if basis == 0 or basis == 4: # 30/360
price = face * (1 - discount * days / 360)
elif basis == 1: # 实际/实际
# 简化版：假设一年365天
price = face * (1 - discount * days / 365)
elif basis == 2: # 实际/360
price = face * (1 - discount * days / 360)
elif basis == 3: # 实际/365
price = face * (1 - discount * days / 365)
else:
raise ValueError(f"不支持的日计数基准: {basis}")
return round(price, 4)
except Exception as e:
return [], [], f"计算失败: {str(e)}"
# 测试示例（对标MATLAB：Price = prdisc('2000-10-14','2001-03-17',100,0.087)）
result = price_discounted_security(Settle='2000-10-14',Maturity='2001-03-17',Face=100,Discount=0.087,Basis=0)
print(f"证券价格: {result}")
# 证券价格: 96.3025
primes 小于等于输入值的质数
小于等于输入值的质数
p = primes(n) 返回包含所有小于或等于n的质数的行向量. p与n具有相同的数据类型.
n — 输入值,标量、实整数值
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.core.sympify import SympifyError
def prime_range_list(input_str):
"""
对标MATLAB的primes函数，生成所有小于等于输入值的质数列表
参数:
input_str: 输入字符串，应能被SymPy解析为整数表达式
返回:
质数列表 (若输入有效)
错误信息 (若输入无效或无法解析)
示例:
prime_range_list("10") 返回 [2, 3, 5, 7]
prime_range_list("5.5") 返回 "输入错误：请输入非负整数，当前输入为5.5"
"""
try:
# 尝试将输入字符串解析为SymPy表达式
expr = sp.sympify(input_str)
# 检查是否为数值类型
if not expr.is_number:
return f"输入错误：'{input_str}' 不是有效数字"
# 检查是否为整数
if not expr.is_integer:
return f"输入错误：请输入非负整数，当前输入为{input_str}"
# 转换为Python整数
n = int(expr)
# 检查是否为非负数
if n < 0:
return f"输入错误：负数 {n} 无效"
# MATLAB primes(n) 当n<2时返回空数组
if n < 2:
return []
# 生成质数列表 (使用primerange包含上限需要+1)
return list(sp.primerange(2, n + 1))
# 处理表达式解析错误
except SympifyError:
return f"解析错误：'{input_str}' 不是有效表达式"
# 处理其他意外错误
except Exception as e:
return f"运行时错误：{str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试正常输入
print(prime_range_list("10"))
# [2, 3, 5, 7]
print(prime_range_list("2"))
# [2]
print(prime_range_list("1"))
# []
# 测试边界值
print(prime_range_list("0"))
# []
print(prime_range_list("1+2"))
# [2, 3]
primorial 质数阶乘
质数阶乘
primorial(n)返回所有小于或等于n的素数的乘积
n -- 输入, 标量, 整数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.core.sympify import SympifyError
def prime_product(input_str):
"""
对标MATLAB的primorial函数，计算质数阶乘
参数:
input_str: 输入字符串，支持两种模式：
- 单个整数n: 计算前n个质数的乘积 (MATLAB默认行为)
- 元组(n, False): 计算所有≤n的质数的乘积
返回:
计算结果数值 (若输入有效)
错误信息 (若输入无效或无法解析)
示例:
prime_product("5") 返回 2310 (2×3×5×7×11)
prime_product("(5, False)") 返回 30 (2×3×5)
"""
try:
expr = sp.sympify(input_str)
nth_flag = True # 默认MATLAB模式：前n个质数的乘积
n = None
# 解析输入参数
if isinstance(expr, tuple):
if len(expr) != 2:
return f"输入错误：需要(数值, 模式)格式，当前长度{len(expr)}"
# 提取数值部分
n_expr = expr[0]
if not n_expr.is_integer or n_expr < 0:
return f"数值错误：{n_expr} 应为非负整数"
# 提取模式标识
mode = expr[1]
if mode == sp.false:
nth_flag = False
elif mode == sp.true:
nth_flag = True
else:
return f"模式错误：{mode} 应为True/False"
n = int(n_expr)
elif expr.is_integer and expr >= 0:
n = int(expr)
else:
return f"输入错误：'{input_str}' 应为整数或元组"
# 执行计算
if nth_flag:
# MATLAB模式：前n个质数的乘积
if n == 0:
return 1 # 空乘积约定
return sp.primorial(n, nth=False)
else:
# 所有≤n的质数的乘积
if n < 2:
return 1
return sp.primorial(n)
except SympifyError:
return f"解析错误：'{input_str}' 不是有效表达式"
except Exception as e:
return f"运行时错误：{str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试前n个质数模式
print(prime_product("3"))
# 6
print(prime_product("5"))
# 30
print(prime_product("0"))
# 1
prmat 到期含息价格
到期含息价格
[Price,AccruInterest] = prmat(Settle,Maturity,Issue,Face,CouponRate,Yield) 返回到期支付利息的证券的价格和应计利息。
Settle - 证券结算日期
Maturity - 证券到期日
Issue - 证券发行日期
Face - 赎回价值（票面价值）
CouponRate - 息票率
Yield - 年产量，年收益
Basis - 用于时间因素计算的天数基础
Price - 证券价格，以数值形式返回
AccruInterest - 担保应计利息
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
from datetime import datetime
def price_maturity(Settle,Maturity,Issue,Face,CouponRate,Yield,Basis=0):
"""计算附息票债券（到期一次还本付息）的价格和应计利息"""
try:
# 解析输入参数
def days_360(start_date, end_date, method='US'):
"""按照30/360（美国版）规则计算两个日期之间的天数"""
d1 = start_date.day
d2 = end_date.day
m1 = start_date.month
m2 = end_date.month
y1 = start_date.year
y2 = end_date.year
# 根据美国版规则调整日
if method == 'US':
if d1 == 31:
d1 = 30
if d2 == 31 and d1 == 30:
d2 = 30
# 计算天数
total_days = (y2 - y1) * 360 + (m2 - m1) * 30 + (d2 - d1)
return total_days
def yearfrac(start_date, end_date, basis):
"""计算两个日期之间的年分数（根据日计数基准）"""
actual_days = (end_date - start_date).days
if basis == 0: # actual/actual (简化为actual/365)
return actual_days / 365.0
elif basis == 1: # 30/360
days = days_360(start_date, end_date)
return days / 360.0
elif basis == 2: # actual/360
return actual_days / 360.0
elif basis == 3: # actual/365
return actual_days / 365.0
else:
raise ValueError(f"不支持的日计数基准: {basis}")
# 将日期字符串转换为datetime对象
date_format = "%Y-%m-%d"
settle = datetime.strptime(Settle, date_format).date()
maturity = datetime.strptime(Maturity, date_format).date()
issue = datetime.strptime(Issue, date_format).date()
# 检查日期顺序
if not (issue <= settle <= maturity):
raise ValueError("日期必须满足: Issue <= Settle <= Maturity")
face = float(Face)
coupon_rate = float(CouponRate)
yield_rate = float(Yield)
basis = int(Basis)
# 计算年分数
total_period = yearfrac(issue, maturity, basis) # 发行日到到期日
accrued_period = yearfrac(issue, settle, basis) # 发行日到结算日
discount_period = yearfrac(settle, maturity, basis) # 结算日到到期日
# 计算总利息和到期总支付
total_interest = face * coupon_rate * total_period
total_payment = face + total_interest
# 计算现值（全价）
present_value = total_payment / (1 + yield_rate * discount_period)
# 计算应计利息
accrued_interest = face * coupon_rate * accrued_period
# 计算净价
clean_price = present_value - accrued_interest
# 四舍五入到两位小数
clean_price = round(clean_price, 4)
accrued_interest = round(accrued_interest, 4)
return clean_price, accrued_interest
except Exception as e:
raise RuntimeError(f"计算失败: {str(e)}")
# 测试示例
if __name__ == "__main__":
# 测试案例：票面利率等于收益率，净价应为100
price, accru = price_maturity(Settle='2002-02-07',Maturity='2002-04-13',Issue='2001-10-11',Face=100,CouponRate=0.0608,Yield=0.0608,Basis=1)
print(f"Price: {price}, Accrued Interest: {accru}")
# Price: 99.9784, Accrued Interest: 1.9591
prtbill 国库券价格
国库券价格
Price = prtbill(Settle,Maturity,Face,Discount) 返回国库券价格。
Settle - 国库券结算日期
Maturity - 国库券到期日，
Face - 国库券的赎回价值（票面价值）
Discount - 国库券贴现率
Price - 国库券价格
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
from datetime import datetime
def price_treasury_bill(Settle,Maturity,Face,Discount):
"""
计算中国短期国库券价格，对标MATLAB的prtbill函数
公式：
Price = Face * (1 - Discount * (DSM / 360))
其中：
DSM = 结算日到到期日的实际天数
注意：中国短期国库券使用实际天数/365日计数基准
参数：
Settle - 结算日
Maturity - 到期日
Face - 面值
Discount - 贴现率（小数形式）
返回：
国库券价格（保留4位小数）
"""
try:
# 将日期字符串转换为datetime对象
date_format = "%Y-%m-%d"
settle = datetime.strptime(Settle, date_format).date()
maturity = datetime.strptime(Maturity, date_format).date()
# 检查日期顺序
if settle >= maturity:
raise ValueError("结算日必须早于到期日")
face = float(Face)
discount_rate = float(Discount)
# 计算结算日到到期日的实际天数
dsm = (maturity - settle).days
# 计算国库券价格
# 如果计算美国国库券价格，这里的分母为360
price = face * (1 - discount_rate * (dsm / 365))
# 四舍五入到4位小数（MATLAB默认精度）
return round(price, 4)
except Exception as e:
raise RuntimeError(f"计算失败: {str(e)}")
# 测试示例
if __name__ == "__main__":
# 测试案例1：MATLAB示例
price = price_treasury_bill(Settle='2002-02-10',Maturity='2002-08-06',Face=1000,Discount=0.0377)
print(f"测试案例1: {price}")
# 981.7181
# 测试案例2：短期国库券（3个月）
price = price_treasury_bill(Settle='2023-01-01',Maturity='2023-04-01',Face=10000,Discount=0.0185)
print(f"测试案例2: {price}")
# 9954.3836
# 测试案例3：长期国库券（1年）
price = price_treasury_bill(Settle='2023-01-01',Maturity='2024-01-01',Face=5000,Discount=0.025)
print(f"测试案例3: {price}")
# 4875.0
# 测试案例4：高贴现率
price = price_treasury_bill(Settle='2023-03-15',Maturity='2023-06-15',Face=100000,Discount=0.05)
print(f"测试案例4: {price}")
# 98739.726
# 测试案例5：边界检查（1天）
price = price_treasury_bill(Settle='2023-12-31',Maturity='2024-01-01',Face=1000,Discount=0.01)
print(f"测试案例5: {price}")
# 999.9726
psi digamma和polygamma函数
digamma和polygamma函数
Y = psi(X) 为数组X的每个元素计算digamma函数,各元素必须为非负实数.
Y = psi(k,X) 计算X的polygamma函数,计算的值是在X处的digamma函数的k阶导数.因此,psi(0,X)是digamma函数,psi(1,X)是trigamma函数,psi(2,X)是tetragamma函数,以此类推.
X是标量,向量,矩阵,多维数组
k是导数的阶,非负整数标量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import polygamma, digamma
def psi_digamma_polygamma_func(input_str):
"""
对标 MATLAB 的 psi 函数，计算双伽玛函数 (digamma) 或多伽玛函数 (polygamma)。
参数:
input_str: 输入的字符串表达式，可以是数值、矩阵、符号表达式，或形如 "(n, x)" 的元组。
返回:
计算结果 (数值、符号表达式、矩阵) 或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 确定阶数 n 和变量 x
if isinstance(expr, tuple) and len(expr) == 2:
# 处理元组输入，如 "(n, x)"
if all(e.is_number for e in expr):
k = int(expr[0])
p = float(expr[1])
result = polygamma(k, p)
else:
result = sp.polygamma(*expr).evalf()
elif expr.is_number:
z = complex(expr)
result = digamma(z)
elif expr.free_symbols:
k = 0
p = expr
result = sp.polygamma(k, p).evalf()
else:
error = True
return result if not error else f"输入错误:{input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例 1: 计算 digamma(5)
print("digamma(5):", psi_digamma_polygamma_func("5"))
# (1.5061176684317998+0j)
# 示例 2: 计算 polygamma(1, 0.5)
print("polygamma(1, 0.5):", psi_digamma_polygamma_func("(1, 0.5)"))
# 4.93480220054468
# 示例 3: 计算 polygamma(1, x)
print("polygamma(1, x):", psi_digamma_polygamma_func("1, x"))
# polygamma(1, x)