负二项式累积分布函数
y=nbincdf(x,R,p)使用相应的成功次数,R和单次试验的成功概率计算x中每个值的负二项式cdf.
x,R和p可以是向量或矩阵,它们都具有相同的大小,也就是y的大小.
x,R或p的标量输入被扩展为与其他输入具有相同维度的常数数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, cdf
from sympy import Expr
from scipy.stats import nbinom
def negative_binomial_cdf(input_str):
"""
对标 MATLAB 的 nbincdf 函数，计算负二项分布的累积分布函数
参数:
input_str (str): 输入表达式，支持三种形式：
- 标量模式："(2, 3, 0.5)" 表示 x=2, R=3, p=0.5
- 矩阵模式："([[1,2],[3,4]], 5, 0.2)" 表示 x矩阵, R=5, p=0.2
返回:
Union[float, Matrix, str]: 计算结果或错误信息
MATLAB 参数对照：
nbincdf(x, R, p) 其中：
- x: 失败次数（标量或矩阵）
- R: 要求成功的次数（正整数）
- p: 单次试验成功概率（0 ≤ p ≤ 1）
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def nbincdf_sci(k, n, p):
"""
计算负二项分布的累积分布函数（CDF）
对标 MATLAB 的 nbincdf 函数
参数:
k (int/float/array): 成功次数
n (int): 需要的失败次数（正整数）
p (float): 每次试验成功的概率 (0 ≤ p ≤ 1)
返回:
float/array: 负二项分布的CDF值
"""
# 参数验证
if not isinstance(n, int) or n <= 0:
raise ValueError(f"参数 n 必须为正整数，当前 n={n}")
if not (0 <= p <= 1):
raise ValueError(f"概率 p 必须在 [0,1] 之间，当前 p={p}")
# 计算CDF（注意scipy的n参数表示成功次数，与MATLAB相反）
return nbinom.cdf(k, n, p)
# 计算逻辑 ======================================================
def nbincdf_sym(x_val, R, p):
# 验证 R 是正整数
if not (R.is_integer and R > 0):
return f"参数 R 必须为正整数，当前 R={R}"
# 验证 p 在概率范围内
if not (0 <= p <= 1):
return f"概率 p 必须在 [0,1] 之间，当前 p={p.evalf(4)}"
# 分布定义 ======================================================
X = NegativeBinomial("X", R, p) # 定义分布：X ~ NegativeBinomial(R, p)
"""计算单个值的 CDF，自动处理符号和数值"""
cdf_expr = cdf(X)(x_val)
# 如果是符号表达式直接返回，数值则求值
if isinstance(cdf_expr, Expr):
return cdf_expr
else:
return cdf_expr.evalf()
# 标量处理
# 参数解析 ======================================================
if isinstance(expr, tuple) and len(expr) == 3:
if all(e.is_number for e in expr):
x, R, p = float(expr[0]), int(expr[1]), float(expr[2])
result = nbincdf_sci(x, R, p)
else:
result = nbincdf_sym(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1：标量计算（数值）
print(negative_binomial_cdf("(2, 3, 0.5)"))
# 0.5
# 示例3：符号计算
x = sp.symbols('x')
print(negative_binomial_cdf(f"({x}, 2, 0.5)"))
# Piecewise((-0.25*0.5**(floor(x) + 1)*(floor(x) + 2)*hyper((1, 1.0*floor(x) + 3.0),
# (floor(x) + 2,), 0.5) + 1.0, x >= 0), (0, True))
负二项式参数估计
[parmhat，parmci]=nbinfit（data，alpha）返回MLE和100（1-alpha）百分比置信区间。默认情况下，alpha=0.05，对应于95%的置信区间。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.stats as stats
import scipy.optimize as optimize
import ast
def parse_input(input_str):
"""
解析输入字符串并转换为数值数组。
参数:
input_str (str): 输入的字符串，可以是列表、元组或SymPy矩阵的字符串形式。
返回:
numpy.ndarray: 转换后的数值数组。若解析失败则返回None。
"""
try:
# 尝试解析为Python字面量（列表、元组或数值）
data = ast.literal_eval(input_str)
if isinstance(data, (list, tuple)):
return np.array(data, dtype=float).flatten()
else:
return np.array([float(data)], dtype=float)
except:
try:
# 尝试作为SymPy表达式解析
expr = sp.sympify(input_str)
if isinstance(expr, sp.MatrixBase):
# 转换SymPy矩阵为numpy数组并展平
return np.array(expr.tolist(), dtype=float).flatten()
elif isinstance(expr, (sp.List, list, tuple)):
# 转换SymPy列表或类似结构为数组
return np.array([float(e) for e in expr], dtype=float)
else:
# 处理单个数值的情况
return np.array([float(expr)], dtype=float)
except:
return None
def negative_binomial_parameter_estimates(input_str):
"""
估计负二项分布的参数r和p，对标Matlab的nbinfit函数。
参数:
input_str (str): 输入数据的字符串形式，如"[3,5,7,2,4]"或"Matrix([1,2,3])"。
返回:
tuple or str: 成功时返回(r, p)的估计值，失败时返回错误信息。
"""
try:
data = parse_input(input_str)
if data is None:
return "错误：无法解析输入数据"
if len(data) == 0:
return "错误：输入数据为空"
# 计算矩估计作为初始猜测
mean = np.mean(data)
var = np.var(data, ddof=0) # 使用样本方差（非无偏估计）
if var > mean:
p_initial = mean / var
r_initial = mean ** 2 / (var - mean)
initial_guess = [r_initial, p_initial]
else:
# 方差不足时使用默认初始值
initial_guess = [1.0, 0.5]
# 定义负对数似然函数
def negative_binomial_nll(params):
r, p = params
# 处理参数超出范围的情况
if r <= 1e-6 or p <= 1e-6 or p >= 1 - 1e-6:
return np.inf
# 计算对数似然，避免无效值
try:
log_pmf = stats.nbinom.logpmf(data, r, p)
return -np.sum(log_pmf)
except:
return np.inf
# 参数约束：r > 0，0 < p < 1
bounds = [(1e-6, None), (1e-6, 1 - 1e-6)]
# 使用L-BFGS-B方法进行优化
res = optimize.minimize(
negative_binomial_nll,
initial_guess,
method='L-BFGS-B',
bounds=bounds
)
if res.success:
r_hat, p_hat = res.x
return (r_hat, p_hat)
else:
return f"优化失败：{res.message}"
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例输入数据（字符串形式）
test_cases = [
"[3,5,7,2,4]",
# 估计参数: r=202.9085, p=0.9798
"Matrix([3,5,7,2,4])",
# 估计参数: r=202.9085, p=0.9798
"5,3,7,8,2",
# 估计参数: r=125.0000, p=0.9615
]
for input_str in test_cases:
print(f"输入: {input_str}")
result = negative_binomial_parameter_estimates(input_str)
if isinstance(result, tuple):
r, p = result
print(f"估计参数: r={r:.4f}, p={p:.4f}\n")
else:
print(f"结果: {result}\n")
负二项式逆累积分布函数
X=nbininv(Y,R,P)返回负二项式cdf的倒数,其中包含相应的成功次数,R和单次试验中的成功概率P.
由于二项式分布是离散的,nbininv返回最小整数X,使得在X处计算的负二项cdf等于或超过Y. Y,R和P可以是向量或矩阵, 它们都具有相同的大小, 也就是X的大小.
将Y,R或P的标量输入扩展为与其他输入具有相同维度的常数数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.stats as stats
def negative_binomial_inverse_cdf(input_str):
"""
计算负二项分布的逆累积分布函数，支持标量和矩阵输入。
参数:
input_str: 输入字符串，格式为三元组 (Y, R, P)，例如 "([0.5], 5, 0.4)"
返回:
计算结果，SymPy矩阵或标量。错误时返回字符串描述。
"""
try:
# 解析输入字符串为Python对象
expr = sp.sympify(input_str)
if not isinstance(expr, (tuple, list)) or len(expr) != 3:
return f"输入错误: 需要三元组 (Y, R, P)，但得到的是 {input_str}"
def eval_element(y_val, r_val, p_val):
"""计算单个元素的逆CDF值"""
try:
y_val = float(y_val)
r_val = float(r_val)
p_val = float(p_val)
if not (0 <= y_val <= 1):
return np.nan
if not (0 < p_val < 1):
return np.nan
if r_val <= 0:
return np.nan
x = stats.nbinom.ppf(y_val, r_val, p_val)
return x
except Exception as e:
print(f"元素计算错误: {e}")
return np.nan
# 标量情况
result = eval_element(*expr)
return sp.sympify(result) if not np.isnan(result) else "无效输入值"
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 标量输入
print(negative_binomial_inverse_cdf("(0.99, 10, 0.5)"))
# 23.0000000000000
负二项式概率密度函数
Y=nbinpdf(X,R,P)使用相应的成功次数,R和单次试验的成功概率P,在X中的每个值处返回负二项式pdf.
X,R和P可以是向量或矩阵,它们都具有相同的大小,也就是Y的大小.X,R或P的标量输入被扩展为与其他输入具有相同维度的常数数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, density
def negative_binomial_pdf(input_str):
"""
计算负二项分布的概率密度函数（PDF）。
参数:
input_str: 输入字符串，格式为三元组(k, r, p)，其中：
k: 失败次数（标量、向量或矩阵）
r: 成功次数（标量、向量或矩阵）
p: 成功概率（标量、向量或矩阵）
返回:
PDF值组成的SymPy矩阵或标量。错误时返回错误信息字符串。
"""
try:
# 解析输入字符串为Python对象
expr = sp.sympify(input_str)
error = False
result = None
# 计算单个元素的PDF
def eval_element(k_val, r_val, p_val):
"""计算负二项分布PDF值"""
if r_val > 0 and (0 < p_val < 1):
X = NegativeBinomial('X', r_val, p_val)
pdf = density(X)(k_val)
return pdf.evalf() # 转换为数值形式
else:
sp.nan
if isinstance(expr, tuple) and len(expr) == 3:
result = eval_element(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{str(e)}"
# 示例用法
if __name__ == "__main__":
# 标量输入
print("标量输入:", negative_binomial_pdf("(2, 3, 0.5)"))
# 0.187500000000000
负二项式随机数
RND=nbinnd(R,P)是从负二项分布中选择的随机数矩阵,具有相应的成功数R和单次试验的成功概率P.
R和P可以是大小相同的向量或矩阵,也就是RND的大小.R或P的标量输入被扩展为与其他输入具有相同维数的常数数组.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.stats as stats
def negative_binomial_random_numbers(input_str):
"""
生成负二项分布的随机数，对标MATLAB的nbinrnd函数。
参数:
input_str: 输入字符串，格式为 (r, p) 或 (r, p, size)，
其中r为成功次数，p为成功概率，size为输出维度。
返回:
SymPy矩阵或数值。错误时返回错误信息字符串。
"""
try:
# 安全解析输入字符串
expr = sp.sympify(input_str)
error = False
result = None
def eval_random_number(r, p):
if r > 0 and (0 < p < 1):
random_number = stats.nbinom.rvs(float(r), float(p))
# 计算随机数
return random_number
else:
return sp.nan
if isinstance(expr, tuple):
# 解析参数
size = None
if len(expr) == 2:
result = eval_random_number(*expr)
elif len(expr) == 3:
n, x, size = expr
if (not size.is_integer) or size <= 0:
return "错误：size必须为正整数"
result = stats.nbinom.rvs(float(n), float(x), size=int(size))
result = sp.Matrix(result).T
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__":
# 标量输入
print("标量输出:", negative_binomial_random_numbers("(5, 0.4)"))
# 7
# 带size参数的向量输出
print("向量输出:", negative_binomial_random_numbers("(3, 0.5, 5)"))
# Matrix([[1, 11, 3, 4, 1]])
负二项均值和方差
M = nbinstat(R,P)返回负二项分布的平均值,其中R成功,单次试验成功的概率为P.
[M,V]=nbinstat(R,P)也返回分布的方差.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, E, variance
def negative_binomial_mean_variance(input_str):
"""
对标Matlab的nbinstat函数，计算负二项分布的均值和方差。
参数:
input_str: 输入表达式字符串，格式为元组 (r, p)，其中r和p可以是标量、矩阵或符号。
返回:
均值矩阵和方差矩阵。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
M, V = None, None
import numpy as np
def nbinstat_simple(r, p):
"""
计算负二项分布的均值和方差（标量输入版本）
负二项分布定义：在获得r次成功前的失败次数
参数:
r (float): 需要的成功次数（正数）
p (float): 每次试验成功的概率 (0 ≤ p ≤ 1)
返回:
tuple: (mean, variance) 两个浮点数
"""
# 参数验证
if r <= 0:
raise ValueError(f"r 必须为正数，当前值: {r}")
if not (0 <= p <= 1):
raise ValueError(f"p 必须在 [0,1] 范围内，当前值: {p}")
# 特殊情况处理
if p == 0:
return np.inf, np.inf # 不可能成功，无限次失败
if p == 1:
return 0.0, 0.0 # 立即成功，无失败
# 计算公式
q = 1 - p
mean = r * q / p
variance = r * q / (p ** 2)
return mean, variance
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
r = int(expr[0])
p = float(expr[1])
M, V = nbinstat_simple(r, p)
else:
error = True
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(negative_binomial_mean_variance("(2, 0.3)"))
# (4.666666666666667, 15.555555555555555)
二项式系数或所有组合
b = nchoosek(n,k) 返回二项式系数, 这是从 n 项中一次取 k 项的组合的数目。n 和 k 必须为非负整数
C = nchoosek(v,k) 返回一个矩阵，其中包含了从向量 v 的元素中一次取 k 个元素的所有组合。矩阵 C 有 k 列和 m!/((m–k)! k!) 行，其中 m 为 length(v)
n — 可选项的数目,非负整数标量
k — 选中项的数目,非负整数标量
b — 二项式系数,非负标量值
C — v 中的所有组合,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import itertools
def nk_binomial_coefficient(input_str):
"""
对标Matlab的nchoosek函数，实现二项式系数计算和组合生成。
参数:
input_str: 输入表达式字符串，可以是元组 (n, k) 或 (vector, k)
返回:
二项式系数值 或 组合矩阵。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
v, n = expr
if isinstance(v, list) and n.is_number:
np_v = np.array(v, dtype=float).ravel()
combinations = list(itertools.combinations(np_v, n))
# 将组合转换为矩阵
result = sp.Matrix(combinations)
else:
b_result = sp.binomial(v, n).evalf()
# 简化表达式
result = sp.combsimp(b_result)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 计算二项式系数
print(nk_binomial_coefficient("(5, 2)"))
# 10.0000000000000
print(nk_binomial_coefficient("(sympify('n'), 2)"))
# n*(n - 1)/2
# 示例2: 生成组合矩阵
print(nk_binomial_coefficient("([1, 2, 3], 2)"))
# Matrix([[1.00000000000000, 2.00000000000000],
# [1.00000000000000, 3.00000000000000],
# [2.00000000000000, 3.00000000000000]])
负对数似然
nll=negloglik（X）返回用于拟合概率分布pd的数据的负对数似然函数的值,负对数似然值越小，表示数据与分布的拟合越好。
mu — 向量X的中心值, 默认值是1.5, 用于计算每个数据点的概率密度
sigma — 向量X的标准差, 代表分布的离散程度, 默认值是0.5, 用于计算每个数据点的概率密度
normal — 正态分布, 数据对称分布在均值附近。
weibull — 威布尔分布, 描述故障，寿命等非对称数据。
X — 输入值, 向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from sympy import symbols, log, exp, Piecewise
def matrix_to_vector(matrix):
"""
将矩阵转换为向量（列表形式）。若输入不是向量，返回None。
参数:
matrix: SymPy矩阵对象
返回:
列表形式的向量，若非向量则返回None。
"""
if not isinstance(matrix, sp.Matrix):
raise ValueError("输入必须是一个 sympy.Matrix 对象")
# 获取矩阵的形状
rows, cols = matrix.shape
# 如果只有一行或一列，则将其转换为列表
if rows == 1:
# 行向量：直接转换为列表
return list(matrix.row(0))
elif cols == 1:
# 列向量：提取每一行的第一个元素并转换为列表
return [matrix[i, 0] for i in range(rows)]
else:
# 如果不是向量，返回空列表
return []
def nll_dist(input_str, distribution='normal'):
"""
计算分布的负对数似然函数，支持正态分布和Weibull分布。
参数:
input_str: 输入数据（字符串形式的向量，如 "[1,2,3]"）
distribution: 分布类型，可选 'normal' 或 'weibull'
params: 已知参数（字典，例如 {'A': 26.5, 'B': 3.27}），若为None则从数据估计
返回:
数值结果（若输入为数值且参数已知）或符号表达式
"""
try:
# 解析输入数据
expr = sp.sympify(input_str)
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is None:
return "错误：输入无法转换为矩阵"
vector = matrix_to_vector(matrix)
if vector is None:
return "错误：输入不是向量"
# 定义符号变量和参数
if distribution == 'normal':
mu, sigma = symbols('mu sigma', real=True, positive=True)
z = symbols('z', real=True)
# 正态分布PDF
pdf = (1 / (sigma * sp.sqrt(2 * sp.pi))) * exp(-(z - mu) ** 2 / (2 * sigma ** 2))
param_symbols = {'mu': mu, 'sigma': sigma}
elif distribution == 'weibull':
A, B = symbols('A B', real=True, positive=True) # A:尺度参数, B:形状参数
z = symbols('z', real=True, nonnegative=True)
# Weibull分布PDF
pdf = (B / A) * (z / A) ** (B - 1) * exp(-(z / A) ** B)
param_symbols = {'A': A, 'B': B}
else:
return "错误：不支持的分布类型"
# 负对数似然函数
neg_log_likelihood = -log(pdf)
# 总负对数似然
total_nll = sum([neg_log_likelihood.subs(z, val) for val in vector])
# 处理参数
# 从数据估计参数
if all(val.is_number for val in vector):
data = np.array([float(v) for v in vector], dtype=float)
if distribution == 'normal':
mu_est = np.mean(data)
sigma_est = np.std(data, ddof=0)
result = total_nll.subs({mu: mu_est, sigma: sigma_est}).evalf()
elif distribution == 'weibull':
# Weibull参数估计（需使用scipy等库进行更精确估计）
from scipy.stats import weibull_min
shape, loc, scale = weibull_min.fit(data, floc=0)
result = total_nll.subs({A: scale, B: shape}).evalf()
else:
result = total_nll # 返回符号表达式
return result
except Exception as e:
return f"错误：{str(e)}"
# 示例使用
if __name__ == "__main__":
# 示例1: 使用已知Weibull参数
data = "[2.5, 3.0, 5.1, 7.2, 9.5]"
params = {'A': 26.5079, 'B': 3.27193}
print(nll_dist(data, distribution='weibull'))
# 11.6157041432461
# 示例2：数值向量
input2 = "[1, 2, 3, 4]"
print(f"示例1输入: {input2}")
print("输出:", nll_dist(input2))
# 6.12204123544711
2的更高次幂的指数
P = nextpowTwo(A)返回对A中每个元素满足的最小的2的幂的指数. 按照惯例,nextpowTwo(0)返回零.
您可以使用nextpowTwo填充传递到fft的信号。当信号长度并非2次幂时，这样做可以加快FFT的运算速度
A — 输入值,由实数组成的标量,向量,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def next_pow_2_e(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
def evaluation_new_pow2(x_val):
"""核心计算函数"""
if x_val == 0:
return 0
else:
return sp.ceiling(sp.log(abs(x_val), 2))
# 处理单个数值
if expr.is_number or expr.free_symbols:
result = evaluation_new_pow2(expr)
# 无法处理的类型（如符号表达式）
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
if __name__ == "__main__":
# 示例1: 单个数值
print("示例1: next_pow_2_e('5') =", next_pow_2_e('5'))
# 3
# 示例2: 零值
print("示例2: next_pow_2_e('0') =", next_pow_2_e('0'))
# 0
# 示例3: 符号输入
print("示例3: next_pow_2_e('x') =", next_pow_2_e('x'))
# ceiling(log(Abs(x))/log(2))
尼柯尔斯图
nicholsplot(tf)
尼柯尔斯图是控制系统频域分析的重要工具，以 分贝（dB）为纵轴（增益）、度（°）为横轴（相位），绘制开环传递函数的频率响应曲线。
电机位置控制系统:
相位裕度 ≈ 45°（曲线离 -180° 线的垂直距离）。
幅值裕度 ≈ 10 dB（曲线离 0 dB 的水平距离）。
谐振峰值 ≈ 2 dB → 瞬态响应轻微振荡。
nicholsplot(10,s*(0.5s+1)*(0.1s+1))
无人机姿态控制器:
相位裕度 > 60° → 高稳定性。
带宽 ≈ 3 rad/s（0 dB 交点频率）→ 快速响应。
无谐振峰值（曲线避开 +3 dB 区域）→ 平滑跟踪。
nicholsplot(2*(s+5),s*(s+1)*(s+10))
机械臂力控制系统:
相位裕度 < 30° → 需增加补偿器。
幅值裕度 ≈ 6 dB → 临界稳定。
谐振峰值 +5 dB → 操作时存在抖动。
nicholsplot(0.8*(s+3),s^2*(s+4))
温度调节系统:
相位裕度 ≈ 70° → 强鲁棒性。
带宽 ≈ 0.5 rad/s → 缓慢调节，避免超调。
无谐振 → 温度变化平稳。
nicholsplot(1.2,s*(2s+1)*(0.2s+1))
一阶惯性系统
分析温度传感器或 RC 电路的动态响应
检查相位裕度（曲线在 -90° 附近是否远离临界点）
nicholsplot([1], [1, 2])
欠阻尼二阶系统
弹簧-质量阻尼系统或 LCR 振荡电路
观察谐振峰（曲线在 -120°~-150° 处的凸起是否接触 +3dB 等幅线）
nicholsplot([1], [1, 0.6, 1])
含微分项的系统
带速度反馈的电机位置控制
验证零点（s=-0.5s=−0.5）是否提升相位裕度（曲线向左移动）
nicholsplot([1, 0.5], [1, 2, 4])
非最小相位系统
飞机俯仰动力学或化学过程控制
右半平面零点（s=1s=1）导致高频段相位滞后（曲线向右卷绕）
nicholsplot([-1,1],[1,4,3])
三阶系统
伺服电机或机械臂关节控制
评估稳定性：曲线在 -180° 时幅值需 <0dB（避免包围临界点）
nicholsplot([2],[1,6,5,0])
叠加尼柯尔斯图
nicholsplot([tf1],[tf2]) 是在同一坐标系中绘制多个开环传递函数的尼柯尔斯曲线，用于直观比较不同系统的稳定性、鲁棒性和动态性能。
汽车悬架系统优化对比: 比较被动悬架 vs 主动悬架
主动悬架曲线（红色）相位滞后减少 25° → 过弯稳定性提升
谐振峰值从 +4dB 降至 +1.2dB → 颠簸路面振动减弱
带宽从 2 rad/s 增至 8 rad/s → 响应速度提高 3 倍
nicholsplot([0.8,(s^2+1.2s+4)*(0.02s+1)], # 被动悬架 (传统弹簧-阻尼)
[3.5*(s+5), s*(s^2+2s+6)*(0.01*s+1)]) # 主动悬架 (带液压作动器)
航空发动机燃油控制系统: 不同飞行高度下的鲁棒性验证
需增加增益调度补偿器补偿高度变化影响
nicholsplot([50,s*(s+1)*(0.1*s^2+0.5s+2)], # 海平面模式 (h=0)
[28,s*(s+0.7)*(0.15*s^2+0.3s+1.8)]) # 巡航高度模式 (h=10km)
医疗输液泵安全验证: 正常状态 vs 管路堵塞故障
故障曲线（橙色）逼近 -180°线 → 相位裕度从 65° 降至 32°
增益交点左移 → 带宽缩减 60%（触发堵塞报警阈值）
曲线进入 +6dB 危险区 → 需启动压力保护机制
nicholsplot([2.4,s*(0.5s+1)*(0.1s+1)], # 正常流动 (标称模型)
[1.2,s*(1.2s+1)*(0.15s+1)]) # 管路堵塞 (流量下降40%)
机器人关节控制器迭代: PID 控制器 vs 自适应模糊控制器
相位穿越频率从 5 rad/s → 12 rad/s → 动作更敏捷
避开 +3dB 谐振区 → 关节运动超调量从 15% 降至 3%
曲线形态更平坦 → 负载扰动敏感度降低
nicholsplot([3*(s+0.5),s^2*(0.02s+1)], # 传统PID
[4.5*(s^2+2s+1), s^2*(0.01*s^2+0.1s+1)]) # 模糊自适应控制器
光伏逆变器并网性能: 晴天 vs 阴天发电工况
阴天曲线下移 6dB → 谐波失真率增加 2.3%
相位裕度保持 >50° → 满足并网稳定标准
0dB 交点频率变化 0.5 rad/s → 需调整锁相环参数
nicholsplot([8.5,(0.05*s^2+0.3s+1)*(s+10)], # 最大功率输出 (晴天)
[3.2,(0.08*s^2+0.4s+1.2)*(s+8)]) # 降额运行 (阴天)
向量范数和矩阵范数
n = norm(v) 返回向量 v 的欧几里德范数. 此范数也称为 2-范数,向量模或欧几里德长度.
n = norm(v,p) 返回广义向量 p 范数.
n = norm(X) 返回矩阵 X 的 2-范数或最大奇异值，该值近似于 max(svd(X)).
n = norm(X,p) 返回矩阵 X 的 p-范数，其中 p 为 1,2 或 Inf
如果 p = 1, 则 n 是矩阵的最大绝对列之和.
如果 p = 2, 则 n 近似于 max(svd(X)). 此值等效于 norm(X)
如果 p = Inf, 则 n 是矩阵的最大绝对行之和.
n = norm(X,"fro") 返回矩阵或数组X的Frobenius范数
v — 输入向量, 向量
X — 输入数组, 矩阵, 数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def norm_vector_matrix(input_str):
"""
对标MATLAB的norm函数，计算向量或矩阵的范数。
参数:
input_str: 输入字符串，格式示例:
- 矩阵: "[[1,2],[3,4]]"
- 带范数类型: "([1,2,3], 2)" 或 "([[1,2],[3,4]], 'fro')"
返回:
数值结果（浮点数） 或 错误信息字符串
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
v = expr[0]
p = str(expr[1])
if isinstance(v, list) and p in ['1', '2', 'inf', 'fro']:
matrix = sp.Matrix(v)
if p == '1':
result = matrix.norm(1).evalf()
elif p == '2':
result = matrix.norm(2).evalf()
elif p == 'inf':
result = matrix.norm(sp.oo).evalf()
elif p == 'fro':
result = matrix.norm('fro').evalf()
else:
error = True
elif isinstance(expr, list):
A = sp.Matrix(expr)
result = A.norm().evalf()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 向量默认范数（2范数）
print("示例1: norm('[1,2,3]') =", norm_vector_matrix('[1,2,3]'))
# 3.74165738677394
# 示例2: 矩阵默认范数（2范数，最大奇异值）
print("示例2: norm('[[1,2],[3,4]]') =", norm_vector_matrix('[[1,2],[3,4]]'))
# 5.47722557505166
# 示例3: 矩阵Frobenius范数
print("示例3: norm('([[1,2],[3,4]], \"fro\")') =",
norm_vector_matrix('([[1,2],[3,4]], "fro")'))
# 5.47722557505166
# 示例4: 向量无穷范数
print("示例4: norm('([1,-3,2], inf)') =", norm_vector_matrix('([1,-3,2], inf)'))
# 3.00000000000000
正态累积分布函数
p = normcdf(x) 返回标准正态分布的累积分布函数 (cdf)，在 x 中的值处计算函数值。
p = normcdf(x,mu) 返回具有均值 mu 和单位标准差的正态分布的 cdf，在 x 中的值处计算函数值。
p = normcdf(x,mu,sigma) 返回具有均值 mu 和标准差 sigma 的正态分布的 cdf，在 x 中的值处计算函数值。
x — 用于计算 cdf 的值,标量值,标量值组成的数组
mu — 均值,0 (默认),标量值,标量值组成的数组
sigma — 标准差,1 (默认),非负标量值,非负标量值组成的数组
p — cdf 值,标量值,标量值组成的数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.stats import norm
import sympy as sp
from sympy.stats import Normal, cdf
def normal_distribution_cdf(input_str):
"""
对标Matlab的normcdf函数，计算正态分布的累积分布函数（CDF）。
参数:
input_str: 输入表达式字符串，可以是标量、矩阵（列表形式）、元组或符号。
返回:
正态分布的CDF值，类型可以是标量、矩阵或符号表达式。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def normal_cdf_sym(x_val, mu=0, sigma=1):
"""计算正态分布的CDF，处理数值和符号输入。"""
W = Normal("W", mu, sigma)
cdf_expr = cdf(W)(x_val)
# 仅当所有参数均为数值时求值
return cdf_expr.evalf()
if isinstance(expr, tuple):
if len(expr) == 3:
x, mu, sigma = expr
elif len(expr) == 2:
x, mu = expr
sigma = 1
else:
error = True
if all(e.is_number for e in expr):
result = norm.cdf(int(x), loc=float(mu), scale=float(sigma))
else:
result = normal_cdf_sym(x, mu, sigma)
elif expr.is_number:
x = expr
mu = 0
sigma = 1
result = norm.cdf(float(x), loc=mu, scale=sigma)
elif expr.free_symbols:
# 输入是标量或符号
result = normal_cdf_sym(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 标量输入
print(normal_distribution_cdf("0"))
# 0.5
# 示例2: 符号输入
print(normal_distribution_cdf("x"))
# 0.5*erf(sqrt(2)*x/2) + 0.5
# 数值计算
print(normal_distribution_cdf("0"))
# 0.5
print(normal_distribution_cdf("1, 0, 2"))
# 0.6914624612740131
2-范数估值
n = normest(S) 返回矩阵 S 的 2-范数估值.
本函数旨在主要用于稀疏矩阵,尽管也适用于大型满矩阵.
S — 输入矩阵,稀疏矩阵,满矩阵.
tol — 相对误差容限,1e-6 (默认),非负实数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def norm_estimate_matrix(input_str):
"""
对标MATLAB的normest函数，估计矩阵的2-范数（最大奇异值）。
参数:
input_str: 输入字符串，支持以下格式:
- 仅矩阵: "[[1,2],[3,4]]"
- 矩阵+容差: "([[1,2],[3,4]], 1e-6)"
返回:
数值估计结果（浮点数） 或 错误信息字符串
"""
try:
def eval_normest(A_np, tol=1e-6, max_iter=100):
# 如果输入是向量，直接返回其2-范数
if A_np.ndim == 1 or A_np.shape[1] == 1:
return np.linalg.norm(A_np, 2)
# 以下是矩阵的幂迭代算法
b = np.random.randn(A_np.shape[1])
b /= np.linalg.norm(b)
current_norm = 0
for _ in range(max_iter):
Ab = A_np @ b
new_b = A_np.T @ Ab
new_norm = np.linalg.norm(Ab)
new_b_norm = np.linalg.norm(new_b)
if new_b_norm == 0:
break
b = new_b / new_b_norm
if abs(new_norm - current_norm) < tol:
break
current_norm = new_norm
return current_norm
# 解析输入
expr = sp.sympify(input_str)
A_sym = None
tol = 1e-6
# 处理输入格式：可能是 (矩阵, tol) 或 仅矩阵
if isinstance(expr, tuple) and len(expr) == 2:
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
tol = float(expr[1])
else:
A_sym = sp.Matrix(expr) if isinstance(expr, list) else None
# 验证矩阵有效性
if A_sym is None:
return "错误：输入不是有效矩阵"
if A_sym.is_symbolic():
return "错误：矩阵包含非数值元素"
# 转换为NumPy数组
A_np = np.array(A_sym.tolist(), dtype=float)
# 调用核心算法
# 展平为向量（如果是列向量）
if A_np.shape[1] == 1:
A_np = A_np.flatten()
return eval_normest(A_np, tol)
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 小型稠密矩阵
print("示例1: normest([[3,0],[0,4]]) =",
norm_estimate_matrix("[[3,0],[0,4]]"))
# 3.999999438364604
# 示例2: 自定义容差
print("示例2: normest(([[1,2],[3,4]], 1e-8)) =",
norm_estimate_matrix("([[1,2],[3,4]], 1e-8)"))
# 5.4649857042187255
# 示例3: 一维输入
print("示例3: normest([1,2,3]) =",
norm_estimate_matrix("[1,2,3]"))
# 3.7416573867739413
y = normpdf(x) 返回标准正态分布的概率密度函数 (pdf)，在 x 中的值处计算函数值。
y = normpdf(x,mu) 返回具有均值 mu 和单位标准差的正态分布的 pdf，在 x 中的值处计算函数值。
y = normpdf(x,mu,sigma) 返回具有均值 mu 和标准差 sigma 的正态分布的 pdf，在 x 中的值处计算函数值。
x — 用于计算 pdf 的值,标量值 | 标量值组成的数组
mu — 均值, 0 (默认) | 标量值 | 标量值组成的数组
sigma — 标准差,1 (默认) | 正标量值 | 正标量值组成的数组
y — pdf 值,标量值 | 标量值组成的数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.stats import norm
import sympy as sp
from sympy.stats import Normal, density
def normal_distribution_pdf(input_str):
"""
对标Matlab的normcdf函数，计算正态分布的累积分布函数（CDF）。
参数:
input_str: 输入表达式字符串，可以是标量、矩阵（列表形式）、元组或符号。
返回:
正态分布的CDF值，类型可以是标量、矩阵或符号表达式。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def normal_pdf_sym(x_val, mu=0, sigma=1):
"""计算正态分布的CDF，处理数值和符号输入。"""
W = Normal("W", mu, sigma)
pdf_expr = density(W)(x_val)
# 仅当所有参数均为数值时求值
return pdf_expr.evalf()
if isinstance(expr, tuple):
if len(expr) == 3:
x, mu, sigma = expr
elif len(expr) == 2:
x, mu = expr
sigma = 1
else:
error = True
if all(e.is_number for e in expr):
result = norm.pdf(int(x), loc=float(mu), scale=float(sigma))
else:
result = normal_pdf_sym(x, mu, sigma)
elif expr.is_number:
x = expr
mu = 0
sigma = 1
result = norm.cdf(float(x), loc=mu, scale=sigma)
elif expr.free_symbols:
# 输入是标量或符号
result = normal_pdf_sym(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 标量输入
print(normal_distribution_pdf("0"))
# 0.5
# 示例2: 符号输入
print(normal_distribution_pdf("x"))
# 0.398942280401433*exp(-x**2/2)
# 数值计算
print(normal_distribution_pdf("0"))
# 0.5
print(normal_distribution_pdf("1, 0, 2"))
# 0.17603266338214976
质因子及其指数
nPrimes(n)返回质因子及对应的指数
n — 正整数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
def number_prime_factors(input_str):
"""
此函数用于将输入字符串解析为整数，并计算该整数的质因数分解结果。
参数:
input_str (str): 输入的字符串，该字符串需能被解析为正整数。
返回:
dict: 若输入合法，返回一个字典，字典的键为质因数，值为该质因数在分解中的指数。
str: 若输入不合法，返回包含错误信息的字符串。
"""
try:
# 解析输入字符串为数值
n = eval(input_str)
# 检查输入是否为整数
if isinstance(n, float):
if not n.is_integer():
return "错误：输入必须为整数"
n = int(n)
elif not isinstance(n, int):
return "错误：输入必须为整数"
# 检查输入是否为正整数
if n < 1:
return "错误：输入必须为正整数"
factors = {} # 存储质因数及其指数
# 处理所有2的因子，使n变为奇数
while n % 2 == 0:
factors[2] = factors.get(2, 0) + 1
n = n // 2
# 处理奇数因子，从3开始，每次递增2
i = 3
max_factor = int(n ** 0.5) # 初始最大因子为当前n的平方根
while i <= max_factor and n > 1:
while n % i == 0:
factors[i] = factors.get(i, 0) + 1
n = n // i
max_factor = int(n ** 0.5) # 更新最大因子为新的n的平方根
i += 2 # 跳过偶数
# 若剩余n为质数且大于2
if n > 1:
factors[n] = 1
return factors
except Exception as e:
return f"错误: {e}"
print(number_prime_factors("200"))
# {2: 2, 5: 2}
print(number_prime_factors("12345"))
# {3: 1, 5: 1, 823: 1}
实数的第n次实根
Y = nthroot(X,N)返回X的第n次实根.如果X中的元素为负数,则N必须为奇数.
X — 输入数组,标量,向量,矩阵
N — 要计算的根,标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from sympy.simplify.simplify import nthroot
def nth_root_surds(input_str):
"""
MATLAB nthroot函数的Python实现，支持矩阵广播机制
参数格式：
"(x_matrix, n_matrix)" 或 "x_matrix"
特殊行为：
- 当x是行向量(1×N)、n是列向量(M×1)时，生成M×N矩阵
- 自动处理标量与矩阵的混合运算
- 严格校验负数偶次根情况
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def nthroot_sci(x, n):
"""
计算实数的第N次实根（对标MATLAB的nthroot函数）- 标量优化版
参数:
x (float): 输入实数
n (float): 根的次数
返回:
float: 对应的实根，或NaN（当实根不存在时）
说明:
- 此版本专门优化标量输入，更简洁高效
- 对于数组输入，请使用向量化版本
"""
# 处理x=0的情况
if x == 0:
if n == 0:
return np.nan # 0^0未定义
return 0.0 if n > 0 else np.inf # 0的正数次根为0，负数次根为无穷大
# 处理n=0的情况
if n == 0:
return np.nan # 任何数的0次根未定义
# 处理x为正数的情况
if x > 0:
return x ** (1 / n)
# 处理x为负数的情况
if n % 2 == 1: # 奇数次根
return -(-x) ** (1 / n)
# 其他情况（负数开偶次根）无实数解
return np.nan
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e) for e in expr)
result = nthroot_sci(*params)
else:
result = nthroot(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 测试用例 ------------------------------------------------------------------
if __name__ == "__main__":
# 原题测试案例
test_input = "(-5,3)"
print(f"输入: {test_input}")
print("输出:")
print(nth_root_surds(test_input))
# -1.7099759466766968
test_input = "x,3"
print(f"输入: {test_input}")
print("输出:")
print(nth_root_surds(test_input))
# x**(1/3)
n次方根
nRoot(z,n)返回底数的根指数次方根的符号化结果
z — 底数
n — 根指数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def nth_root_symbolic(input_str):
"""
对输入的字符串表达式进行符号化处理，计算其n次方根。
参数:
input_str: 输入的字符串，格式应为包含两个元素的元组，如 "(z, n)"，其中z是底数，n是根指数
返回:
如果输入格式正确，返回计算得到的n次方根的符号化结果；否则返回错误信息
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为包含两个元素的元组
if isinstance(expr, tuple) and len(expr) == 2:
# 计算z的n次方根，并将结果进行数值计算
result = sp.root(*expr).evalf()
else:
# 若输入格式不符合要求，标记错误
error = True
# 如果没有错误，返回计算结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 捕获并返回计算过程中出现的错误信息
return f"错误: {e}"
# 示例输入，计算4的2次方根
input_example1 = "(4, 2)"
result1 = nth_root_symbolic(input_example1)
print(f"输入: {input_example1}, 结果: {result1}")
# 2.00000000000000
# 示例输入，使用变量进行符号化计算
x = sp.Symbol('x')
input_example3 = f"({x}, 3)"
result3 = nth_root_symbolic(input_example3)
print(f"输入: {input_example3}, 结果: {result3}")
# x**0.333333333333333
非均匀快速傅里叶变换
Y = nufft(X,t) 使用采样点 t 返回 X 的非均匀离散傅里叶变换 (NUDFT)。
如果 X 是向量，则 nufft 返回该向量的变换。
如果 X 是矩阵，则 nufft 将 X 的各列视为向量，并返回每列的变换。
如果 X 是一个多维数组，则 nufft 将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的变换。
Y = nufft(X,t,f) 使用采样点 t 计算查询点 f 的 NUDFT。要指定 f 而不指定采样点，请使用 nufft(X,[],f)。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from pynufft import NUFFT
def nufft_response(input_str):
"""
对标MATLAB的nufft函数实现非均匀快速傅里叶变换。
参数:
input_str: 输入表达式，格式为元组，其中：
- 第一个元素为数据（多维列表、SymPy矩阵、numpy数组）
- 第二个元素为非均匀采样坐标（归一化到[0, 1)范围）
- 第三个元素（可选）为输出尺寸s（默认为数据形状）
- 第四个元素（可选）为是否执行正向变换（默认为True）
返回:
- numpy数组或错误信息字符串
示例:
>>> data = np.array([1.0, 2.0, 3.0])
>>> t = np.array([[0.1], [0.5], [0.9]]) # 注意二维坐标输入
>>> nufft_response((data, t, (3,)))
array([...])
"""
try:
expr = sp.sympify(input_str, evaluate=False)
if not isinstance(expr, tuple) or len(expr) < 2:
return "错误：输入必须为元组，至少包含数据和坐标"
# 解析输入参数
data_part = expr[0]
t_part = expr[1]
s_param = expr[2] if len(expr) >= 3 else None
forward = expr[3] if len(expr) >= 4 else True
# 转换数据
# 转换数据(复数)和坐标(浮点)
data = np.array(data_part, dtype=complex) if isinstance(data_part, list) else None
t = np.array(t_part, dtype=np.float64) if isinstance(t_part, list) else None
if data is None:
return "错误：数据格式不支持"
# 转换坐标并确保二维结构
if t is None:
return "错误：坐标格式不支持"
# 确保坐标是二维数组 (N, d)
if t.ndim == 1:
t = t.reshape(-1, 1) # 转换为(N, 1)
# 归一化坐标到[0,1)并转换到[-π, π)
om = 2 * np.pi * t - np.pi
# 确定输出尺寸
if s_param is None:
s_param = data.shape
else:
# 确保s_param是元组格式
if isinstance(s_param, (list, tuple)):
s_param = tuple(s_param)
else:
s_param = (int(s_param),) # 标量转换为单元素元组
# 初始化NUFFT参数
ndim = om.shape[1] # 现在om保证是二维数组
Kd = tuple(2 * s for s in s_param) # 默认网格尺寸为2倍输出尺寸
Jd = (6,) * ndim # 邻近点数
# 创建NUFFT对象
nufft = NUFFT()
nufft.plan(om, s_param, Kd, Jd)
# 执行变换
if forward:
result = nufft.forward(data.flatten())
else:
result = nufft.adjoint(data.flatten())
return result.reshape(s_param)
except Exception as e:
return f"错误：{str(e)}"
# 正确使用的测试用例
data = [1.0, 2.0, 3.0]
t = [[0.1], [0.5], [0.9]] # 二维坐标输入
result = nufft_response(f"{data}, {t}, 3")
print("NUFFT结果:", result)
# [0.73477517+1.53991759j 5.98673999+0.00310206j 0.73477517-1.53991759j]
矩阵的零空间
Z = null(A) 返回A的零空间的标准正交基.
Z = null(A,tol) 也指定容差. 小于或等于tol的A的奇异值被视为零, 这会影响Z中的列数.
Z = null(A,"rational") 返回 A 的零空间的有理基, 它通常不是正交基. 如果A是具有小整数元素的小矩阵,则Z的元素是小整数的比率.此方法在数值上不如 null(A) 准确.
A — 输入矩阵,矩阵
tol — 奇异值容差,标量
Z — 零空间基向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def null_space_matrix(input_str):
"""
MATLAB null函数的Python实现，支持符号计算和数值计算
参数:
input_str: 输入字符串，支持格式：
- 纯矩阵（如"[[1,2],[3,6]]"）
- 矩阵+计算方式（如"(Matrix([[1,2],[3,6]]), 'rational')"）
- 矩阵+数值容差（如"(Matrix([[1,2],[3,6]]), 1e-5)"）
返回:
SymPy矩阵 或 错误信息字符串
示例:
>>> null_space_matrix("([[1,2],[3,6]], 'rational')")
Matrix([[-2], [1]])
>>> null_space_matrix("([[1,2],[3,6]], 1e-5)")
Matrix([[0.89442719], [-0.4472136]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def null_space(A, rtol=1e-5):
"""
使用SVD的数值方法计算零空间
参数:
A: numpy数组格式的矩阵
rtol: 奇异值判断阈值
返回:
零空间基向量组成的numpy数组（列向量形式）
"""
u, s, vh = np.linalg.svd(A)
# 识别接近零的奇异值（根据机器精度调整阈值）
null_mask = (s <= rtol)
# 提取对应的右奇异向量并转置为列向量
return np.compress(null_mask, vh, axis=0).T
# 处理带参数的输入，格式：(矩阵, 计算方式)
if isinstance(expr, tuple) and len(expr) == 2:
matrix_part = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
if matrix_part is None:
return f"输入错误: 无法解析矩阵部分 {input_str}"
# 符号计算模式
if str(expr[1]) == "rational":
result = matrix_part.nullspace()
# 数值计算模式
elif expr[1].is_number:
A = np.array(matrix_part.tolist(), dtype=float)
result = null_space(A, rtol=float(expr[1]))
# 处理无参数的纯矩阵输入
elif isinstance(expr, list):
A = np.array(expr, dtype=float)
result = null_space(A)
else:
error = True
# 结果格式处理
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"计算错误: {str(e)}"
# 示例测试代码
if __name__ == "__main__":
# 测试1：符号计算
print("符号计算测试:")
test_input1 = "([[1,2],[3,6]], 'rational')"
print(f"输入: {test_input1}")
print("输出:", null_space_matrix(test_input1))
# Matrix([[-2], [1]])
# 测试2：数值计算
print("\n数值计算测试:")
test_input2 = "([[1,2],[3,6]], 1e-5)"
print(f"输入: {test_input2}")
print("输出:", null_space_matrix(test_input2))
# Matrix([[-0.894427190999916], [0.447213595499958]])
提取分子和分母
[N，D]=numden(A)将A转换为有理形式,其中分子和分母是具有整数系数的相对素数多项式.
函数返回表达式有理形式的分子和分母.
如果A是符号矩阵或数字矩阵,那么N是分子的符号矩阵,D是分母的符号矩阵.N和D都是与A大小相同的矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Matrix, sympify
def numerator_denominator(input_str):
"""
MATLAB numden函数的Python实现，支持矩阵和符号表达式
参数:
input_str: 输入字符串，可以是：
- 符号表达式（如"x/y"）
- 矩阵（如"[[1/x, 2], [3, 4/y]]"）
返回:
(分子, 分母) 元组，类型可能是：
- 对于符号表达式：两个SymPy表达式
- 对于矩阵：两个SymPy矩阵
- 错误信息字符串
示例:
>>> numerator_denominator("(x**2 + 1)/y")
(x**2 + 1, y)
>>> numerator_denominator("[[1/x, 2], [3, 4/y]]")
(Matrix([
[1, 2],
[3, 4]]), Matrix([
[x, 1],
[1, y]]))
"""
try:
expr = sympify(input_str)
# 处理矩阵情况
if isinstance(expr, list):
sp_matrix = sp.Matrix(expr)
# 初始化分子分母矩阵
num_matrix = []
den_matrix = []
# 遍历矩阵每个元素
for element in sp_matrix:
num, den = sp.fraction(element)
num_matrix.append(num)
den_matrix.append(den)
# 重构矩阵保持原始形状
original_shape = sp_matrix.shape
return (
Matrix(num_matrix).reshape(*original_shape),
Matrix(den_matrix).reshape(*original_shape)
)
else:
# 对普通表达式直接分解
return sp.fraction(sp.cancel(expr))
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试符号表达式
print("符号表达式测试:")
test1 = "(x**2 + 1)/y"
print(f"输入: {test1}")
num, den = numerator_denominator(test1)
print(f"分子: {num}")
# x**2 + 1
print(f"分母: {den}\n")
# y
# 测试矩阵输入
print("矩阵输入测试:")
test2 = "[[1/x, 2], [3, 4/y]]"
print(f"输入: {test2}")
m_num, m_den = numerator_denominator(test2)
print("分子矩阵:")
print(m_num)
# Matrix([[1, 2],
# [3, 4]])
print("分母矩阵:")
print(m_den)
# Matrix([[x, 1],
# [1, y]])
纳托尔的布莱克曼-哈里斯窗
w = nuttallwin(L,sym=1) 返回一个长度为L个点的修正的纳托尔的布莱克曼-哈里斯窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import scipy.signal as signal
import sympy as sp
def nuttallwin_window(input_str):
"""
MATLAB nuttallwin函数的Python实现
参数:
params: 输入参数，支持以下格式：
- 单个整数M（窗口长度）
- 元组(M, sym_flag)：
M: 窗口长度（正整数）
sym_flag: 对称性标志（0表示周期性窗口，1表示对称窗口）
返回:
list: 生成的Nuttall窗口系数列表
或错误信息字符串
示例:
>>> nuttallwin_window(5)
[0.0, 0.2235, 0.8705, 0.8705, 0.2235]
>>> nuttallwin_window((6, 0))
[0.0003628, 0.061334, 0.52923, 1.0, 0.52923, 0.061334]
"""
try:
expr = sp.sympify(input_str)
error = False
sym = True
if isinstance(expr, tuple) and len(expr) == 2:
M = int(expr[0])
sym_flag = expr[1]
# 验证sym_flag有效性
if sym_flag not in (0, 1):
raise ValueError("sym_flag 必须为 0 或 1")
sym = bool(sym_flag)
elif expr.is_number:
# 处理单个整数输入
M = int(expr)
else:
error = True
# 验证M有效性
if not isinstance(M, int):
raise TypeError("窗口长度M必须为整数")
if M <= 0:
raise ValueError("窗口长度必须为正整数")
# 生成窗口
window = signal.windows.nuttall(M, sym=sym)
return list(window) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
# 测试基本功能
print("基本测试:")
print("nuttallwin_window(5) =>", nuttallwin_window(5))
# [0.0003628000000000381, 0.22698240000000006, 1.0, 0.22698240000000006, 0.0003628000000000381]
奈奎斯特图
nyquistplot(sys1,sys2...sysN)控制系统中用于分析系统稳定性的频域工具，它在复平面上绘制开环传递函数在频率从负无穷到正无穷变化时的轨迹。
该图以美国工程师哈里·奈奎斯特（Harry Nyquist）命名，是控制系统分析与设计中的重要工具。
直流电机位置控制系统
分析电机位置控制系统的稳定性
确定系统在不同增益下的稳定裕度
评估系统对负载变化的鲁棒性
预测超调和振荡行为
nyquistplot(5,s*(s+1)*(s+5))
温度控制系统
评估工业加热炉的温度控制性能
分析系统对温度扰动的响应特性
确定防止温度超调的最佳控制器参数
预测系统在环境温度变化时的稳定性
nyquistplot(2(s+0.5), s*(s+1)*(s+2)*(s+0.1))
飞机俯仰控制系统
分析飞行控制系统在俯仰轴上的稳定性
评估飞机在湍流中的姿态保持能力
确定自动驾驶模式下控制参数的合理范围
预测系统在极端飞行条件下的稳定性边界
nyquistplot(0.5(s+0.2)*(s+2), s^2*(s+1)*(s+5))
机械振动控制系统
评估主动减振系统的性能
分析系统在共振频率附近的稳定性
确定抑制机械振动的有效控制策略
预测系统在负载变化时的振动抑制效果
nyquistplot(0.8*(s^2+0.2s+1), s*(s^2+0.1s+4))
电源稳压系统
分析开关电源的稳定性
评估输出电压在负载突变时的恢复能力
确定防止振荡的最佳补偿网络参数
预测系统在输入电压波动时的稳定性
nyquistplot(10*(s+100), s*(s+50)*(s+200))
稳定性分析（条件稳定系统）
观察轨迹是否包围(-1,0)点。当K=5时轨迹包围(-1,0)点两次（闭环不稳定），K=1时轨迹不包围（闭环稳定）。用于设计航天器姿态控制器的增益范围。
nyquistplot([K,2K],[1,2,-3,0])
电机控制系统
确保工业机器人关节电机在负载变化时保持稳定（增益裕度>6dB）
nyquistplot([10],[1,5.5,2.5,0])
非最小相位系统（化工过程控制）
分析精馏塔温度控制系统的反常行为（右半平面零点导致反向响应）。
nyquistplot([-2,2],[1,3,2])
鲁棒性分析（飞行控制系统）
保证战斗机在气动参数摄动时（如湍流中）仍能稳定飞行（PM>40°, GM>6dB）。
nyquistplot([20,2],[1,12,20,0])
多曲线奈奎斯特图
将两个或多个传递函数的奈奎斯特图绘制在同一坐标系中（例如 nyquistplot([4s²+5s+1, 1], [3s²+2s+5, s+1/3])），这种可视化方法具有重要的工程意义和分析价值。
不同积分时间的 PI 控制器
nyquistplot([s+5,s^2+2s],[s+3,s^2+2s],[s+1,s^2+2s])
不同阻尼比的二阶系统
nyquistplot([1,s^2+0.4s+1],[1,s^2+s+1],[1,s^2+1.6s+1],[1,s^2+2s+1])
标称系统与鲁棒控制器对比
nyquistplot([s+0.5, s^3+3s^2+2s+0.5],[2s^2+s+0.2, s^3+3.5s^2+2.2s+0.6])