锯齿波或三角波
x = sawtooth(t)为时间数组t的元素生成周期为2π的锯齿波. sawtooth类似于正弦函数,但会创建峰值为-1和1的锯齿波.锯齿波定义为在2π的倍数处为-1,而在所有其他时间处以斜率为1/π 随时间呈现线性增加.
x = sawtooth(t,xmax) 生成修正三角波,其每个周期的最大值位置由xmax控制, 默认值1.
t — 时间数组,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
def sawtooth_wave(input_str):
"""
根据输入生成锯齿波或三角波，对标MATLAB的sawtooth函数。
参数:
input_str: 字符串输入，可以是以下形式：
- 时间变量，如 "t"
- 时间数组，如 "[0, 0.1, 0.2]"
- 元组形式的时间和占空比，如 "(t, 0.5)"（占空比必须在0到1之间）
返回:
tuple: (x_values, y_values) 分别为时间数组和对应的波形值。
错误时返回错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
duty = 1.0 # 默认占空比为1（锯齿波）
x_values = None
# 处理元组输入（时间和占空比）
if isinstance(expr, tuple):
if len(expr) != 2:
raise ValueError("输入元组必须包含两个元素：时间和占空比")
t_expr, duty_expr = expr
# 解析占空比
try:
duty = float(duty_expr.evalf())
if not (0 <= duty <= 1):
raise ValueError("占空比必须在0到1之间")
except:
raise ValueError("无效的占空比")
# 解析时间表达式
matrix = sp.Matrix(t_expr) if isinstance(t_expr, list) else None
if matrix is not None:
x_values = np.array(matrix).astype(float).ravel()
elif isinstance(t_expr, sp.Expr):
# 生成默认时间序列（0到2π，500个点）以匹配MATLAB的周期
var = next(iter(t_expr.free_symbols)) if t_expr.free_symbols else None
t = np.linspace(0, 2 * np.pi, 500, endpoint=False)
x_values = sp.lambdify(var, t_expr, 'numpy')(t)
else:
raise ValueError("无法解析时间表达式")
# 处理非元组输入
else:
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is not None:
x_values = np.array(matrix).astype(float).ravel()
elif isinstance(expr, sp.Expr):
# 生成默认时间序列（0到2π，500个点）
var = next(iter(expr.free_symbols)) if expr.free_symbols else None
t = np.linspace(0, 2 * np.pi, 500, endpoint=False)
x_values = sp.lambdify(var, expr, 'numpy')(t)
else:
raise ValueError("无法解析输入表达式")
# 生成波形（不再额外乘以2π）
y_values = signal.sawtooth(x_values, width=duty)
return x_values, y_values
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：基本锯齿波（width=1）
x, y = sawtooth_wave("t")
plt.figure(figsize=(10, 4))
plt.plot(x, y)
plt.title("Sawtooth Wave (duty=1)")
plt.grid(True)
plt.show()
# 示例2：复杂表达式
x, y = sawtooth_wave("t*pi/2*sin(t/5)")
plt.figure(figsize=(10, 4))
plt.plot(x, y)
plt.title("Custom Waveform: t*pi/2*sin(t/5)")
plt.grid(True)
plt.show()
舒尔分解
[U,T] = schur(A,mode) 返回酉矩阵U, 舒尔矩阵T. 并满足 A = U*T*U'.
如果mode为"real", U, T = schur(A,mode)(其中 A 为实矩阵),返回实数拟三角舒尔矩阵T. 如果mode为"complex",则返回复数三角舒尔矩阵.
如果A为复矩阵,则不管mode的值是多少,schur都返回复数舒尔形式.
t — 时间数组,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.linalg import schur as scipy_schur
def schur_decomposition(input_expr):
"""
执行舒尔分解（Schur Decomposition），对标 MATLAB 的 schur 函数。
参数:
input_expr: 输入可以是以下形式之一：
- SymPy 矩阵
- 列表或嵌套列表表示的矩阵
- 元组 (矩阵, 'real') 或 (矩阵, 'complex') 指定输出类型
返回:
若成功，返回 (Z, T)，其中：
Z: 酉矩阵（SymPy 矩阵）
T: 上三角矩阵（SymPy 矩阵）
若失败，返回错误信息字符串。
示例:
>>> schur_decomposition([[1, 2], [3, 4]])
(Matrix(...), Matrix(...))
"""
try:
# 解析输入
if isinstance(input_expr, tuple) and len(input_expr) == 2:
matrix_part, output_type = input_expr
output_type = str(output_type).lower()
if output_type not in ['real', 'complex']:
return "错误: 分解类型必须为 'real' 或 'complex'"
else:
matrix_part = input_expr
output_type = 'real' # 默认实数分解
# 转换为 SymPy 矩阵
A = sp.Matrix(matrix_part) if isinstance(matrix_part, list) else None
if A is None:
return "错误: 输入不是有效矩阵"
# 检查矩阵是否为数值矩阵
if not all(val.is_number for val in A):
return "错误: 矩阵包含符号变量，仅支持数值矩阵"
# 转换为 NumPy 数组进行舒尔分解
A_np = np.array(A, dtype=float)
# 调用 SciPy 的舒尔分解
T_np, Z_np = scipy_schur(A_np, output=output_type)
# 转换回 SymPy 矩阵
Z = sp.Matrix(Z_np)
T = sp.Matrix(T_np)
return (Z, T)
except Exception as e:
return f"错误: {e}"
# --------------------- 示例用法 ---------------------
if __name__ == "__main__":
# 示例 1: 实数分解 (默认)
matrix = [[1, 2], [3, 4]]
Z, T = schur_decomposition(matrix)
print("示例1 实数分解:")
print("Z =\n", Z)
# Matrix([[-0.824564840132394, -0.565767464968992],
# [0.565767464968992, -0.824564840132394]])
print("T =\n", T)
# Matrix([[-0.372281323269014, -1.00000000000000],
# [0, 5.37228132326901]])
print("验证 Z*T*Z.T ≈ A:\n", Z * T * Z.T)
# Matrix([[1.00000000000000, 2.00000000000000],
# [3.00000000000000, 4.00000000000000]])
# 示例 2: 复数分解
input_tuple = ([[0, -1], [1, 0]], 'complex')
Z, T = schur_decomposition(input_tuple)
print("\n示例2 复数分解:")
print("Z =\n", Z)
# Matrix([[-0.707106781186547*I, -0.531446838683691 + 0.466437839002273*I],
# [-0.707106781186547, -0.466437839002273 - 0.531446838683691*I]])
print("T =\n", T)
# Matrix([[1.0*I, -2.2677956401069e-16 + 3.09244860014577e-17*I],
# [0, 2.77555756156289e-17 - 1.0*I]])
角的正割(以弧度为单位)
Y = sec(X) 返回X的元素的正割. sec函数按元素处理数组.该函数同时接受实数和复数输入.
对于X的实数值,sec(X)返回区间 [-∞, - 1] 和 [1, ∞] 内的实数值.
对于X的复数值,sec(X)返回复数值
X — 输入角(以弧度为单位),标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def secant_of_radians(input_str):
"""
计算角的正割(以弧度为单位)，对标MATLAB的sec函数
参数:
input_str: 输入字符串，可以是:
- 数值 (如 "0", "pi/3")
- 矩阵 (如 "[[1, 0], [pi/4, pi/6]]")
- 包含符号的表达式 (如 "t", "2*x + pi/2")
返回:
计算结果 (标量/矩阵/符号表达式) 或错误信息字符串
特征:
- 支持符号计算
- 支持矩阵元素级运算
- 自动数值求值
- 完善的错误处理
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
if expr.is_number:
z = complex(expr)
result = 1 / np.cos(z)
elif expr.free_symbols:
# 符号运算后数值求值
result = sp.sec(expr).evalf()
else:
error = True
# 处理无法识别的输入类型
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 初始化美化打印
sp.init_printing(use_unicode=True)
# 示例1：基本数值计算
print("示例1：标量计算")
print("sec(0) =", secant_of_radians("0")) #
# (1+0j)
print("sec(pi/3) =", secant_of_radians("pi/3"))
# (2.0000000000000004+0j)
print("sec(pi/2) =", secant_of_radians("pi/2"))
# (1.633123935319537e+16+0j)
# 示例3：符号运算
print("\n示例3：符号表达式")
t = secant_of_radians("t")
print("sec(t) =", t)
# sec(t)
expr = secant_of_radians("2*x + pi/3")
print("复合表达式:", expr)
# sec(2*x + pi/3)
参量的正割(以度为单位)
Y = secd(X)返回X的元素的正割(以度为单位).
X — 以度为单位的角,标量值,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def secant_degree_radians(input_str):
"""
计算角度的正割（以度为单位），对标MATLAB的secd函数
参数:
input_str: 输入字符串，可以是：
- 数值角度 (如 "30", "45")
- 矩阵 (如 "[[0, 60], [90, 180]]")
- 符号表达式 (如 "theta", "2*x + 45")
返回:
计算结果 (数值/矩阵/符号表达式) 或错误信息字符串
特征:
- 自动将度数转换为弧度
- 支持符号计算
- 支持矩阵元素级运算
- 完善的错误处理
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
# 定义核心计算函数
def sec_degree_sym(x):
"""实际计算函数，处理单个元素"""
# 转换为弧度：degree * π / 180
x_rad = x * sp.pi / 180
return sp.sec(x_rad).evalf()
def sec_degree_num(x):
"""
计算以度数为单位的角度的正割值（secant）
等效于 MATLAB 的 secd 函数
参数:
x : 角度（度数），可以是标量或NumPy数组
返回:
正割值，当 cos(radians) = 0 时返回无穷大（inf）
"""
# 将角度转换为弧度
x_rad = x * np.pi / 180
# 计算余弦值
cos_vals = np.cos(x_rad)
# 计算正割 = 1 / 余弦，避免除零警告
with np.errstate(divide='ignore', invalid='ignore'):
sec_vals = np.where(np.abs(cos_vals) > np.finfo(float).eps,
1.0 / cos_vals,
np.inf)
return sec_vals
# 处理元组输入的特殊情况
if isinstance(expr, tuple):
raise ValueError("不支持元组输入，请使用数值、矩阵或符号表达式")
# 处理标量或符号表达式
if expr.is_number:
z = complex(expr)
return sec_degree_num(z)
else:
return sec_degree_sym(expr)
# 处理未知类型
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 初始化美化打印
sp.init_printing(use_unicode=True)
# 示例1：基本数值计算
print("示例1：标量计算")
print("secd(0) =", secant_degree_radians("0"))
# (1+0j)
print("secd(60) =", secant_degree_radians("60"))
# (1.9999999999999996+0j)
print("secd(90) =", secant_degree_radians("90"))
# (inf+0j)
# 示例3：符号运算
print("\n示例3：符号表达式")
theta = secant_degree_radians("theta")
print("secd(theta) =", theta)
# sec(pi*theta/180)
expr = secant_degree_radians("2*x + 45")
print("复合表达式:", expr)
# sec(pi*(x/90 + 1/4))
双曲正割
Y = sech(X) 返回X的元素的双曲正割.sech函数按元素处理数组.该函数同时接受实数和复数输入.所有的角度都以弧度表示.
X — 输入角(以弧度为单位),标量值,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def secant_hyperbolic_radians(input_str):
"""
计算双曲正割函数(sech)，对标MATLAB的sech函数
参数:
input_str: 输入字符串，可以是：
- 数值 (如 "0", "1")
- 矩阵 (如 "[[0, 1], [2, 3]]")
- 符号表达式 (如 "x", "2*y + 1")
返回:
计算结果 (数值/矩阵/符号表达式) 或错误信息字符串
特征:
- 支持符号计算
- 自动处理矩阵元素级运算
- 数值计算自动转换为浮点数
- 兼容NumPy数组输入
"""
try:
# 解析输入并替换特殊字符
expr = sp.sympify(input_str)
# 检查元组输入类型
if isinstance(expr, tuple):
raise ValueError("不支持元组输入，请使用数值、矩阵或符号表达式")
# 处理标量和符号表达式
if expr.is_number:
# 数值优化处理
z = complex(expr)
return 1 / np.cosh(z)
elif expr.free_symbols:
# 符号运算
return sp.sech(expr).evalf()
# 处理未知类型
raise TypeError("无法识别的输入类型")
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
sp.init_printing(use_unicode=True) # 美化打印
# 示例1：标量计算
print("示例1：标量计算")
print("sech(0) =", secant_hyperbolic_radians("0"))
# (1+0j)
print("sech(1) =", secant_hyperbolic_radians("1"))
# (0.6480542736638855+0j)
# 示例3：符号运算
print("\n示例3：符号表达式")
x = secant_hyperbolic_radians("x")
print("sech(x) =", x)
# sech(x)
expr = secant_hyperbolic_radians("2*y + 1")
print("复合表达式:", expr)
# sech(2*y + 1)
以列表形式返回序列
Sequence(f,x,start,end,step)依据输入字符串构建一个序列，并依据输入的不同参数对序列进行切片操作.
f -- 符号表达式
x -- 符号标量
start, end, step -- 标量, 切片参数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def squence_list(input_str):
"""
此函数的作用是依据输入字符串构建一个序列，并依据输入的不同参数对序列进行切片操作。
参数:
input_str (str): 包含序列定义和切片参数的字符串。
返回:
list 或者 str: 如果输入无误，返回切片后的序列列表；若输入有误，返回错误信息。
"""
try:
# 把输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入表达式的长度是否小于 4，若小于 4 则判定为输入错误
if len(expr) < 4:
error = True
# 若输入表达式长度为 4，构建一个完整的序列
elif len(expr) == 4:
# 构建序列，第一个参数是序列的通项公式，后面三个参数分别是变量、起始值、结束值
s = sp.sequence(expr[0], (expr[1], expr[2], expr[3]))
# 获取序列的所有元素
result = s[:]
# 若输入表达式长度为 5，构建序列并按指定步长切片
elif len(expr) == 5:
s = sp.sequence(expr[0], (expr[1], expr[2], expr[3]))
# 获取步长
step = expr[4]
# 按步长对序列进行切片
result = s[::step]
# 若输入表达式长度为 6，构建序列，按指定步长切片，再从指定起始索引开始取元素
elif len(expr) == 6:
s = sp.sequence(expr[0], (expr[1], expr[2], expr[3]))
step = expr[4]
# 获取起始索引
start_index = expr[5]
my_list = s[::step]
result = my_list[start_index:]
# 若输入表达式长度为 7，构建序列，按指定步长切片，从指定起始索引开始取指定数量的元素
elif len(expr) == 7:
s = sp.sequence(expr[0], (expr[1], expr[2], expr[3]))
step = expr[4]
start_index = expr[5]
# 获取要取的元素数量
num_elements = expr[6]
my_list = s[::step]
if start_index > 0:
result = my_list[start_index:start_index + num_elements]
else:
result = my_list[start_index - num_elements:start_index][::-1]
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示范代码
if __name__ == "__main__":
# 测试输入，构建一个从 1 到 10 的整数序列
input1 = "(i, i, 1, 10)"
print(squence_list(input1))
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# 测试输入，构建一个从 1 到 10 的整数序列，步长为 2
input2 = "(i, i, 1, 10, 2)"
print(squence_list(input2))
# [1, 3, 5, 7, 9]
# 测试输入，构建一个从 1 到 10 的整数序列，步长为 2，从索引 1 开始取元素
input3 = "(i, i, 1, 10, 2, 1)"
print(squence_list(input3))
# [3, 5, 7, 9]
# 测试输入，构建一个从 1 到 10 的整数序列，步长为 2，从索引 1 开始取 2 个元素
input4 = "(i, i, 1, 10, 2, 1, 2)"
print(squence_list(input4))
# [3, 5]
皮瑟级数
series(f,var)用f的Puiseux级数展开逼近f,直到在点var＝0处的五阶.如果不指定var,那么series将使用由f确定默认变量.
f —— 近似输入, 符号表达式,符号函数
var —— 符号变量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def series_puiseux(input_str):
"""
计算给定表达式在指定变量处的Puiseux级数展开（支持分数指数）。
参数：
input_str (str): 输入的表达式字符串，格式为"表达式, 变量"，例如"sqrt(x),x"。
如果未指定变量且表达式只有一个变量，则自动选择该变量。
返回：
sympy.Expr 或 str: 级数展开结果，若出错则返回错误信息字符串。
"""
try:
# 尝试将输入的字符串转换为 sympy 表达式，evaluate=False 表示不进行求值
expr = sp.sympify(input_str, evaluate=False)
# 初始化错误标志
error = False
# 如果表达式是一个元组且长度为 2，说明输入格式为"表达式, 变量"
if isinstance(expr, tuple) and len(expr) == 2:
# 尝试判断第一个元素是否为函数，如果不是则直接使用该元素作为表达式
f = expr[0]
# 第二个元素作为展开变量
variable = expr[1]
# 如果表达式有自由变量
elif expr.free_symbols:
# 直接使用表达式
f = expr
# 获取表达式中的自由变量
variables = f.free_symbols
# 若有自由变量，选择第一个作为展开变量
variable = tuple(variables)[0]
else:
# 若表达式没有自由变量，设置错误标志
error = True
# 默认展开的项数为 5
n = 5
# 生成级数展开并移除无穷小量
series_exp = f.series(variable, n=n).removeO()
# 如果没有错误则返回展开结果，否则返回错误信息
return series_exp if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{e}"
# 示例代码
if __name__ == "__main__":
# 示例1：展开sin(x)到5阶
input1 = "sin(x), x"
print(f"输入：{input1}")
print("输出：", series_puiseux(input1))
# -x**3/6 + x
# 示例2：展开sqrt(x)到5阶
input2 = "sqrt(x), x"
print(f"\n输入：{input2}")
print("输出：", series_puiseux(input2))
# sqrt(x)
# 示例3：多变量未指定
input3 = "x*y + y"
print(f"\n输入：{input3}")
print("输出：", series_puiseux(input3))
# x*y + y
移动数组维度
B = shiftdim(A,n) 将数组 A 的维度移动 n 个位置。当 n 为正整数时，shiftdim 向左移动维度；当 n 为负整数时，向右移动维度。
A — 输入数组, 向量 | 矩阵
n — 位置的数量, 整数
B — 输出数组, 向量 | 矩阵
var —— 符号变量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def shift_arr_dim(input_str):
"""
实现类似 MATLAB shiftdim 函数的功能
:param arr: 输入的数组，可以是 sympy 矩阵或 numpy 数组
:param n: 移位的位数，正数表示向左移位，负数表示向右移位
:return: 移位后的数组
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
if isinstance(expr, tuple) and isinstance(expr[0], list) and expr[1].is_integer:
arr = np.array(expr[0], dtype=object)
n = int(expr[1])
# 处理 n 为正数的情况（循环左移）
if n > 0:
ndim = arr.ndim
if ndim == 0:
pass # 标量无需处理
else:
n_shift = n % ndim # 处理 n 大于维度数的情况
if n_shift != 0:
new_order = list(range(n_shift, ndim)) + list(range(n_shift))
arr = np.transpose(arr, new_order)
# 处理 n 为负数的情况（添加前导单一维度）
elif n < 0:
for _ in range(-n):
arr = np.expand_dims(arr, axis=0)
else:
error = True
if not error:
if arr.ndim == 2:
# 如果输入是 sympy 矩阵且结果为二维，将结果转换回 sympy 矩阵
result = sp.Matrix(arr.tolist())
else:
result = arr
return result
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例使用
# 创建一个符号数组
# 向左移位 1 位
shifted_left = shift_arr_dim("[[x, y],[z,x+y]],1")
print("\n向左移位 1 位后的数组:")
print(shifted_left)
# Matrix([[x, z],
# [y, x + y]])
# 向左移位 1 位
shifted_left = shift_arr_dim("[[3, 4],[5,6]],1")
print("\n向左移位 1 位后的数组:")
print(shifted_left)
# Matrix([[3, 5],
# [4, 6]])
# 向右移位 1 位
shifted_right = shift_arr_dim("[[x, y],[z,x+y]],-1")
print("\n向右移位 1 位后的数组:")
print(shifted_right)
# [[[x y]
# [z x + y]]]
对数组执行符号函数
Y=sign(x)返回一个与x大小相同的数组Y,其中Y的每个元素为:
如果x的对应元素大于0,则为1.
如果x的相应元素等于0,则为0.
如果x的对应元素小于0,则为1.
如果x是复数, 返回x/abs(x).
a — 输入数组,标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def sign_function(input_str):
"""
符号函数计算，对标 MATLAB 的 sign 函数
定义:
对于实数 x:
sign(x) = 1 if x > 0
0 if x = 0
-1 if x < 0
对于复数 z = a + bi:
sign(z) = z / |z| (当 z ≠ 0 时)
0 (当 z = 0 时)
参数:
input_str (str): 输入表达式，可以是：
- 数值 (如 "2.5", "-3")
- 复数 (如 "1+2j")
- 符号表达式 (如 "x", "y**2")
- 矩阵 (如 "[[1, -2], [0, 3+4j]]")
返回:
[类型] 数值/符号表达式/矩阵 或 错误信息字符串
示例:
>>> sign_function("5")
1.0
>>> sign_function("-x")
-sign(x)
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 情况2：处理数值和符号输入
if expr.is_number or expr.free_symbols:
# 数值类型转换
# 符号表达式保持符号形式
return sp.sign(expr).evalf()
# 其他无法处理的情况
return f"输入错误: 不支持的类型 {type(expr)}"
except Exception as e:
return f"计算错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1：实数计算
print("示例1 实数计算:")
print(sign_function("5"))
# 1.00000000000000
print(sign_function("-3.2"))
# -1.00000000000000
print(sign_function("0"))
# 0
# 示例2：复数计算
print("\n示例2 复数计算:")
print(sign_function("3+4j"))
# 0.6 + 0.8*I
print(sign_function("0+0j"))
# 0
# 示例3：符号运算
print("\n示例3 符号运算:")
print(sign_function("x"))
# sign(x)
print(sign_function("-y**2"))
# -sign(y**2)
双曲正弦积分
如果x是浮点数，则sinhint(x)返回浮点结果.
Y = sinhint(x)返回算术表达式
x是标量,符号变量,表达式,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.special import shichi
def shi_hyper_sine_integral(input_str):
"""
计算双曲正弦积分函数 Shi(z)，对标 MATLAB 的 sinhint 函数
定义:
Shi(z) = ∫₀ᶻ (sinh(t)/t) dt
参数:
input_str (str): 输入表达式，可以是：
- 数值 (如 "0", "3.14")
- 符号表达式 (如 "x", "2*t+1")
- 矩阵 (如 "[[0, 1], [2, 3]]")
返回:
[类型] 数值/符号表达式/矩阵 或 错误信息字符串
示例:
>>> shi_hyper_sine_integral("0")
0.0
>>> shi_hyper_sine_integral("x")
Shi(x)
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 处理数值和符号输入
if expr.is_number:
z = complex(expr)
sin_arr, cos_arr = shichi(z)
return sin_arr
elif expr.free_symbols:
return sp.Shi(expr).evalf()
return f"输入错误: 无法处理类型 {type(expr)}"
except Exception as e:
return f"计算错误: {str(e)}"
if __name__ == "__main__":
# 示例 1：数值输入
print("示例1 数值计算:")
print(shi_hyper_sine_integral("0"))
# 0j
print(shi_hyper_sine_integral("1.0"))
# (1.0572508753757286+0j)
# 示例 2：符号计算
print("\n示例2 符号运算:")
print(shi_hyper_sine_integral("x"))
# Shi(x)
print(shi_hyper_sine_integral("2*t + 1"))
# Shi(2*t + 1)
正弦积分
sinint（X）返回X的正弦积分函数.
根据其参数,sinint返回浮点或精确的符号结果.
X是指定为符号数,变量,表达式或函数,或指定为符号,变量,公式或函数的矢量或矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import sici
def si_sine_integral(input_str):
"""
计算正弦积分函数 Si(z)，对标 MATLAB 的 sinint 函数
定义:
Si(z) = ∫₀ᶻ (sin(t)/t) dt
参数:
input_str (str): 输入表达式，可以是：
- 数值 (如 "0", "3.14")
- 符号表达式 (如 "x", "2*t+1")
- 矩阵 (如 "[[0, 1], [2, 3]]")
返回:
[类型] 数值/符号表达式/矩阵 或 错误信息字符串
示例:
>>> si_sine_integral("0")
0.0
>>> si_sine_integral("x")
Si(x)
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
result = None
error = False
# 处理元组输入
if isinstance(expr, tuple):
error = True
# 处理数值和符号输入
if expr.is_number:
z = complex(expr)
result = sici(z)[0]
elif expr.free_symbols:
result = sp.Si(expr.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：数值输入
print("示例1 数值计算:")
print(si_sine_integral("0"))
# 0j
print(si_sine_integral("1.0"))
# (0.946083070367183-0j)
# 示例 2：符号计算
print("\n示例2 符号运算:")
print(si_sine_integral("x"))
# Si(x)
print(si_sine_integral("2*t + 1"))
# Si(2.0*t + 1.0)
代数简化
S＝simplify(expr)执行expr的代数简化. 如果expr是符号向量或矩阵,则此函数简化了expr的每个元素.
expr — 输入表达式, 符号表达式, 符号函数, 符号向量, 符号矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def simplify_function(input_str):
"""
简化输入的数学表达式或矩阵中的每个元素，对标Matlab的simplify函数。
参数:
input_str (str): 输入的数学表达式或矩阵的字符串表示。
返回:
SymPy对象或错误信息字符串: 简化后的表达式或矩阵，若出错则返回错误信息。
示例:
>>> simplify_function("x + x + 2*y - y")
2*x + y
>>> simplify_function("[[x + x, 2*y], [sin(pi/2), sqrt(4)]]")
Matrix([[2*x, 2*y], [1, 2]])
>>> simplify_function("3 + 4")
7
>>> simplify_function("(1, 2)")
'输入错误: (1, 2)'
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def eval_expression(expr_item):
"""内部函数，用于简化单个表达式元素。"""
return sp.simplify(expr_item)
# 检查输入是否为元组（视为无效输入）
if isinstance(expr, tuple):
error = True
else:
# 处理普通标量或符号表达式
result = eval_expression(expr)
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 简化符号表达式
print(simplify_function("x + x + y - y/2"))
# 2*x + y/2
# 简化数值表达式
print(simplify_function("sqrt(9) + cos(pi)")) # 2
简化符号有理表达式
simplifyFraction(expr)简化有理表达式expr,使得分子和分母没有共同的除数.
simplifyFraction(expr,'Expand',true)将所得简化分数的分子和分母展开为无需分解的多项式.
expr — 输入数字,向量,矩阵,数组,符号数,符号变量,符号数组,符号函数,符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def simplify_fraction(input_str):
"""
简化符号有理表达式，对标Matlab的simplifyFraction函数，支持矩阵输入。
参数:
input_str (str): 输入的有理表达式或矩阵的字符串表示。
返回:
SymPy对象或str: 简化后的表达式或矩阵，错误时返回错误信息。
示例:
>>> simplify_fraction("(x^2 - 1)/(x - 1)")
x + 1
>>> simplify_fraction("[[(x^2 - 1)/(x - 1), 6/3], [2*y/(y+y), (a+b)^2/(a+b)]]")
Matrix([[x + 1, 2], [1, a + b]])
>>> simplify_fraction("4/2 + 6/3")
4
>>> simplify_fraction("invalid expr")
错误:...
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
return sp.cancel(expr)
except Exception as e:
return f"错误:{e}"
print(simplify_fraction("((y+1)^2*(x^2-1))/((x+1)*(x-1)^2)"))
# (y**2 + 2*y + 1)/(x - 1)
以弧度为单位的参量的正弦.
Y = sin(X)返回X的每个元素的余弦.sin函数按元素处理数组.该函数同时接受实数和复数输入.
对于X的实数值,sin(X)返回区间[-1, 1]内的实数值.
对于X的复数值,sin(X)返回复数值.
X是标量，向量，数组，矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def sine_trig_function(input_str):
"""
对标 MATLAB 的 sin 函数，计算以弧度为单位的正弦值
参数:
input_str: 输入的数学表达式字符串，可以是数值、符号、矩阵等
返回:
计算结果（数值、符号表达式或矩阵），若输入不合法则返回错误信息
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 处理元组输入（视为非法）
if isinstance(expr, tuple):
error = True
elif expr.is_number or expr.free_symbols:
result = sp.sin(expr).evalf()
# 其他无法处理的类型
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例 1: 数值输入
print("示例 1:", sine_trig_function("0"))
# 0
print("示例 1:", sine_trig_function("pi/2"))
# 1.00000000000000
# 示例 2: 符号输入
x = sp.symbols('x')
print("示例 2:", sine_trig_function("x"))
# sin(x)
print("示例 2:", sine_trig_function("x + pi"))
# -sin(x)
辛格函数
y = sinc(x) 返回数组 y，其元素是输入 x 的元素的 sinc。输出 y 与 x 的大小相同.
x是标量值,向量,矩阵或多维数组.
y是标量值,向量，矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def sinc_sampling_function(input_str):
"""
计算输入表达式的 sinc 采样函数值，对标 MATLAB 的 sinc 函数。
MATLAB 定义：sinc(x) = sin(πx)/(πx)，且在 x=0 时值为 1
参数:
input_str (str): 输入表达式字符串，可以是：
- 数值（如 "0"）
- 符号表达式（如 "x"）
- 矩阵（如 "[[0, 1], [2, 3]]"）
返回:
sympy 表达式/矩阵/浮点数 或 str: 计算结果或错误信息
"""
try:
# 将输入字符串转换为 SymPy 对象
expr = sp.sympify(input_str)
def sinc_numpy_style(x):
if x == 0:
return 1
else:
return sp.sin(sp.pi * x) / (sp.pi * x)
# 情况 2：处理数值输入
if expr.is_number:
z = complex(expr)
return np.sinc(z)
# 情况 3：处理符号表达式
elif expr.free_symbols:
return sinc_numpy_style(expr).evalf()
# 其他无法处理的情况
else:
return f"输入错误: 不支持的类型 {type(expr)}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例 1：数值输入
print("示例 1:")
print(sinc_sampling_function("0"))
# (1+0j)
print(sinc_sampling_function("1.0"))
# (3.8981718325193755e-17-0j)
# 示例 2：符号输入
print("\n示例 2:")
x = sp.symbols('x')
print(sinc_sampling_function("x"))
# 0.318309886183791*sin(pi*x)/x
print(sinc_sampling_function("2*x + 1"))
# 0.318309886183791*sin(pi*(2*x + 1))/(2.0*x + 1.0)
以度为单位的参量的正弦.
Y = sind(X)返回X的每个元素的余弦.sind函数按元素处理数组.该函数同时接受实数和复数输入.
X是标量，向量，数组，矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def sind_sine_degree(input_str):
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 定义度数转弧度的正弦计算函数
def sine_deg(x_val):
return sp.sin(x_val * sp.pi / 180).evalf()
# 处理数值或符号表达式
if expr.is_number:
z = complex(expr)
result = np.sin(z * np.pi / 180)
elif expr.free_symbols:
result = sine_deg(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:", sind_sine_degree("30"))
# (0.49999999999999994+0j)
# 示例 2: 符号输入
x = sp.symbols('x')
print("示例 2:", sind_sine_degree("x"))
# sin(pi*x/180)
正弦拟合
p = sinfit(x,y,n) 返回次数为n的正弦拟合 p(x) 的系数
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是正弦拟合的次数, 正整数标量.
Y = a1*sin(b1*x+c1)
Y = a1*sin(b1*x+c1)+a2*sin(b2*x+c2)
Y = a1*sin(b1*x+c1)+...+a3*sin(b3*x+c3)
Y = a1*sin(b1*x+c1)+...+a8*sin(b8*x+c8)
以此类推，最大到 sin8
# 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 sin_fit_nonlinear(input_str):
"""
使用非线性最小二乘法进行正弦曲线拟合，支持1到8阶正弦项。
参数:
x_data (array_like): 自变量数据数组
y_data (array_like): 因变量数据数组
n (int): 正弦项的数量（1到8）
maxfev (int): 最大函数评估次数
返回:
sympy.Expr: 拟合的正弦表达式
str: 错误信息（如果发生错误）
示例:
x = np.linspace(0, 10, 100)
y = 2.5 * np.sin(3 * x + 1) + 0.5 * np.sin(5 * x - 2)
expr = sin_fit_nonlinear(x, y, n=2)
print(expr) # 输出类似 2.5000*sin(3.0000*x + 1.0000) + 0.5000*sin(5.0000*x - 2.0000)
"""
try:
expr = sp.sympify(input_str)
error = False
maxfev = 100000
# 定义正旋拟合函数，从一阶到八阶
if isinstance(expr, tuple):
if len(expr) == 3:
x, y, n = expr[0], expr[1], int(expr[2])
elif len(expr) == 2:
x, y = expr
n = 1
else:
error = True
if isinstance(x, list):
x_data = np.array(x, dtype=float).ravel()
if isinstance(y, list):
y_data = np.array(y, dtype=float).ravel()
# 输入验证
if n < 1 or n > 8:
return None, "阶数n必须在1到8之间"
if len(x_data) != len(y_data):
return None, "x_data和y_data长度不一致"
else:
error = True
# 定义各阶正弦函数
func_dict = {
1: lambda x, *p: p[0] * np.sin(p[1] * x + p[2]),
2: lambda x, *p: sum(p[i * 3] * np.sin(p[i * 3 + 1] * x + p[i * 3 + 2]) for i in range(2)),
# ... 类似定义到n=8
}
# 生成初始猜测
a_guess = (np.max(y_data) - np.min(y_data)) / 2
b_guess = 2 * np.pi / (np.max(x_data) - np.min(x_data))
initial_guess = [a_guess, b_guess, 0] * n
params, _ = curve_fit(func_dict[n], x_data, y_data,
p0=initial_guess, maxfev=maxfev,
bounds=(0, [np.inf, np.inf, 2 * np.pi] * n))
# 构建表达式
terms = []
for i in range(n):
a, b, c = params[i * 3: (i + 1) * 3]
terms.append(f"{a:.4f}*sin({b:.4f}*x + {c:.4f})")
return sp.sympify("+".join(terms))
except Exception as e:
return f"错误:{e}"
# 生成测试数据
x = [1, 2, 3, 4, 5, 6, 7, 8]
y = [9, 10, 11, 12, 13, 14, 15, 16]
# 进行拟合
expr = sin_fit_nonlinear(f"{x},{y},2")
if expr:
print("拟合表达式:", expr)
# 5.1845*sin(1.0957*x)
else:
print("错误:")
双曲正弦
Y = sinh(X)返回X的每个元素的双曲正弦.sinh函数按元素处理数组.该函数同时接受实数和复数输入.
所有的角度都以弧度表示.
X是标量,向量,数组,矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def sin_hyperbolic_sine(input_str):
"""
计算输入的双曲正弦值（对标MATLAB sinh函数）
参数:
expr: 输入数据，支持以下类型：
- 数值（int/float）
- SymPy符号表达式
- SymPy矩阵
- 嵌套列表（自动转换为SymPy矩阵）
返回:
SymPy对象（数值/矩阵/表达式）或错误信息字符串
示例:
1. 标量计算
>>> sin_hyperbolic_sine(0)
0
2. 符号计算
>>> x = sp.Symbol('x')
>>> sin_hyperbolic_sine(x)
sinh(x)
3. 矩阵计算
>>> sin_hyperbolic_sine([[0, 1], [2, 3]])
Matrix([
[0, sinh(1)],
[sinh(2), sinh(3)]])
4. 错误处理
>>> sin_hyperbolic_sine("invalid")
"错误: 不支持的类型: "
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理标量或符号表达式输入
if expr.is_number:
z = complex(expr)
result = np.sinh(z)
elif expr.free_symbols:
result = sp.sinh(expr).evalf()
# 其他类型错误
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 以下是测试用例
if __name__ == "__main__":
# 标量测试
print("标量测试:", sin_hyperbolic_sine(0))
# 0j
# 符号测试
x = sp.Symbol('x')
print("符号测试:", sin_hyperbolic_sine(x))
# sinh(x)
准确计算sin(x*pi)
Y = sinpi(X).该函数同时接受实数和复数输入.
X是标量,向量,数组,矩阵.
对于奇数,sinpi(n/2) 为+1或-1.
对于整数,sinpi(n) 为零
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def sin_pi_sine(input_str):
"""
准确计算 sin(X*π) 的符号表达式和数值计算
对标MATLAB sinpi函数特性：
1. 精确处理符号表达式 (如 sin(pi*x))
2. 数值计算时不直接展开π，避免浮点误差
3. 支持矩阵元素级计算
参数:
input_str (str): 输入表达式，可以是：
- 数值 (如 "0.5")
- 符号表达式 (如 "x")
- 矩阵 (如 "[[1, 2], [x, 1/2]]")
返回:
float: 数值计算结果
MutableDenseMatrix: 矩阵计算结果
sp.Expr: 符号表达式结果
str: 错误信息
示例:
>>> sin_pi_sine("0.5")
1.0
>>> sin_pi_sine("x")
sin(pi*x)
>>> sin_pi_sine("[[1/2, 0], [x, 1/6]]")
Matrix([[1.0, 0], [sin(pi*x), 0.5]])
"""
try:
# 符号化输入（自动识别矩阵和表达式）
expr = sp.sympify(input_str)
# 处理标量数值/符号
if expr.is_number:
z = complex(expr)
result = np.sin(np.pi * z)
elif expr.free_symbols:
result = sp.sin(sp.pi * expr)
return result
except Exception as e:
return f"计算错误: {str(e)}"
# ----------------------------
# 测试案例
# ----------------------------
if __name__ == "__main__":
# 案例1: 数值计算
print(sin_pi_sine("0.5"))
# (1+0j)
print(sin_pi_sine("1/6"))
# (0.49999999999999994+0j)
# 案例2: 符号计算
print(sin_pi_sine("x"))
# sin(pi*x)
# 案例4: 复合表达式
print(sin_pi_sine("x + 1/2"))
# sin(pi*(x + 1/2))
求解优化问题或方程问题
sol = solve(prob) 求解优化问题或方程问题prob
prob — 优化问题或方程问题, 线性或非线性方程组,包括等式方程及不等式方程
sol — 解
将球面坐标转换为笛卡尔坐标
[x,y,z] = sph2cart(azimuth,elevation,r) 将球面坐标数组 azimuth,elevation和r的对应元素转换为笛卡尔坐标,即 xyz 坐标
azimuth — 方位角,标量,向量,矩阵
elevation — 仰角,标量,向量,矩阵
r — 半径,标量,向量,矩阵
x,y,z — 笛卡尔坐标,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def spherical_to_cartesian(input_str):
"""
球坐标转笛卡尔坐标 (对标 MATLAB sph2cart 函数)
参数:
input_str: 字符串形式的输入，支持三种格式：
1. 单独方位角 (az)：返回 (x, y) 平面坐标
2. 方位角 (az) 和仰角 (el)：假定 r=1
3. 方位角 (az), 仰角 (el), 半径 (r)
返回:
(x, y, z) 元组，元素类型与输入保持一致（标量/矩阵）
实现特性:
- 支持符号计算
- 自动广播标量到矩阵维度
- 输入验证和错误处理
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str)
error = False
result = None
# 定义转换函数
def compute_coords(az_val, el_val, r_val):
x = r_val * sp.cos(el_val) * sp.cos(az_val)
y = r_val * sp.cos(el_val) * sp.sin(az_val)
z = r_val * sp.sin(el_val)
return x, y, z
def spherical_to_cartesian(az, el, r):
"""
将球面坐标转换为笛卡尔坐标，支持标量、向量和多维数组输入。
参数:
az: 方位角（弧度），可以是标量或NumPy数组
el: 仰角（弧度），可以是标量或NumPy数组
r: 半径，可以是标量或NumPy数组
返回:
x, y, z: 笛卡尔坐标，与输入具有相同的维度
"""
# 将输入转换为NumPy数组以支持广播
az = np.asarray(az)
el = np.asarray(el)
r = np.asarray(r)
# 执行坐标转换
x = r * np.cos(el) * np.cos(az)
y = r * np.cos(el) * np.sin(az)
z = r * np.sin(el)
return x, y, z
if isinstance(expr, tuple) and len(expr) == 3:
if all(e.is_number for e in expr):
params = tuple(float(e) for e in expr)
result = spherical_to_cartesian(*params)
elif any(e.free_symbols for e in expr):
result = compute_coords(*expr)
else:
error = True
return result if not error else f"输入错误:{input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1: 标量输入
print("示例1 输入: (pi/4, pi/6, 2)")
print("输出:", spherical_to_cartesian("(pi/4, pi/6, 2)"))
# (1.2247448713915892, 1.224744871391589, 1.0)
# 示例2: 符号输入
a, e = sp.symbols('a e')
print("\n示例2 输入: (a, e, 1)")
print("输出:", spherical_to_cartesian(f"({a}, {e}, 1)"))
# (cos(a)*cos(e), sin(a)*cos(e), sin(e))
球面坐标3D绘图.
球面坐标θ和φ的半径绘制函数.
SphericalPlot3D(L,r1,r2,options)
L是具有两个变量的过程或表达式，或三个此类过程或表达式的列表
r1,r2是表达式变量的范围.
1. 地球形状近似（轴对称椭球）:
描述地球的 两极扁平化（赤道半径 > 极半径）
赤道（θ=π/2）膨胀约 21.3 km
SphericalPlot3D(6371*(1+0.003*(1.5cos(theta)^2-0.5)), theta=[0,@pi])
2. 声学指向性模式（扬声器设计）:
优化扬声器阵列的声波覆盖
抑制特定方向干扰（如舞台反馈啸叫）
输出特征:主瓣沿 z 轴 + 4个侧瓣
SphericalPlot3D(abs(sinc(3*sin(theta)*cos(phi))), theta=[0,@pi], phi=[0,2@pi])
3. 等离子体湍流（核聚变研究）:
模拟托卡马克装置中磁流体不稳定性
预测等离子体破裂位置
输出特征:螺旋波纹结构（类似扭曲模）
SphericalPlot3D(1+0.1*(sin(theta+2phi)+sin(2theta+2phi)+sin(3theta+2phi)),theta=[0,@pi], phi=[0,2@pi])
4. 天线辐射方向图（通信工程）:
设计 5G 基站天线覆盖范围
优化卫星通信链路
输出特征: 轮胎状辐射模式（最大增益在赤道面）
SphericalPlot3D(abs(cos((@pi/2)*cos(theta))/sin(theta)),theta=[0,@pi])
5. 液滴振动模式（轴对称基础模态）:
赤道膨胀 + 两极扁平化的振荡
SphericalPlot3D(1+0.2*(1.5*cos(theta)^2-0.5),theta=[0,@pi])
6. 核物质分布（球形原子核）:
沿旋转轴压缩的椭球形状
SphericalPlot3D(1.2*(1-0.05*(1.5*cos(theta)^2-0.5)),theta=[0,@pi])
多曲面球面坐标3D图
内球提供基准形状,外曲面突出相对变形,清晰显示表面与内部的对应
1. 行星大气环流模型:
内球：地球固体表面
外曲面：大气温度异常分布
5个经向波动 + 10个纬向波动表示大气环流模式
角度限制聚焦北半球中高纬度区域
SphericalPlot3D(6371,6371+8*sin(5phi)*sin(10theta), theta=[0,0.7@pi],phi=[0,1.8@pi])
2. 核聚变装置等离子体约束
内球：真空室边界
外曲面：等离子体密度分布
sin(5ϕ)sin(10θ) 模拟磁流体不稳定性
角度限制避免极区复杂物理现象
SphericalPlot3D(1,1.5+0.2*sin(@phi)*sin(10theta), theta=[0,0.7@pi],phi=[0,1.8@pi])
3. 脑皮层神经活动映射:
内球：基准脑皮层表面
外曲面：神经活动强度
5×10 振荡模式模拟不同功能区域
角度限制对应可观测的脑区
SphericalPlot3D(0.9,1+0.15*sin(5phi)*sin(10theta), theta=[0.2@pi,0.7@pi],phi=[0.1@pi,1.8@pi])
4. 多波束天线方向图
内球：参考单位增益
外曲面：实际辐射强度
5个方位角波束 + 10个仰角波束
φ范围限制表示天线机械旋转范围
SphericalPlot3D(1,1.5+0.3*sin(5phi)*sin(10theta), theta=[0,0.7@pi],phi=[0,1.8@pi])
5. 地幔对流模型
内球：地核-地幔边界
外曲面：地幔对流速度
5个经向对流单元 + 10个纬向结构
θ限制表示上地幔区域（深0-670km）
SphericalPlot3D(0.9,1+0.25*sin(5phi)*sin(10theta), theta=[0,0.7@pi],phi=[0,1.8@pi])
6. 蛋白质分子表面
内球：分子核心
外曲面：范德华表面
5×10 波动模拟结合口袋结构
角度限制聚焦功能活性区域
SphericalPlot3D(0.8,1+0.4*sin(5phi)*sin(10theta), theta=[0.1,0.6@pi],phi=[0.2,1.6@pi])
SpotPlot(list1, list2)
通过笛卡尔坐标系展示两个变量之间关系的图表。每个数据点由横坐标（X轴）和纵坐标（Y轴）的值共同定位，直观反映变量间的关联性、分布模式或趋势。
电商运营：价格与销量关系
确定最优定价区间（如$20时销量180件，总收入$3600 vs $30时销量150件，总收入$4500）
发现异常点：若$40时突然销量200件，可能需排查数据错误或促销活动影响
SpotPlot([10,20,30,40,50],[200,180,150,120,80])
医疗研究：药物剂量与疗效
确定安全剂量窗口（100–150mg）
警示毒性风险：>150mg时疗效骤降
SpotPlot([0,50,100,150,200],[10,40,75,80,30])
气候科学：CO₂浓度与温度
量化碳排放对变暖的影响（如CO₂每升50ppm，温度+0.8℃）
支持气候政策制定（如控制CO₂<450ppm以维持温度增幅<2℃）
SpotPlot([300,350,400,450],[14.0,14.5,15.2,15.8])
教育评估：学习时间与成绩
建议最佳学习时长（15-20小时）
识别低效学生：如投入25小时成绩仅82分，需学习方法干预
SpotPlot([5,10,15,20,25],[60,75,85,88,82])
工业制造：设备温度与故障率
设置温度报警阈值（80℃）
预测设备寿命：100℃时故障率15次/月，需提前更换
SpotPlot([60,70,80,90,100],[1,2,3,8,15])
金融分析：广告投入与ROI
优化预算：>30万时边际收益递减
止损点：40万后ROI<150%，应削减投入
SpotPlot([10,20,30,40,50],[120,150,170,160,140])
结构秩
r = sprank(A) 计算稀疏矩阵 A 的结构秩.
A — 输入矩阵,稀疏矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import maximum_bipartite_matching
def sparse_rank_matrix(input_str):
"""
计算稀疏矩阵的结构秩，对标 MATLAB 的 sprank 函数
参数:
input_str: 字符串形式的矩阵输入，例如 "[[1, 0], [0, 1]]"
返回:
结构秩数值（整数）或错误信息字符串
实现原理:
1. 通过二分图匹配计算最大匹配数
2. 使用 SciPy 的 maximum_bipartite_matching 算法
3. 结构秩 = 最大匹配数
注意:
- 不支持符号矩阵
- 自动忽略复数元素
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查转换结果是否有效（sympify 转换失败返回元组）
if isinstance(expr, tuple):
error = True
else:
# 转换为 SymPy 矩阵
a = sp.Matrix(expr) if isinstance(expr, list) else None
if a is not None:
# 检查是否包含符号或复数
contains_symbols = any(element.free_symbols for element in a)
contains_complex = any(element.has(sp.I) for element in a)
if contains_symbols or contains_complex:
error = True
else:
# 转换为 SciPy 稀疏矩阵
A = csr_matrix(np.array(a, dtype=float))
try:
# 正确调用方式（SciPy≥1.4.0）
matching = maximum_bipartite_matching(A)
# 修正：直接统计匹配结果中不等于-1的元素个数
return int(np.sum(matching != -1))
except Exception as e:
result = sp.Matrix(A.toarray()).rank()
return result
else:
error = True # 非矩阵输入
# 生成错误信息
if error:
if contains_symbols:
return f"错误: 不支持符号矩阵输入"
elif contains_complex:
return f"错误: 不支持复数矩阵输入"
else:
return f"输入错误: {input_str}"
return result
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 满秩对角矩阵
m1 = "[[1,0,2,0],[2,0,4,0]]"
print(f"示例1 输入: {m1}")
print(f"示例1 输出: {sparse_rank_matrix(m1)}\n")
# 2
# 示例2: 秩亏缺但结构满秩
m2 = "[[1, 1], [1, 1]]"
print(f"示例2 输入: {m2}")
print(f"示例2 输出: {sparse_rank_matrix(m2)}\n")
# 2
# 示例3: 结构秩小于数值秩
m3 = "[[0, 1], [1, 0]]"
print(f"示例3 输入: {m3}")
print(f"示例3 输出: {sparse_rank_matrix(m3)}\n")
# 2
# 示例4: 非方阵
m4 = "[[1, 0, 0], [0, 1, 0]]"
print(f"示例4 输入: {m4}")
print(f"示例4 输出: {sparse_rank_matrix(m4)}\n")
# 2
创建稀疏矩阵
r = sprank(A) 计算稀疏矩阵 A 的结构秩.
S = sparse(A) 通过挤出任何零元素将满矩阵转换为稀疏格式. 如果矩阵包含许多零, 将矩阵转换为稀疏存储空间可以节省内存.
S = sparse(m,n) 生成 m×n 全零稀疏矩阵.
S = sparse(i,j,v) 根据 i、j 和 v 三元组生成稀疏矩阵 S, 以便 S(i(k),j(k)) = v(k). max(i)×max(j) 输出矩阵为 length(v) 个非零值元素分配了空间. 如果输入 i, j 和 v 为向量或矩阵, 则它们必须具有相同数量的元素. 参量 v 和/或 i 或 j 其中一个参量可以使标量
A — 输入矩阵, 满矩阵, 稀疏矩阵
i,j — 下标对组（指定为单独的参量）, 标量, 向量, 矩阵
v — 值, 标量, 向量, 矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.sparse import csr_matrix
def sparse_matrix(input_str):
"""
对标 MATLAB 的 sparse 函数，创建稀疏矩阵
支持的输入格式:
1. 三数组格式: (行索引列表, 列索引列表, 值列表) 例如: "([0,1], [0,1], [1,2])"
2. 全矩阵转换: 嵌套列表或 SymPy 矩阵 例如: "[[1,0], [0,2]]"
3. 指定维度空矩阵: (行数, 列数) 例如: "(2,3)"
参数:
input_str: 输入字符串，描述矩阵的格式
返回:
SymPy 矩阵对象或错误信息字符串
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 情况1：处理三数组格式 (i, j, v)
if isinstance(expr, tuple):
if len(expr) >= 3 and all(isinstance(item, list) for item in expr[:3]):
# 提取前三个元素作为行索引、列索引和值
i = np.array(expr[0], dtype=int)
j = np.array(expr[1], dtype=int)
v = np.array(expr[2], dtype=float)
# 创建稀疏矩阵（自动处理重复索引的求和）
sparse_mat = csr_matrix((v, (i, j)))
result = sparse_mat.toarray()
# 情况2：处理指定维度的空矩阵 (m, n)
elif len(expr) == 2 and all(e.is_integer for e in expr):
rows, cols = map(int, expr)
sparse_mat = csr_matrix((rows, cols))
result = sparse_mat.toarray()
else:
error = True
# 情况3：处理全矩阵转换
elif isinstance(expr, list):
# 将 SymPy 矩阵转换为 NumPy 数组
dense_array = np.array(expr, dtype=float)
# 转换为稀疏格式后再转回密集数组（用于验证）
sparse_mat = csr_matrix(dense_array)
result = sparse_mat.toarray()
else:
error = True
# 返回 SymPy 矩阵或错误信息
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：三数组格式创建对角矩阵
print("示例1:")
print(sparse_matrix("([0, 1], [0, 1], [1, 2])"))
# Matrix([[1.00000000000000, 0],
# [0, 2.00000000000000]])
# 示例2：全矩阵转换
print("\n示例2:")
print(sparse_matrix("[[1, 0], [0, 2]]"))
# Matrix([[1.00000000000000, 0],
# [0, 2.00000000000000]])
# 示例3：创建 2x3 空矩阵
print("\n示例3:")
print(sparse_matrix("(2, 3)"))
# Matrix([[0, 0, 0],
# [0, 0, 0]])
# 示例4：重复索引自动求和
print("\n示例4:")
print(sparse_matrix("([0, 0], [0, 0], [1, 2])"))
# Matrix([[3.00000000000000]])
三次方样条数据插值
s = spline(x,y,xq) 返回与 xq 中的查询点对应的插值 s 的向量。s 的值由 x 和 y 的三次样条插值确定。
pp = spline(x,y) 返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp。
x - x 坐标, 向量
y - x 坐标处的函数值, 向量 | 矩阵
xq - 查询点, 标量 | 向量 | 矩阵
s - 查询点位置的插值, 标量 | 向量 | 矩阵
pp - 分段多项式,结构体
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.interpolate import CubicSpline
def spline_cubic(input_str):
"""
对标MATLAB的spline函数，实现三次样条插值
参数:
input_str: 输入参数字符串，格式示例：
"x, y" -> 返回样条函数
"y" -> 自动生成x坐标并返回样条函数
"x, y, xq" -> 返回xq处的插值结果
"y, xq" -> 自动生成x坐标并返回xq处插值结果
返回:
样条函数对象或插值结果数组，错误时返回字符串信息
示例:
>>> spline_cubic("[1,2,3], [4,5,6]")
>>> spline_cubic("[1,2,3], [4,5,6], 2.5")
array(5.0)
"""
try:
# 符号解析与类型转换 ------------------------------------------------
def convert(param):
"""将SymPy对象转换为数值数组"""
if isinstance(param, list):
matrix = sp.Matrix(param)
if matrix.is_symbolic():
raise ValueError(f"Symbolic variable {param} not allowed")
else:
return np.array(matrix.tolist(), dtype=float).flatten()
elif param.is_Number:
return np.array([float(param)])
elif param.free_symbols:
raise ValueError(f"Symbolic variable {param} not allowed")
else: # 处理包含符号的表达式
try:
return np.array([float(param.evalf())])
except:
raise ValueError(f"Cannot convert {param} to numeric")
# 解析输入字符串
expr = sp.sympify(input_str, evaluate=False) # 强制转换为元组
params = [convert(p) for p in expr] if isinstance(expr, tuple) else [convert(expr)]
# 参数处理与插值计算 ------------------------------------------------
def validate_length(a, b):
if len(a) != len(b):
raise ValueError(f"Length mismatch: {len(a)} vs {len(b)}")
# 处理不同参数数量
if len(params) == 1: # spline(y)
y = params[0]
x = np.arange(len(y))
return CubicSpline(x, y, bc_type='not-a-knot')
elif len(params) == 2: # spline(x,y) 或 spline(y,xq)
a, b = params
if len(a) == len(b): # x,y模式
return CubicSpline(a, b, bc_type='not-a-knot')
else: # y,xq模式
y, xq = a, b
x = np.arange(len(y))
return CubicSpline(x, y)(xq)
elif len(params) == 3: # spline(x,y,xq)
x, y, xq = params
validate_length(x, y)
return CubicSpline(x, y, bc_type='not-a-knot')(xq)
else:
raise ValueError("Invalid number of parameters (1-3 required)")
except Exception as e:
return f"Error: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 标准插值测试
print(spline_cubic("[1,2,3], [4,5,6], 2"))
# [5.]
print(spline_cubic("[0,1,2], [0,1,4], 1.5"))
# [2.25]
# 自动生成x坐标测试
print(spline_cubic("[3,1,4], [0.5, 1.5]"))
# [1.375 1.875]
平方根
B=sqrt（X）返回数组X中每个元素的平方根.对于X中为负或复数的元素，sqrt（X）会产生复数结果.
X —— 输入数组,标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def square_root(input_str):
"""
计算元素的平方根，对标 MATLAB 的 sqrt 函数。
参数:
input_str: 字符串形式的输入，可以是标量或矩阵，例如 "4" 或 "[[1, 4], [9, 16]]"
返回:
计算结果（SymPy 矩阵或数值），错误时返回字符串描述错误信息
示例:
>>> square_root("4")
2.0
>>> square_root("[[1, 4], [9, 16]]")
Matrix([
[1.0, 2.0],
[3.0, 4.0]])
>>> square_root("[[a, 4], [9, b]]")
Matrix([
[sqrt(a), 2.0],
[ 3.0, sqrt(b)]])
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查转换结果是否有效（sympify 转换失败返回元组）
if isinstance(expr, tuple):
error = True
else:
# 标量处理：直接计算平方根并数值化
result = sp.sqrt(expr).evalf()
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例 1: 标量输入
print("示例1 输入: '4'")
print("示例1 输出:", square_root("4"), end="\n\n")
# 2.00000000000000
# 示例 2: 复数结果
print("示例2 输入: '-1'")
print("示例2 输出:", square_root("-1"), end="\n\n")
# 1.0*I
矩阵平方根
X = sqrtm(A)返回矩阵A的主平方根,即X*X=A.
X是每个特征值都有非负实部的唯一平方根.如果A有任何具有负实部的特征值,则会产生复结果.
如果A是单数,那么A可能没有平方根.
A —— 输入矩阵,方阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy import linalg as dense_linalg
def square_root_matrix(input_str):
"""
计算矩阵的平方根，对标 MATLAB 的 sqrtm 函数。
参数:
input_str: 字符串形式的矩阵输入，例如 "[[1, 0], [0, 4]]"。
返回:
如果输入合法且为方阵，返回 SymPy 矩阵；否则返回错误信息字符串。
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查转换后的表达式是否有效（如果转换失败，sympify 可能返回元组）
if isinstance(expr, tuple):
error = True
else:
# 转换为 SymPy 矩阵
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is not None:
# 检查矩阵是否包含符号或虚数单位 I
contains_symbols = A.is_symbolic() or any(element.has(sp.I) for element in A)
# 仅处理方阵
if A.is_square:
if contains_symbols:
# 符号矩阵：使用 SymPy 的矩阵幂运算
result = A ** sp.Rational(1, 2)
else:
# 数值矩阵：转换为 NumPy 数组并用 SciPy 计算平方根
N_A = np.array(A, dtype=float)
scipy_result = dense_linalg.sqrtm(N_A)
# 将结果转换回 SymPy 矩阵并四舍五入到合理精度
result = sp.Matrix(scipy_result.round(10))
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: 数值方阵
matrix1 = "[[1, 0], [0, 4]]"
print("示例1 输入:", matrix1)
print("示例1 输出:", square_root_matrix(matrix1))
# Matrix([[1.00000000000000, 0],
# [0, 2.00000000000000]])
# 示例 2: 符号矩阵
a, b = sp.symbols('a b')
matrix2 = f"[[{a}, 0], [0, {b}]]"
print("\n示例3 输入:", matrix2)
print("示例3 输出:", square_root_matrix(matrix2))
# Matrix([[sqrt(a), 0],
# [0, sqrt(b)]])
# 示例 3: 复数矩阵
matrix3 = "[[-1, 0], [0, -1]]"
print("\n示例4 输入:", matrix3)
print("示例4 输出:", square_root_matrix(matrix3))
# Matrix([[1.0*I, 0],
# [0, 1.0*I]])
格式化距离矩阵
ZOut = squareform(yIn) 将yIn（m个观测值的长度为m（m–1）/2的成对距离向量）转换为ZOut（对角线上有零的m乘m对称矩阵）
yIn是数值向量,逻辑向量.
ZOut是数值矩阵,逻辑矩阵
# 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_square_form(input_str):
"""
对标MATLAB的squareform函数，实现压缩向量与距离方阵的相互转换。
参数:
input_str (str): 输入字符串，可以是以下两种形式：
1. 压缩距离向量：如 "[0.5, 1.2, 2.3]"
2. 对称距离方阵：如 "[[0, 0.5, 1.2], [0.5, 0, 2.3], [1.2, 2.3, 0]]"
返回:
sympy.Matrix: 根据输入类型自动转换：
- 输入为向量时返回对称方阵
- 输入为方阵时返回压缩向量
如果输入格式错误，返回描述性错误信息
示例:
>>> distance_square_form("[1, 2, 3]")
Matrix([
[0, 1, 2],
[1, 0, 3],
[2, 3, 0]])
>>> distance_square_form("[[0, 1, 2], [1, 0, 3], [2, 3, 0]]")
Matrix([1.0, 2.0, 3.0])
"""
try:
# 将输入字符串解析为SymPy对象
expr = sp.sympify(input_str)
error = False
result = None
# 处理列表输入
if isinstance(expr, list):
arr = np.array(expr, dtype=float).ravel()
else:
error = True
# 检查数组维度有效性
if arr.ndim not in (1, 2):
raise ValueError("输入必须是向量(1D)或矩阵(2D)")
# 使用scipy进行格式转换
result = distance.squareform(arr)
# 转换为SymPy矩阵返回
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except ValueError as ve:
return f"数值错误: {str(ve)}"
except sp.SympifyError:
return "格式错误: 无法解析输入字符串"
except Exception as e:
return f"处理错误: {str(e)}"
input_str = "[1, 2, 3, 4, 5, 6]"
print(distance_square_form(input_str))
# Matrix([[0, 1.00000000000000, 2.00000000000000, 3.00000000000000],
# [1.00000000000000, 0, 4.00000000000000, 5.00000000000000],
# [2.00000000000000, 4.00000000000000, 0, 6.00000000000000],
# [3.00000000000000, 5.00000000000000, 6.00000000000000, 0]])
方波
x = square(t) 为时间数组t的元素生成周期为2π的方波.square类似于正弦函数,但会创建值为-1和1的方波.
x = square(t,duty) 生成指定占空比为duty的方波.占空比是方波为正的信号周期的百分比.
t — 时间数组,向量,矩阵
duty — 占空比,0到1之间
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy import signal
def square_wave(input_str):
"""
生成方波信号，对标MATLAB square函数
参数:
input_expr: 输入表达式，支持以下形式：
- 时间向量 (数值数组/矩阵)
- 符号表达式 (如 't'，需包含单个符号变量)
- 元组 (时间向量, 占空比)
duty (可选): 占空比 (0到1之间，默认0.5)
freq (可选): 信号频率 (单位Hz，默认1)
sample_points (可选): 符号表达式采样点数 (默认500)
返回:
tuple: (时间向量, 方波信号) 或 错误信息字符串
示例:
>>> # 数值输入
>>> t = [0, 0.1, 0.2, 0.3]
>>> square_wave(t)
(array([0. , 0.1, 0.2, 0.3]), array([ 1., 1., 1., -1.]))
>>> # 符号表达式输入
>>> square_wave('t', duty=0.3)
(array([0. , 0.002,...,1.0]), array([1, 1,...,-1]))
>>> # 矩阵输入
>>> square_wave([[0, 0.5], [1, 1.5]])
(array([0. , 0.5, 1. , 1.5]), array([ 1., 1., -1., -1.]))
"""
try:
# 统一输入处理
expr = sp.sympify(input_str)
error = False
duty = 0.5
freq = 1
sample_points = 500
# 处理元组输入 (时间向量, 占空比)
if isinstance(expr, tuple) and len(expr) == 2:
t, custom_duty = expr
duty = float(custom_duty)
# 参数校验
if not 0 <= duty <= 1:
raise ValueError("占空比必须在0到1之间")
expr = t
# 矩阵/数组处理
if isinstance(expr, list):
t_values = np.array(expr).astype(float).ravel()
# 符号表达式处理
elif expr.free_symbols:
symbols = expr.free_symbols
if len(symbols) != 1:
raise ValueError("符号表达式必须包含单个变量")
var = symbols.pop()
t_values = np.linspace(0, 1 / freq, sample_points)
expr_func = sp.lambdify(var, expr, 'numpy')
t_values = expr_func(t_values)
# 数值处理
elif expr.is_number:
t_values = np.array(expr, dtype=float)
else:
error = True
# 生成方波 (MATLAB的2π周期对应频率f=1/(2π))
y_values = signal.square(2 * np.pi * freq * t_values, duty=duty)
return t_values, y_values if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 处理各种输入类型
print(square_wave("[[0,0.5],[1,1.5]]"))
# (array([0. , 0.5, 1. , 1.5]),
# array([ 1., -1., 1., -1.]))
删除长度为 1 的维度
B = squeeze(A) 返回一个数组，其元素与输入数组 A 相同，但删除了长度为 1 的维度。例如，如果 A 是 3×1×2 数组，则 squeeze(A) 返回 3×2 矩阵。
如果 A 是行向量、列向量、标量或没有长度为 1 的维度的数组，则 squeeze 返回输入 A。
A — 输入数组,多维数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def squeeze_arr_dim(input_str):
"""
该函数的主要功能是将输入的字符串转换为数组，然后去除数组中长度为 1 的维度。
如果输入不是有效的矩阵格式，或者在处理过程中出现异常，会返回相应的错误信息。
参数:
input_str (str): 输入的字符串，预期为表示矩阵的字符串
返回:
numpy.ndarray 或 str: 如果处理成功，返回去除长度为 1 的维度后的数组；如果出现错误，返回错误信息字符串
"""
try:
# 这个函数可能用于处理字符串中的一些格式问题，使其符合后续处理的要求
# 使用 sympy 的 sympify 函数将修正后的字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 初始化结果变量，用于存储最终处理后的数组
result = None
# 初始化错误标志，用于标记是否出现错误
error = False
# 如果是矩阵，该函数会返回矩阵对象；否则返回 None
# 使用海象运算符（Python 3.8 及以上版本支持）将结果赋值给 matrix 变量
if isinstance(expr,list):
# 如果是矩阵，将矩阵对象转换为 NumPy 数组，数据类型为 object
arr = np.array(expr, dtype=object)
# 使用 NumPy 的 squeeze 函数去除数组中长度为 1 的维度
result = np.squeeze(arr)
else:
# 如果输入不是矩阵，将错误标志设置为 True
error = True
# 如果没有出现错误，返回处理后的数组；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 如果在处理过程中出现异常，捕获异常并返回错误信息
return f"错误:{e}"
print(squeeze_arr_dim("[[[1,2]]]"))
# [1 2]
print(squeeze_arr_dim("[[[x,y]]]"))
# [x y]
状态空间模型
A,B,C,D,dt = ss(sys) 将线性系统转换为空间系统形式。始终创建一个新系统，即使sys已经是一个状态空间系统.
A - 装态数, 状态数是矩阵A的行(或列)数. 如果矩阵A是一个4×4的矩阵，因此有四个状态.
B - 输入数：输入数是矩阵B的列数. 比如矩阵B是一个4×2的矩阵,因此有两个输入.
C - 输出数：输出数是矩阵C的行数. 比如矩阵C是一个3×4的矩阵,因此有三个输出.
D - 矩阵D代表直接传递矩阵（Direct transmission matrix）,也称为前馈矩阵（feedforward matrix）.它描述了系统输入对系统输出的直接影响,而不经过系统的状态.
dt - 采样时间
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import control.matlab as ct
def convert_to_numpy_or_float(expr):
"""
将表达式转换为 NumPy 数组或浮点数。
:param expr: 输入的表达式
:return: 转换后的 NumPy 数组或浮点数
"""
if isinstance(expr, list):
return np.array(expr, dtype=float)
elif expr.is_number:
return float(expr)
else:
# 如果既不是列表也不是数字，抛出类型错误
raise TypeError(f"不支持的类型: {type(expr).__name__}")
def state_space_system(input_str):
"""构建状态空间系统
参数格式示例：
- 连续系统： "[[-1,-2],[3,-4]],[[5],[7]],[[6, 8]],[[9]]"
- 离散系统： "[[-1.5,-2],[1,0]],[[0.5],[0]],[0,1],0,0.2"
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查表达式是否为元组
if not isinstance(expr, tuple):
return f"输入错误: 期望输入为元组，但得到 {type(expr).__name__}"
# 检查元组长度是否至少为 4
if len(expr) < 4:
return f"输入错误: 元组长度至少为 4，但得到 {len(expr)}"
# 提取 A, B, C, D 矩阵
A = convert_to_numpy_or_float(expr[0])
B = convert_to_numpy_or_float(expr[1])
C = convert_to_numpy_or_float(expr[2])
D = convert_to_numpy_or_float(expr[3])
# 检查是否有采样时间
if len(expr) > 4:
try:
dt = float(expr[4])
except ValueError:
return f"输入错误: 采样时间必须是数字，但得到 {expr[4]}"
# 创建连续时间状态空间系统
sys1 = ct.ss(A, B, C, D, dt)
result = sp.Matrix(sys1.A), sp.Matrix(sys1.B), sp.Matrix(sys1.C), sp.Matrix(sys1.D), sys1.dt
else:
# 创建离散时间状态空间系统
sys1 = ct.ss(A, B, C, D)
result = sp.Matrix(sys1.A), sp.Matrix(sys1.B), sp.Matrix(sys1.C), sp.Matrix(sys1.D)
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 测试案例
if __name__ == "__main__":
test_cases = [
"([[-1,-2],[3,-4]],[[5],[7]],[[6, 8]],[[9]])",
# A: Matrix([[-1.00000000000000, -2.00000000000000], [3.00000000000000, -4.00000000000000]])
# B: Matrix([[5.00000000000000], [7.00000000000000]])
# C: Matrix([[6.00000000000000, 8.00000000000000]])
# D: Matrix([[9.00000000000000]])
"([[-1.5,-2],[1,0]],[[0.5],[0]],[0,1],0,0.2)",
# A: Matrix([[-1.50000000000000, -2.00000000000000], [1.00000000000000, 0]])
# B: Matrix([[0.500000000000000], [0]])
# C: Matrix([[0, 1.00000000000000]])
# D: Matrix([[0]])
]
for case in test_cases:
print(f"测试输入: {case}")
result = state_space_system(case)
if isinstance(result, tuple):
if len(result) == 4:
A, B, C, D = result
else:
A, B, C, D, _ = result
print("A:", A)
print("B:", B)
print("C:", C)
print("D:", D)
else:
print("错误信息:", result)
print("─" * 50)
移位正弦积分
ssinint(X)返回移位后的正弦积分函数. ssiint(X)=sinint(X)-pi/2.
X —— 输入,符号数,符号变量,符号表达式,符号函数,符号向量,符号矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Si, pi
import numpy as np
from scipy.special import sici
def shifted_sine_integral(input_str):
"""
对标 MATLAB 的 ssinint 函数，计算位移正弦积分 Si(x) - π/2。
参数:
input_str: 输入的字符串，可以表示数值、符号表达式或矩阵。
返回:
- 对于数值输入，返回浮点数结果。
- 对于符号输入，返回符号表达式。
- 对于矩阵输入，返回各元素对应的位移正弦积分结果的 SymPy 矩阵。
- 错误时返回错误消息字符串。
示例:
>>> shifted_sine_integral('2')
-0.110538260536646 # 示例值，实际结果可能不同
>>> shifted_sine_integral('x')
Si(x) - pi/2
>>> shifted_sine_integral('[[1, x], [2, y]]')
Matrix([
[Si(1) - pi/2, Si(x) - pi/2],
[Si(2) - pi/2, Si(y) - pi/2]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理单个数值或符号
if expr.is_number:
z = complex(expr)
"""计算位移正弦积分 Si(x) - π/2，仅处理数值标量输入。"""
si, _ = sici(z) # 仅取正弦积分部分
result = si - np.pi / 2
elif expr.free_symbols:
"""计算位移正弦积分 Si(x) - π/2，根据输入类型返回数值或符号表达式。"""
result = Si(expr) - pi / 2
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 数值输入
print(shifted_sine_integral('2'))
# (0.03461665000779801-0j)
# 符号输入
print(shifted_sine_integral('x'))
# Si(x) - pi/2
stairs(Y) 绘制 Y 中元素的阶梯图.
如果Y为矩阵，则stairs为每个矩阵列绘制一个线条
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.interpolate import interp1d
def stair_wave(input_str):
"""
对标 MATLAB 的 stairs 函数，生成阶梯图数据。
参数:
input_str: 输入的字符串，可以表示：
- 数值矩阵（视为 Y 值，X 自动生成）
- 符号表达式（如 'sin(x)'，需配合 x_range 使用）
x_range: 符号表达式时 X 的取值范围，格式为 (start, stop, num_points)
返回:
- 成功时返回 (x_out, y_out) 的 NumPy 数组
- 错误时返回错误消息字符串
示例:
>>> x, y = stair_wave('[[1, 2, 3], [4, 5, 6]]') # 输入二维矩阵
>>> x, y = stair_wave('sin(x)', x_range=(0, 2*np.pi, 100)) # 符号表达式
"""
try:
# 尝试将输入解析为 SymPy 表达式或矩阵
expr = sp.sympify(input_str)
x_out, y_out = None, None
error = False
# Continuous data
NUM = 1000 # modify the number of points you want
if isinstance(expr, list):
# 如果 expr 是一个矩阵，则将其转换为 NumPy 数组并展平
x = np.array(expr, dtype=float).ravel()
# 这里 y 的计算暂时写为 np.sin(x)，实际使用时 y 需要根据具体需求输入，例如 0.5*cos(x)
y = np.sin(x)
# 对 y 进行累加操作
y = np.cumsum(y)
# 初始化 x=0 时的 y 值
y0 = np.array([0])
# 通过插值计算 y=0 时对应的 x 值
x0 = np.interp([0], y, x)
# 将 x0 与原 x 数组拼接
x = np.concatenate([x0, x])
# 将 y0 与原 y 数组拼接
y = np.concatenate([y0, y])
# 使用 'next' 方式进行插值，创建插值函数
funct = interp1d(x, y, kind='next')
# 生成连续的 x 值
x_cont = np.linspace(x[0], x[-1], NUM)
# 通过插值函数计算对应的连续 y 值
y_cont = funct(x_cont)
elif expr.free_symbols:
# 如果输入是符号表达式，默认生成 0 到 4π 之间的 50 个 x 值
x_values = np.linspace(0, 4 * np.pi, 50)
# 如果 expr 是一个 SymPy 表达式，则执行相应的操作
if variables := expr.free_symbols:
# 将 SymPy 表达式转换为可以在 NumPy 中计算的函数
expression_func = sp.lambdify(tuple(variables), expr, modules=['numpy'])
# 计算 y
y = expression_func(x_values)
else:
# 如果没有自由变量，直接将 expr 作为 y
y = expr
# 对 y 进行累加操作
y = np.cumsum(y)
# 初始化 x=0 时的 y 值
y0 = np.array([0])
# 通过插值计算 y=0 时对应的 x 值
x0 = np.interp([0], y, x_values)
# 将 x0 与原 x 数组拼接
x_values = np.concatenate([x0, x_values])
# 将 y0 与原 y 数组拼接
y = np.concatenate([y0, y])
# 使用 'next' 方式进行插值，创建插值函数
funct = interp1d(x_values, y, kind='next')
# 生成连续的 x 值
x_cont = np.linspace(x_values[0], x_values[-1], NUM)
# 通过插值函数计算对应的连续 y 值
y_cont = funct(x_cont)
else:
# 如果 expr 不是矩阵也不是 SymPy 表达式，则抛出错误或者执行其他操作
error = True
return x_cont, y_cont if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 示例 : 符号表达式输入
x, y = stair_wave('sin(x)')
print("符号表达式采样:", x[:5], y[:5])
# [0. 0.01257895 0.0251579 0.03773685 0.0503158 ]
# [0. 0.25365458 0.25365458 0.25365458 0.25365458]
插入标准缺失值
B = standardizeMissing(A,indicator) 用 A 中的标准缺失值替换 indicator 中指定的值，并返回一个标准化数组
A — 输入数据,向量,矩阵
indicator — 非标准缺失值指示符,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from ast import literal_eval # 用于安全解析元组
def standardize_missing(input_str):
"""
对标 MATLAB 的 standardizeMissing 函数，将指定值替换为标准缺失值 NaN
参数:
input_str: 输入的字符串，格式应为 "(matrix, indicator)"，其中：
- matrix: 数值矩阵（支持嵌套列表或 SymPy 矩阵）
- indicator: 要替换的值（数值或列表）
返回:
SymPy 矩阵 或 错误信息字符串
示例:
>>> standardize_missing("([[1, 0], [3, -999]], -999)")
Matrix([[1, 0], [3, nan]])
>>> standardize_missing("([5, 2, 888, 4], [2, 888])")
Matrix([[5, nan, nan, 4]])
"""
try:
# 安全解析输入字符串为Python对象
parsed = literal_eval(input_str.strip())
# 验证输入格式是否为 (matrix, indicator)
if not (isinstance(parsed, tuple) and len(parsed) == 2):
return "输入错误: 必须为 (matrix, indicator) 格式的元组"
matrix, indicator = parsed
# 转换输入矩阵为SymPy矩阵
sp_matrix = sp.Matrix(matrix) if isinstance(matrix, list) else None
if sp_matrix is None:
return "输入错误: 无法解析矩阵数据"
# 验证矩阵是否为纯数值型
if sp_matrix.is_symbolic():
return "错误: 矩阵包含符号元素，仅支持数值矩阵"
# 转换指示符为NumPy数组形式
indicators = np.array([indicator]).flatten() if np.isscalar(indicator) else np.array(indicator)
# 转换为NumPy数组进行操作
np_matrix = np.array(sp_matrix.tolist(), dtype=float)
# 创建掩码矩阵（True表示需要替换的位置）
mask = np.isin(np_matrix, indicators)
# 替换目标值为NaN
np_matrix[mask] = np.nan
# 转换回SymPy矩阵
if sp_matrix.shape[1] == 1:
return sp.Matrix(np_matrix).T
else:
return sp.Matrix(np_matrix)
except SyntaxError:
return "语法错误: 输入格式应为有效的Python元组表达式"
except ValueError as ve:
return f"数值转换错误: {ve}"
except Exception as e:
return f"未处理的异常: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试用例1：替换单个值
test1 = "([[1, 0], [3, -999]], -999)"
print(f"测试1输入: {test1}")
print("测试1输出:", standardize_missing(test1))
# Matrix([[1.00000000000000, 0],
# [3.00000000000000, nan]])
# 测试用例2：替换多个值
test2 = "([5, 2, 888, 4], [2, 888])"
print(f"\n测试2输入: {test2}")
print("测试2输出:", standardize_missing(test2))
# Matrix([[5.00000000000000, nan, nan, 4.00000000000000]])
标准差
S = std(A) 返回数组元素的标准偏差,即分布范围的度量.默认情况下,标准偏差是为展开的阵列计算的,否则是在指定的轴上计算的.
S = std(A,dim). 计算标准偏差的一个或多个轴.默认情况是计算展平阵列的标准偏差.
dim=0,计算以列为单位的标准差.
dim=1,计算以行为单位的标准差.
A - 输入数组,向量,矩阵
dim - 沿其运算的维度
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def standard_deviation_array(input_str, ddof=1):
"""
对标 MATLAB 的 std 函数，精确处理向量和矩阵输入
参数:
input_str: 输入字符串，格式支持：
- "matrix" : 直接输入矩阵
- "(matrix, dim)" : 矩阵+维度（1=列方向，2=行方向）
ddof: 自由度调整（默认1对应无偏估计）
返回:
- 标量（向量输入或整体计算）
- 列表（矩阵按列/行计算）
- 错误信息字符串
示例:
>>> standard_deviation_array("[1, 2, 3, 4]") # 向量输入
1.118033988749895
>>> standard_deviation_array("[[1,2],[3,4]]") # 矩阵默认按列
[1.0, 1.0]
>>> standard_deviation_array("([[1,2],[3,4]], 2)") # 按行计算
[0.7071, 0.7071]
"""
try:
# 解析输入
parsed = sp.sympify(input_str)
matrix_data, dim = (parsed if isinstance(parsed, tuple)
else (parsed, None))
# 验证矩阵有效性
sp_matrix = sp.Matrix(matrix_data) if isinstance(matrix_data, list) else None
if not sp_matrix or sp_matrix.is_symbolic():
return "无效矩阵输入" if not sp_matrix else "仅支持数值矩阵"
# 转换为 NumPy 数组
np_matrix = np.array(sp_matrix.tolist(), dtype=float)
original_shape = np_matrix.shape
# 自动识别向量并展平
if np_matrix.ndim == 2 and min(np_matrix.shape) == 1:
np_matrix = np_matrix.ravel()
# 计算维度处理
if dim is None:
# 智能维度判断
if np_matrix.ndim == 1:
axis = None # 向量计算整体
else:
axis = 0 # 矩阵默认按列
else:
axis = int(dim) - 1 # MATLAB 转 NumPy 轴
# 执行标准差计算
result = np.std(np_matrix, axis=axis, ddof=ddof)
# 标准化输出格式
if isinstance(result, np.ndarray):
return np.around(result.tolist(), 4) if result.ndim == 1 else result.tolist()
else:
return float(np.around(result, 4))
except Exception as e:
return f"错误: {str(e)}"
# 验证测试
if __name__ == "__main__":
# 测试1: 行向量输入
print(standard_deviation_array("[[1, 2, 3, 4]]"))
# 1.291
# 测试2: 列向量输入
print(standard_deviation_array("[[1], [2], [3], [4]]"))
# 1.291
# 测试3: 4x3矩阵按列
test_data = "[[1,5,9], [2,6,10], [3,7,11], [4,8,12]]"
print(standard_deviation_array(test_data))
# [1.291 1.291 1.291]
# 测试4: 同上矩阵按行
print(standard_deviation_array(f"({test_data}, 2)"))
# [4. 4. 4. 4.]
阶跃响应图
stepplot(sys) 可视化系统对单位阶跃输入（Unit Step Input）的动态响应.
单位阶跃输入在时间 t = 0时从 0 突然跳变到 1（例如，一个开关的瞬时开启）。
阶跃响应图显示了系统输出随时间的变化，帮助我们分析系统的稳定性、速度、准确性等性能指标。
一阶系统（无振荡）:
类似 RC 电路充电过程（如电容器电压响应）
阶跃响应特性：
无超调（单调上升）
上升时间：约 1.15 秒（从 10% 到 90% 稳态值）
调节时间：约 4 秒（±2% 准则）
稳态值：2.5（因为 G(0) = 5/2 = 2.5G(0)=5/2=2.5)
stepplot(5,s+2)
二阶欠阻尼系统（有振荡）:
弹簧-质量-阻尼系统的振荡响应（如汽车悬架）
阶跃响应特性：
超调量：约 16.3%（因阻尼比 \zeta = 0.316ζ=0.316)
上升时间：约 0.55 秒
峰值时间：约 1.05 秒
调节时间：约 4 秒（±2% 准则）
稳态值：1（因为 G(0) = 10/10 = 1G(0)=10/10=1)
stepplot(10,s^2+2s+10)
二阶临界阻尼系统（最快无振荡）:
门关闭装置的理想响应（快速且无反弹）
阶跃响应特性：
无超调（单调上升）
上升时间：约 1.17 秒
调节时间：约 2.4 秒（±2% 准则）
稳态值：1（因为 G(0) = 9/9 = 1G(0)=9/9=1)
stepplot(9,s^2+6s+9)
高阶系统（含振荡和延迟效应）:
带滤波器的电机位置控制系统
阶跃响应特性：
超调量：约 8%
调节时间：约 8 秒
稳态值：0.4（因为 G(0) = 2/5 = 0.4G(0)=2/5=0.4)
stepplot(s+2,s^3+5s^2+9s+5)
质量-弹簧-阻尼系统（机械工程）
汽车减震器设计、建筑抗震分析。
阶跃响应特性：
输入：突加 1 N 的恒定力（阶跃力）
输出：物体位移
超调量 ≈ 15%（轻微振荡）
调节时间 ≈ 4 秒
stepplot(1,0.5s^2+1.2s+8)
水箱液位控制系统（过程控制）
化工生产中的液位调节，需避免溢出。
阶跃响应特性：
输入：进水流量阶跃增加（如 0→2 m³/s）
输出：液位高度（稳态值 = 1 m）
上升时间 ≈ 11 秒（缓慢上升，无超调）
stepplot(2,5s+1)
直流电机位置控制（机电系统）
机器人关节定位，要求精确无振荡。
阶跃响应特性：
输入：目标角度阶跃（如 0→90°）
输出：实际转子角度
稳态值 = 90°（无稳态误差）
上升时间 ≈ 1.2 秒（单调上升）
stepplot(3,s*(s+4))
RLC带通滤波器（电子工程）:
通信系统中提取高频信号。
阶跃响应特性：
输入：电压阶跃（0→1 V）
输出：电阻两端电压
初始峰值 ≈ 0.8 V（快速响应）
稳态值 = 0（高通特性）
stepplot(0.6s,s^2+0.7s+1)
疫情传播模型（流行病学）
预测隔离政策对感染峰值的影响。
阶跃响应特性：
输入：突现 1 名感染者（阶跃输入）
输出：累计感染比例
峰值感染率 ≈ 20%（超调）
稳定时间 ≈ 30 天
stepplot(0.05,s^2+0.3s+0.05)
股票市场反应模型（金融）
量化交易中事件驱动策略的回测。
阶跃响应特性：
输入：消息公布（阶跃事件）
输出：股价变动
初始过冲 ≈ +15%（投机泡沫）
稳态值 = +100%（最终合理估值）
stepplot(-0.2s+1,s^2+0.5s+1)
卫星姿态控制（航空航天）:
太空望远镜精密指向控制。
阶跃响应特性：
输入：目标角度阶跃（如 0→10°）
输出：实际俯仰角
稳态误差 = 0（积分作用）
超调 ≈ 35%（需优化阻尼）
stepplot(100,s^2*(s+10))
机械减震系统
观察车身振动的超调量、稳定时间，优化舒适性。
汽车悬架系统，输入为路面阶跃冲击（如突然遇到障碍物），输出为车身垂直位移。
分母 s² + 2s + 5 对应质量-弹簧-阻尼系统
s² 项：质量惯性（m）
2s 项：阻尼系数（c）
5：弹簧刚度（k）
stepplot([1],[1,2,5])
温度控制系统
分析水温上升的延迟和稳态时间，避免过热风险。
电热水壶的加热系统，输入为阶跃电压（突然通电），输出为水温。
分母 3s + 1 表示一阶惯性：
3：热容（C）与热阻（R）的乘积（时间常数 τ = RC）
stepplot([1],[3,1])
直流电机调速
检查转速是否快速跟踪指令（上升时间短），且无振荡（超调小）。
电机转速控制，输入为阶跃电压，输出为转速。
分母 s² + 4s + 8 包含： 电感（L）、电阻（R）、反电动势（K）的耦合效应
stepplot([1],[1,4,8])
化学反应釜浓度控制
评估浓度超调是否导致副产物超标，调整进料策略。
反应物输入流量阶跃增加，输出为产物浓度。
三阶系统（s³项）反映多级混合/反应延迟。
stepplot([1],[1,3,3,1])
生态系统种群模型
预测种群是否稳定或震荡，指导生态干预。
突然引入食物资源（阶跃输入），输出为物种数量（如鱼类）。
低阻尼（0.5s项小）易引发振荡（数量暴涨后崩溃）。
stepplot([1],[1,0.5,0.1])
叠加阶跃响应图
stepplot([sys1],[sys2]...[sysN]) 是在同一坐标系中绘制多个系统对单位阶跃输入的响应曲线，用于直观比较不同系统的动态性能。
不同阻尼比的弹簧系统:
质量相同的弹簧系统，阻尼系数递增。
stepplot([1,s^2+0.5s+1], # ζ=0.25 (欠阻尼)
[1,s^2+s+1], # ζ=0.5 (最佳阻尼)
[1,s^2+2s+1]) # ζ=1.0 (临界阻尼)
电路滤波器性能对比
电子滤波器对电压阶跃的响应
stepplot([1,s+1], # 一阶低通
[s,s^2+s+1], # 带通
[1,0.1s^2+0.3s+1]) # 二阶快速响应
温度控制器增益对比
恒温箱加热系统（时间常数5秒）
stepplot([0.5,5s+1], # 低增益
[1,5s+1], # 中增益
[2,5s+1]) # 高增益
无人机高度控制器
stepplot([4,s^2+0.8s+4], # 标准
[6,s^2+1.2s+6], # 高刚度
[3,s^2+0.4s+3]) # 低阻尼
机械臂关节定位精度
stepplot([1,s^2+3s], # PD控制
[s+1,s^3+3s^2+2s], # PID控制
[0.5,s^2+2s+1]) # 简化模型
流线图
StreamPlot(expression)绘制一系列平滑的、永不交叉的曲线（即“流线”）来描绘二维向量场的图像。
1：流体力学 - 一个点源（水龙头滴下一滴水）
想象在原点有一个水龙头，水均匀地向四面八方流出去。
画图范围：x从-2到2，y从-2到2（注意避开原点x=0, y=0，因为那里分母为零）。
它像什么：所有流线都是从原点（源点）直接向外发射的直线。
StreamPlot(x/(x^2+y^2),y/(x^2+y^2),x=[-2,2],y=[-2,2])
2：流体力学 - 一个涡旋（水池排水时形成的漩涡）
想象原点是一个排水口，水在旋转着流下去。
画图范围：x从-2到2，y从-2到2（同样避开原点）。
它像什么：所有流线都是围绕原点旋转的圆圈，形成一个涡旋。
StreamPlot(x/(x^2+y^2),-y/(x^2+y^2),x=[-2,2],y=[-2,2])
3：电磁学 - 一个偶极子（一个正电荷和一个负电荷靠在一起）
就像一根磁铁或一个电荷对的磁场/电场。
画图范围：x从-2到2，y从-1到1。
它像什么：流线从正电荷（源）出发，弯曲地连接并终止于负电荷（汇）。
StreamPlot((x+0.5)/(((x+0.5)^2+y^2 )^1.5)-(x-0.5)/(((x-0.5)^2+y^2)^1.5),y/(((x+0.5)^2+y^2)^1.5)-y/(((x-0.5)^2+y^2)^1.5),x=[-2,2],y=[-1,1])
4：理论模型 - saddle点（鞍点）
这是一个非常基础但又重要的场结构，像一个山口。
画图范围：x从-2到2，y从-2到2。
它像什么：流线沿着x轴向外散开，沿着y轴向内收拢。原点(0,0)这个点就叫做saddle点（鞍点），它既不是源也不是汇，是不稳定点。
StreamPlot(x,-y,x=[-2,2],y=[-2,2])
5：最简单的函数 F(z) = z
表达式：F(z) = z = x + i*y
对应的向量场：( u, v ) = ( x, y )
流线图分析：
这是一个标准的源（Source）。
所有点的向量都从原点指向外。
离原点越远，向量长度（模）越大。
这类似于一个水龙头向平静的水面滴水形成的流场。
StreamPlot(z,z=[-2-2@i,2+2@i])
6：F(z) = i * z （乘以虚数单位 i）
表达式：F(z) = i * (x + i*y) = -y + i*x
对应的向量场：( u, v ) = ( -y, x )
流线图分析：
这是一个纯粹的涡旋（Vortex）。
所有流线都是围绕原点的同心圆。
向量方向永远垂直于该点与原点的连线。
这类似于水池排水时形成的漩涡。
StreamPlot(z*1@i,z=[-2-2@i,2+2@i])
7：F(z) = z^2
表达式：F(z) = (x + i*y)^2 = (x^2 - y^2) + i*(2*x*y)
对应的向量场：( u, v ) = ( x^2 - y^2, 2*x*y )
流线图分析：
这个场看起来会比前两个复杂。
沿着 x 轴 (y=0) 和 y 轴 (x=0) 看，场的方向很好判断。
你会看到两个“源”或“汇”挤压在一起形成的复杂图案，流线是双曲线形状。它在原点是一个鞍点（Saddle Point）。
StreamPlot(z^2,z=[-2-2@i,2+2@i])
8：F(z) = 1 / z
表达式：F(z) = 1 / (x + i*y) = (x - i*y) / (x^2 + y^2)
对应的向量场：( u, v ) = ( x / (x^2+y^2), -y / (x^2+y^2) )
流线图分析：
这是源和涡旋的叠加！
实部 x/(r^2) 代表一个源。
虚部 -y/(r^2) 代表一个涡旋。
最终的流线是从原点出发，同时向外且旋转的螺旋线。
这在物理上对应着一个偶极子（Dipole）的场，但带有旋转效应。原点 z=0 是奇点。
StreamPlot(1/z,z=[-2-2@i,2+2@i])
9：F(z) = exp(z) （指数函数）
表达式：利用欧拉公式 exp(x + i*y) = exp(x) * (cos(y) + i*sin(y))
对应的向量场：( u, v ) = ( exp(x)*cos(y), exp(x)*sin(y) )
流线图分析：
这个场在 y 方向是周期性的（因为包含 cos(y) 和 sin(y)）。
沿着任意一条水平线 (y = constant)，场的强度 |F| = exp(x) 随着 x 增大而指数增长。
流线图会呈现出周期性重复的、从左向右发散的图案。
StreamPlot(exp(z),z=[-2-2@i,2+2@i])
根据常见的子表达式重写符号表达式
[r，sigma] = subexpr(expr)根据公共子表达式重写符号表达式expr,用符号变量sigma替换该公共子表达式. 输入表达式expr不能包含变量sigma.
expr —— 包含常见子表达式的长表达式,符号表达,符号功能
r —— 用缩写替换普通子表达式的表达式,符号表达,符号功能
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def subexpression_common_elimination(input_str):
"""
对标 MATLAB 的 subexpr 函数，通过公共子表达式重写符号表达式
Args:
input_str: 可解析为矩阵或符号表达式的字符串
Returns:
成功时返回元组 (替换规则列表, 简化后的表达式)
失败时返回错误描述字符串
Examples:
>>> subexpression_common_elimination("[[x + y, (x + y)**2], [exp(x + y), 1/(x + y)]]")
([(x0, x + y)], Matrix([
[ x0, x0**2],
[exp(x0), 1/x0]]))
>>> subexpression_common_elimination("sin(a) + sqrt(a) + a**2")
([(x0, a)], sin(x0) + sqrt(x0) + x0**2)
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if expr.free_symbols:
R, sigma = sp.cse(expr)
result = R, sigma
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
# 矩阵测试用例
test_cases = [
"sin(a) + sqrt(a) + a^2",
# [sqrt(a) + a**2 + sin(a)]
]
for case in test_cases:
print(f"\n输入: {str(case)[:50]}...")
result = subexpression_common_elimination(str(case)) # 统一转换为字符串
if isinstance(result, tuple):
print("替换规则:", result[0])
print("简化表达式:", result[1])
else:
print("错误输出:", result)
两个子空间之间的角度
theta = subspace(A,B) 计算 A 和 B 的列指定的两个子空间之间的角度. 如果 A 和 B 是单位长度的列向量, 则此角度与 acos(abs(A'*B)) 相同
A, B —— 输入矩阵, 向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.linalg import svd
def subspace_angle(A, B):
"""
计算两个子空间之间的最大主角度（单位：弧度）
参数:
A, B (np.ndarray): 数值矩阵，每一列代表子空间的基向量
返回:
float: 两个子空间之间的夹角（弧度制）
算法说明:
1. 对两个矩阵进行QR分解获得正交基
2. 计算正交基的乘积矩阵的奇异值
3. 最大主角度由最小奇异值的反余弦确定
"""
# 处理空矩阵情况
if A.size == 0 or B.size == 0:
raise ValueError("输入矩阵不能为空")
# 进行QR分解获得正交基（reduced模式）
Q_A, _ = np.linalg.qr(A, mode='reduced')
Q_B, _ = np.linalg.qr(B, mode='reduced')
# 计算正交基的乘积
C = np.dot(Q_A.T, Q_B)
# 计算奇异值（降序排列）
_, s, _ = svd(C)
# 处理数值误差：将大于1的值截断为1
s = np.clip(s, 0.0, 1.0)
# 最大主角度由最小奇异值确定
return np.arccos(s.min()) if s.size > 0 else np.pi / 2
def subspace_arch_matrix(input_str):
"""
计算两个矩阵列空间之间的最大主角度（对标MATLAB subspace函数）
参数:
input_str (str): 输入字符串，格式应为包含两个矩阵的元组
示例：'([[1,0], [0,1]], [[1,1], [1,1]])'
返回:
float: 子空间夹角（弧度）或错误信息字符串
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
# 验证输入格式是否为二元组
if not (isinstance(expr, tuple) and len(expr) == 2):
return f"输入格式错误：需要两个矩阵组成的元组，但收到 {input_str}"
# 转换矩阵
matrix_A = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
matrix_B = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
# 验证矩阵有效性
if matrix_A is None or matrix_B is None:
return "输入包含无效矩阵，请检查格式是否正确（例如：[[1,0],[0,1]]）"
# 转换为数值矩阵
A = np.array(matrix_A.tolist(), dtype=float)
B = np.array(matrix_B.tolist(), dtype=float)
# 验证矩阵维度
if A.shape[0] != B.shape[0]:
return f"矩阵行数不匹配：矩阵A有{A.shape[0]}行，矩阵B有{B.shape[0]}行"
return subspace_angle(A, B)
except Exception as e:
return f"计算错误：{str(e)}"
# --------------------------
# 示范用例
# --------------------------
if __name__ == "__main__":
# 测试用例1：相同子空间
test_case1 = '([[1,0], [0,1]], [[1,0], [0,1]])'
print(f"测试用例1: {test_case1}")
angle = subspace_arch_matrix(test_case1)
print(f"结果: {angle:.4f} 弧度 ({np.degrees(angle):.1f} 度)\n")
# 0.0000 弧度 (0.0 度)
# 测试用例2：正交子空间
test_case2 = '([[1,0], [0,0]], [[0,0], [0,1]])' # A的列空间是x轴，B的列空间是y轴
print(f"测试用例2: {test_case2}")
angle = subspace_arch_matrix(test_case2)
print(f"结果: {angle:.4f} 弧度 ({np.degrees(angle):.1f} 度)\n")
# .0000 弧度 (0.0 度)
# 测试用例3：45度夹角
test_case3 = '([[1], [0]], [[1], [1]])' # A为x轴，B为y=x线
print(f"测试用例3: {test_case3}")
angle = subspace_arch_matrix(test_case3)
print(f"结果: {angle:.4f} 弧度 ({np.degrees(angle):.1f} 度)\n")
# 0.7854 弧度 (45.0 度)
# 测试用例4：三维空间中的斜交子空间
test_case4 = '([[1,0], [0,1], [0,0]], [[1,1], [1,1], [0,0]])' # 列空间平面与另一平面夹角
print(f"测试用例4: {test_case4}")
angle = subspace_arch_matrix(test_case4)
print(f"结果: {angle:.4f} 弧度 ({np.degrees(angle):.1f} 度)\n")
# 输出: 0.0000 弧度 (0.0 度)
三维曲面法向量
vect1, -vect1 = surfnorm(function,varlist,point) 计算并得到在特定点处的曲面的正负法向量,这些向量表示曲面在该点处的垂直方向,分别朝向曲面内外.
这些向量对于理解曲面的几何性质,绘制曲面的法线图以及进行物理模拟都非常有用.
function -- 曲面函数
varlist -- 曲面函数的变量
point -- 曲面上的特定坐标点
vect1, vect2 -- 向量.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def surface_normal_expression(input_str):
"""
计算三维曲面在指定点处的法向量表达式
参数:
input_str (str): 输入字符串，应为格式 ('函数', (变量1,变量2,变量3), (点坐标1,点坐标2,点坐标3)) 的合法 SymPy 表达式
示例: '(x**2 + y**2 + z**2, (x,y,z), (1,0,0))'
返回:
tuple: 包含正负两个法向量的元组，格式为 (正向单位法向量, 负向单位法向量)
或返回错误信息字符串
注意:
1. 法向量通过计算函数梯度后归一化得到
2. 梯度向量指向函数增长最快的方向，即曲面的正向法向量
3. 若输入格式错误或计算过程中出现异常，返回错误描述
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查是否为三元组结构 (函数, 变量元组, 点坐标元组)
if isinstance(expr, tuple) and len(expr) == 3:
# 解析表达式各部分
func = expr[0] # 曲面函数表达式
variables = expr[1] # 变量元组 (如 (x,y,z))
point_values = expr[2] # 求值点坐标元组 (如 (1,0,0))
# 创建变量到点值的映射字典
substitution_dict = {var: val for var, val in zip(variables, point_values)}
# 计算各变量的偏导数（梯度分量）
gradient_components = [sp.diff(func, var) for var in variables]
# 在指定点代入求值梯度向量
gradient_vector = sp.Matrix([comp.subs(substitution_dict) for comp in gradient_components])
# 计算梯度向量的模长
gradient_magnitude = gradient_vector.norm()
# 处理零梯度异常情况
if gradient_magnitude == 0:
return "错误：梯度为零向量，无法计算法向量"
# 计算单位法向量（正向）
unit_normal_vector = gradient_vector / gradient_magnitude
# 返回正负两个方向的法向量
result = (unit_normal_vector, -unit_normal_vector)
else:
error = True
return result if not error else f"输入格式错误: 需要 (函数, 变量元组, 点坐标元组) 格式，但收到 {input_str}"
except Exception as e:
return f"计算错误: {str(e)}"
# --------------------------
# 示范用例
# --------------------------
if __name__ == "__main__":
# 测试用例1：标准球面函数
test_case1 = '(x**2 + y**2 + z**2, (x,y,z), (1,0,0))'
print("测试用例1:", test_case1)
print("结果:", surface_normal_expression(test_case1))
# (Matrix([[1],
# [0],
# [0]]),
# Matrix([[-1],
# [ 0],
# [ 0]]))
# 测试用例2：平面函数
test_case2 = '(2*x + 3*y + 4*z, (x,y,z), (5,-2,1))'
print("\n测试用例2:", test_case2)
print("结果:", surface_normal_expression(test_case2))
# (Matrix([[2*sqrt(29)/29],
# [3*sqrt(29)/29],
# [4*sqrt(29)/29]]),
# Matrix([[-2*sqrt(29)/29],
# [-3*sqrt(29)/29],
# [-4*sqrt(29)/29]]))
# 测试用例3：变量与点值数量不匹配
test_case3 = '(x**2 + y**2, (x,y,z), (1,2))'
print("\n测试用例4:", test_case3)
print("结果:", surface_normal_expression(test_case3))
# (Matrix([[ sqrt(5)/5],
# [2*sqrt(5)/5],
# [ 0]]),
# Matrix([[ -sqrt(5)/5],
# [-2*sqrt(5)/5],
# [ 0]]))
奇异值分解
U, S, V = svd(A)
U - 左奇异向量组成的矩阵
S - 奇异值组成的对角矩阵
V - 右奇异向量组成的矩阵
A - 输入, 矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.linalg import svd as scipy_svd
def singular_values_matrix(input_str):
"""
执行完整的奇异值分解，对标MATLAB的svd函数
参数:
input_str: 字符串形式的矩阵输入
返回:
当输入为数值矩阵时：
U: 左奇异向量矩阵 (m x m)
S: 奇异值对角矩阵 (m x n)
V: 右奇异向量矩阵 (n x n)
当输入包含符号时返回符号解
错误时返回字符串提示
"""
try:
# 符号解析输入
expr = sp.sympify(input_str)
# 输入有效性验证
if isinstance(expr, tuple):
return "错误：不接受元组输入，请输入单个矩阵"
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is None:
return "错误：输入不是有效矩阵"
# 符号矩阵的特殊处理
try:
# 尝试符号计算SVD
U, S, V = A.singular_value_decomposition()
return (U.evalf(), sp.diag(*S).evalf(), V.evalf())
except NotImplementedError:
# 符号方法失败时转数值计算
np_A = np.array(A, dtype=float)
U_np, S_np, Vt_np = scipy_svd(np_A, full_matrices=True)
return (sp.Matrix(U_np).evalf(),
sp.Matrix(np.diag(S_np)).evalf(),
sp.Matrix(Vt_np.T).evalf())
except Exception as e:
return f"分解错误: {str(e)}"
# 示范代码
if __name__ == "__main__":
# 示例1：标准数值矩阵
case1 = "[[1, 2], [3, 4], [5, 6]]"
print("案例1 (3x2数值矩阵):")
result1 = singular_values_matrix(case1)
if isinstance(result1, tuple):
print("U:\n", result1[0])
# Matrix([[0.883461017698525, 0.229847696400071],
# [0.240782492132547, 0.524744818760294],
# [-0.401896033433432, 0.819641941120516]])
print("S:\n", result1[1])
# Matrix([[0.514300580658644, 0, 0, 0],
# [0, 0, 0, 0],
# [0, 0, 0, 0],
# [0, 0, 0, 9.52551809156511]])
print("V:\n", result1[2])
# Matrix([[-0.784894453267052, 0.619629483829340],
# [0.619629483829340, 0.784894453267052]])
else:
print(result1)
# 示例2：符号矩阵测试
case2 = "[[x, 0], [0, y]]"
print("\n案例2 (符号矩阵):")
print(singular_values_matrix(case2))
# (Matrix([[x/(x*conjugate(x))**0.5, 0],
# [ 0, y/(y*conjugate(y))**0.5]]),
# Matrix([[(x*conjugate(x))**0.5, 0, 0, 0],
# [ 0, 0, 0, 0],
# [ 0, 0, 0, 0],
# [ 0, 0, 0, (y*conjugate(y))**0.5]]),
# Matrix([[1.0, 0],
# [ 0, 1.0]]))
添加数据后修改SVD
[U1，S1，V1] = svdappend(U,S,V,D),其中A=U*S*V'是一个现有的奇异值分解(SVD),计算[A D]的SVD,而不显式形成A或[A D].结果相等,四舍五入误差
U,S,V —— 现有SVD因素,矩阵
D —— 要追加的新列或行,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.linalg import svd
def singular_values_append(input_str):
"""
在追加新数据列后更新SVD分解，对标MATLAB的svdappend功能。
参数:
input_str: 字符串形式的元组，包含四个元素：(U, S, V, D)，其中：
- U: 原左奇异向量矩阵
- S: 原奇异值对角矩阵
- V: 原右奇异向量矩阵
- D: 要追加的新数据列（必须与原矩阵列数相同）
返回:
更新后的U_new, S_new, V_new（SymPy矩阵）或错误信息。
"""
try:
expr = sp.sympify(input_str)
error_msg = None
# 输入格式验证
if not isinstance(expr, tuple) or len(expr) != 4:
error_msg = "输入必须为四元组 (U, S, V, D)"
else:
U, S, V, D = expr
U = sp.Matrix(U) if isinstance(U, list) else None
S = sp.Matrix(S) if isinstance(S, list) else None
V = sp.Matrix(V) if isinstance(V, list) else None
D = sp.Matrix(D) if isinstance(D, list) else None
if None in (U, S, V, D):
error_msg = "无效的矩阵输入"
else:
# 转换为NumPy数组进行数值计算
U_np = np.array(U, dtype=float)
S_np = np.array(S, dtype=float)
D_np = np.array(D, dtype=float)
V_np = np.array(V, dtype=float)
# 验证维度一致性
m, k = U_np.shape
s_dim = S_np.shape
if s_dim[0] != s_dim[1] or s_dim[0] != k:
error_msg = "S必须是方阵且与U的列数一致"
elif D_np.shape[0] != m:
error_msg = "D的行数必须与U的行数一致"
elif V_np.shape[1] != k:
error_msg = "V的列数必须与U的列数一致"
if error_msg:
return error_msg
# 增量SVD核心算法
# 步骤1：计算新列在现有列空间的投影
proj = U_np.T @ D_np
# 步骤2：计算正交补分量
D_ortho = D_np - U_np @ proj
# 步骤3：对正交分量进行QR分解
Q, R = np.linalg.qr(D_ortho, mode='reduced')
# 步骤4：构造中间矩阵
K = np.block([[S_np, proj],
[np.zeros((R.shape[0], S_np.shape[1])), R]])
# 步骤5：计算中间矩阵的SVD
Uk, Sk, Vkt = svd(K, full_matrices=False)
# 步骤6：构造新的U和S
U_new = np.hstack((U_np, Q)) @ Uk
S_new = np.diag(Sk)
# 步骤7：构造新的V
V_ext = np.eye(V_np.shape[0] + 1) # 扩展单位矩阵
V_ext[:V_np.shape[0], :V_np.shape[1]] = V_np
V_new = V_ext @ Vkt.T
# 转换回SymPy矩阵并格式化
return (
sp.Matrix(U_new).evalf(chop=True),
sp.Matrix(S_new).evalf(chop=True),
sp.Matrix(V_new).evalf(chop=True)
)
except Exception as e:
return f"错误: {str(e)}"
# 示范代码
if __name__ == "__main__":
# 示例1：有效输入
U = [[1, 0], [0, 1], [0, 0]]
S = [[5, 0], [0, 3]]
V = [[1, 0], [0, 1]]
D = [[2], [3], [4]]
input_str = f"({U}, {S}, {V}, {D})"
result = singular_values_append(input_str)
print("示例1结果:")
if isinstance(result, tuple):
print("U_new:\n", result[0])
# Matrix([[0.707608901615485, 0.700874167360754, 0.0898055893633361],
# [0.493471643807718, -0.581132003549788, 0.647125436996504],
# [0.505742403909574, -0.413595207874290, -0.757075706195349]])
print("S_new:\n", result[1])
# Matrix([[6.23714970083448, 0, 0],
# [0, 4.39365160364841, 0],
# [0, 0, 2.18947235541766]])
print("V_new:\n", result[2])
# Matrix([[0.567253421479376, 0.797598706709882, 0.205085004021905],
# [0.237354401037559, -0.396798874358104, 0.886686834015394],
# [0.788597770810977, -0.454298312084543, -0.414399082418255]])
else:
print(result)
奇异值和向量的子集
s = svds(A) 返回一个向量,其中包含矩阵A的六个最大的奇异值.当使用svd计算所有奇异值的计算量很大时(例如对于大型稀疏矩阵而言),可以使用此函数
s = svds(A,k) 返回k个最大奇异值
A — 输入, 矩阵
k — 要计算的奇异值个数,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.sparse import linalg as sparse_linalg
def singular_values_subset(input_str):
"""
计算矩阵的奇异值子集，对标MATLAB的svds函数。
参数:
input_str: 输入的矩阵或元组（矩阵和奇异值数量k）的字符串表示。
返回:
如果输入有效，返回奇异值的SymPy矩阵（按降序排列）；否则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
U, S, Vt = None, None, None
# 处理元组输入（矩阵和k）
if isinstance(expr, tuple) and len(expr) == 2:
A_expr, k_expr = expr[0], expr[1]
A = sp.Matrix(A_expr) if isinstance(A_expr, list) else None
if A is not None and k_expr.is_integer:
k = int(k_expr)
N_A = np.array(A, dtype=float)
m, n = N_A.shape
max_k = min(m, n) - 1
if max_k <= 0:
max_k = 1
k = min(k, max_k)
U, S, Vt = sparse_linalg.svds(N_A, k)
else:
error = True
# 处理单个矩阵输入，自动确定k
else:
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is not None:
N_A = np.array(A, dtype=float)
m, n = N_A.shape
max_k = min(m, n) - 1
if max_k <= 0:
max_k = 1
k = min(6, max_k)
U, S, Vt = sparse_linalg.svds(N_A, k)
else:
error = True
if error:
return f"输入错误: {input_str}"
else:
# 将奇异值按降序排列
S_sorted = np.sort(S)[::-1]
return sp.Matrix(S_sorted).evalf()
except Exception as e:
return f"错误: {e}"
# 示范代码
if __name__ == "__main__":
# 示例1：输入矩阵和k值
input1 = "(([[1, 2], [3, 4], [5, 6]]), 1)"
result1 = singular_values_subset(input1)
print(f"示例1结果:\n{result1}")
# Matrix([[9.52551809156511]])
# 示例2：输入矩阵，自动确定k
input2 = "[[1, 2], [3, 4]]"
result2 = singular_values_subset(input2)
print(f"示例2结果:\n{result2}")
# Matrix([[5.46498570421904]
计算低秩矩阵草图的SVD
[U,S,V] = svdsketch(A) 返回输入矩阵A的低秩矩阵草图的奇异值分解(SVD).
[U,S,V] = svdsketch(A,tol) 指定矩阵草图的容差。svdsketch 使用 tol 以自适应方式确定矩阵草图逼近的秩。随着容差变大，矩阵草图中使用的 A 的特征越来越少。
矩阵草图是一种低秩逼近,仅反映A的最重要特征(最大容差). 与使用svds相比,它能够更快地计算大型矩阵的部分SVD.
A — 输入, 矩阵, 指定为稀疏矩阵或满矩阵。A 通常为大型稀疏矩阵，但不总是这样。svdsketch 最适合对特征数量相对较少的秩亏矩阵进行操作。
tol — 矩阵草图容差, eps(class(A))^(1/4) (默认) | 实数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from sklearn.utils.extmath import randomized_svd
def sklearn_randomized_svd(input_str):
"""
对标MATLAB的svdsketch函数，计算低秩矩阵近似SVD
参数:
A: 输入矩阵，支持SymPy矩阵/列表/NumPy数组
tol: 近似容差 (默认1e-10)
max_rank: 最大允许秩 (可选)
random_state: 随机种子 (可选)
返回:
U: 左奇异向量矩阵
S: 奇异值对角矩阵
V: 右奇异向量矩阵
rank: 实际使用的秩
示例:
>>> A = np.random.rand(100, 50)
>>> U, S, V, r = svdsketch(A, tol=1e-5)
"""
try:
# 定义容差，用于后续判断奇异值的截断条件
tol = 1e-10
# 最大秩，初始设为 None，后续会根据矩阵情况确定
max_rank = None
# 随机状态，用于随机化过程的可重复性，初始设为 None
random_state = None
# 错误标志，用于标记是否出现错误情况
error = False
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 检查 expr 是否为元组
if isinstance(expr, tuple):
# 如果元组长度为 4，分别将元素赋值给 A, tol, max_rank, random_state
if len(expr) == 4:
A, tol, max_rank, random_state = expr
# 如果元组长度为 3，分别将元素赋值给 A, tol, max_rank
elif len(expr) == 3:
A, tol, max_rank = expr
# 如果元组长度为 2，分别将元素赋值给 A, tol
elif len(expr) == 2:
A, tol = expr
# 若元组长度不符合上述情况，标记错误
else:
error = True
# 若 expr 不是元组，直接将其赋值给 A
else:
A = expr
# 检查 A 是否为有效的矩阵，如果不是则返回 None
sp_A = sp.Matrix(A) if isinstance(A, list) else None
if sp_A is None:
raise ValueError("无效的矩阵输入")
elif sp_A.is_symbolic():
raise ValueError("需要数值矩阵输入")
else:
np_A = np.array(sp_A, dtype=float)
# 矩阵维度检查
if np_A.ndim != 2:
raise ValueError("输入必须是二维矩阵")
# 自动确定目标秩
m, n = np_A.shape
if max_rank is None:
max_rank = min(m, n) - 1
# 使用自适应随机SVD
rank = min(max_rank, min(m, n))
U, S, Vt = randomized_svd(np_A, n_components=rank,
random_state=random_state)
# 根据容差调整实际秩
sigma_norm = np.linalg.norm(S) # 奇异值的F范数
if sigma_norm > 0:
cumulative = np.cumsum(S[::-1] ** 2)[::-1] # 降序累积平方和
for k in range(len(cumulative)):
if np.sqrt(cumulative[k]) <= tol * sigma_norm:
rank = k
break
rank = max(1, min(rank, max_rank))
# 截断到实际秩
U = U[:, :rank]
S = S[:rank]
Vt = Vt[:rank, :]
if not error:
return U, np.diag(S), Vt.T, rank
except Exception as e:
print(f"计算错误: {e}")
return None, None, None, 0
# 示例代码
if __name__ == "__main__":
# 示例3：数值矩阵测试
try:
C = "[[2, 0], [0, 1]]"
print(sklearn_randomized_svd(C))
# (array([[ 1.00000000e+00],
# [-2.16840434e-19]]),
# array([[2.]]),
# array([[1.],
# [0.]]),
# 1)
except Exception as e:
print(f"示例3 错误捕获: {e}")
求解关于X的西尔维斯特方程AX+XB=C
X=sylvester(A,B,C)将解X返回到sylvester方程.
输入A是一个m乘m矩阵,输入B是一个n乘n矩阵,C和X都是m乘n矩阵.
A,B,C — 输入矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy import linalg as dense_linalg
def solve_sylvester_equation(A, B, C):
"""
解西尔维斯特方程 AX + XB = C。
参数:
A: 方程中的矩阵A，支持 SymPy 矩阵、列表或 NumPy 数组。
B: 方程中的矩阵B，支持类型同上。
C: 方程中的矩阵C，支持类型同上。
返回:
SymPy 矩阵: 解矩阵X（若成功）。
str: 错误信息（若输入无效或计算失败）。
示例:
>>> A = [[1, 0], [0, 2]]
>>> B = [[3, 1], [1, 4]]
>>> C = [[5, 6], [7, 8]]
>>> solve_sylvester_equation(A, B, C)
Matrix([
[ 1.0, 1.0],
[-1.5, 2.0]])
"""
try:
# 转换为 SymPy 矩阵
A_mat = sp.Matrix(A) if isinstance(A, list) else None
B_mat = sp.Matrix(B) if isinstance(B, list) else None
C_mat = sp.Matrix(C) if isinstance(C, list) else None
if None in (A_mat, B_mat, C_mat):
return "错误: 输入无法转换为有效矩阵"
# 符号检查
if A_mat.is_symbolic() or B_mat.is_symbolic() or C_mat.is_symbolic():
return "错误: 暂不支持符号运算"
# 维度验证
if not A_mat.is_square or not B_mat.is_square:
return "错误: A 和 B 必须是方阵"
if C_mat.rows != A_mat.rows or C_mat.cols != B_mat.rows:
return f"错误: C的维度应为({A_mat.rows}, {B_mat.rows})，实际为({C_mat.rows}, {C_mat.cols})"
# 转换为 NumPy 数组
A_np = np.array(A_mat, dtype=float)
B_np = np.array(B_mat, dtype=float)
C_np = np.array(C_mat, dtype=float)
# 解方程
X_np = dense_linalg.solve_sylvester(A_np, B_np, C_np)
return sp.Matrix(X_np)
except ValueError as ve:
return f"数值错误: {ve}"
except Exception as e:
return f"未知错误: {e}"
if __name__ == "__main__":
# 示例 1: 标准数值解
A = [[1, 0], [0, 2]]
B = [[3, 1], [1, 4]]
C = [[5, 6], [7, 8]]
print("示例1解:\n", solve_sylvester_equation(A, B, C))
# Matrix([[1.00000000000000, 1.00000000000000],
# [1.17241379310345, 1.13793103448276]])
对称近似最小度置换
对称近似最小度置换,对于对称正定矩阵A,命令p=symamd(S)返回置换向量p,因此S(p,p)倾向于比S具有更稀疏的乔列斯基因子.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import reverse_cuthill_mckee
def symamd_permutation_vector(input_str):
"""
对标 MATLAB 的 symamd 函数（注意：实际使用 RCM 算法实现带宽优化）。
该函数接受矩阵输入，返回行/列置换向量以优化稀疏性模式。
参数:
input_str: 矩阵的字符串表示（如 "[[1, 2], [3, 4]]"）或 SymPy 矩阵。
返回:
置换向量 (numpy数组) 或错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, list):
matrix = sp.Matrix(expr)
# 检查矩阵是否包含符号
if matrix.is_symbolic():
return "错误: 矩阵包含符号变量，暂不支持符号运算"
# 转换为 numpy 数组并创建稀疏矩阵
try:
np_matrix = np.array(matrix, dtype=float)
graph = csr_matrix(np_matrix)
# 使用反向 Cuthill-McKee 算法获取置换
# 注意：这里实际使用的是 RCM 算法而非精确的 AMD 算法
result = reverse_cuthill_mckee(graph)
except ValueError as ve:
return f"数值转换错误: {ve}"
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示范代码
if __name__ == "__main__":
# 示例 1: 数值矩阵
matrix_str = "[[1, 0, 4], [0, 2, 5], [4, 5, 3]]"
permutation = symamd_permutation_vector(matrix_str)
print(f"示例1置换向量: {permutation}")
# [1 2 0]
# 示例 2: 非对称矩阵
non_sym_matrix = "[[1, 2], [3, 4]]"
print(symamd_permutation_vector(non_sym_matrix))
# [1 0]
级数的乘积
F = symprod(f,k,a,b) 返回表达式 f 指定的级数的乘积, 该级数取决于符号变量 k. k 的值范围从 a 到 b. 如果您未指定 k，symprod 将使用 symvar 确定的变量. 如果f是常量, 则默认变量为x.
F = symprod(f,k) 返回表达式 f 指定的级数的乘积, 该级数取决于符号变量 k. k 的值从 1 开始, 上限未指定. 乘积 F 以 k 的形式返回, 其中 k 表示上限. 此乘积 F 与不定乘积不同. 如果您未指定 k, symprod 将使用 symvar 确定的变量. 如果 f 是常量, 则默认变量为 x.
f — 定义级数项的表达式, 符号表达式, 符号函数, 符号向量, 符号矩阵
k — 乘积指数,符号变量
a — 乘积指数的下限, 数, 符号数, 符号变量, 符号表达式, 符号函数
b — 乘积指数的上限, 数, 符号数, 符号变量, 符号表达式, 符号函数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def symvar_product_series(input_str):
"""
计算符号乘积级数，类似 MATLAB 的 symprod 函数
参数:
input_str: 输入字符串，支持以下格式：
- 普通符号表达式 (例如 "k")
- 矩阵的嵌套列表 (例如 "[[k, 1], [2, 3]]")
- 元组格式，可包含：
* (表达式, 变量) (例如 "(k, k)")
* (表达式, 变量, 下限) (例如 "(k, k, 1)")
* (表达式, 变量, 下限, 上限) (例如 "(k, k, 1, 5)")
* (矩阵, 变量, 下限, 上限) (例如 "(Matrix([[k]]), k, 1, n)")
返回:
SymPy 表达式/矩阵 或 错误消息
示例:
>>> symvar_product_series("k")
Product(k, (k, 1, n))
>>> symvar_product_series("(k, k, 1, 5)")
120
>>> symvar_product_series("([[k, 1], [2, 3]], k, 1, n)")
Matrix([
[Product(k, (k, 1, n)), 1],
[ 2**n, 3**n]])
"""
try:
# 尝试解析输入字符串为Python对象
expr = sp.sympify(input_str)
error = False
result = None
def evaluation_product_series(f_expr, k=None, a=None, b=None):
"""
内部函数：计算乘积级数
参数:
f_expr: 乘积表达式
k: 乘积变量
a: 下限，默认为1
b: 上限，默认为符号n
"""
# 自动确定乘积变量
if f_expr.free_symbols and (k is None):
k = list(f_expr.free_symbols)[0] # 默认取第一个符号
elif k is None:
k = sp.symbols('x') # 默认变量
# 设置默认上下限
a = 1 if a is None else a
b = sp.symbols('n') if b is None else b
return sp.product(f_expr, (k, a, b)).evalf()
# 处理元组输入的情况
if isinstance(expr, tuple):
# 矩阵处理分支 (格式: (矩阵, 变量, 下限, 上限))
if len(expr) >= 4 and isinstance(expr[0], list):
F = sp.Matrix(expr[0])
k_var = expr[1]
lower = expr[2]
upper = expr[3]
# 对矩阵每个元素计算乘积
result = sp.Matrix(F.shape[0], F.shape[1],
lambda i, j: evaluation_product_series(
F[i, j], k_var, lower, upper))
# 普通表达式处理分支 (格式: (表达式, 变量, 下限, 上限))
elif len(expr) >= 4 and expr[0].free_symbols:
f_expr = expr[0]
k_var = expr[1]
lower = expr[2]
upper = expr[3]
result = evaluation_product_series(f_expr, k_var, lower, upper)
# 简写格式处理 (例如 (矩阵, 变量))
elif isinstance(expr[0], list):
F = sp.Matrix(expr[0])
k_var = expr[1]
result = sp.Matrix(F.shape[0], F.shape[1],
lambda i, j: evaluation_product_series(F[i, j], k_var))
# 简写格式处理 (例如 (表达式, 变量))
elif expr[0].free_symbols:
f_expr = expr[0]
k_var = expr[1]
result = evaluation_product_series(f_expr, k_var)
else:
error = True
# 处理矩阵输入
elif isinstance(expr, list):
F = sp.Matrix(expr)
result = sp.Matrix(F.shape[0], F.shape[1],
lambda i, j: evaluation_product_series(F[i, j]))
# 处理普通符号表达式
elif expr.free_symbols:
result = evaluation_product_series(expr)
# 处理纯数值输入
elif expr.is_number:
result = evaluation_product_series(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(symvar_product_series("k"))
# factorial(n)
# 具体数值范围乘积
print(symvar_product_series("(k, k, 1, 5)"))
# 120.000000000000
# 矩阵乘积
mat_expr = symvar_product_series("([[k, 1], [2, 3]], k, 1, n)")
print(mat_expr)
# Matrix([[factorial(n), 1.00000000000000],
# [2.0**n, 3.0**n]])
符号和(有限或无限)
F = symsum(f,k,a,b)从下界a到上界b,返回级数F相对于求和索引k的符号定和.如果不指定k，symsum将使用symvar确定的变量作为求和索引.如果f是常数，那么默认变量是x.
F = symsum(f,k)返回级数F相对于求和索引k的不定和(反差).f自变量定义级数,使得不定和f满足关系式f(k+1)-f(k)=f(k).如果不指定k,symsum将使用symvar确定的变量作为求和索引.如果f是常数,那么默认变量是x。
f —— 定义级数项的表达式,符号表达式,符号函数,符号向量,符号矩阵,符号数
k —— 求和指数,符号变量
a —— 求和指数的下界,数字,符号数,符号变量,符号表达式,符号函数
b —— 求和指数的上界,数字,符号数,符号变量,符号表达式,符号函数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def symsum(input_str):
"""
对标 MATLAB 的 symsum 函数，计算符号级数和
参数:
input_str: 输入字符串，支持以下格式：
- "f" : 自动选择变量，从1到无穷大
- "(f, k)" : 从k=1到无穷大
- "(f, k, a)" : 从k=a到无穷大
- "(f, k, a, b)" : 从k=a到k=b
- "[[f1, f2], [f3, f4]]" : 矩阵输入
返回:
- 符号表达式结果
- SymPy 矩阵（矩阵输入时）
- 错误信息字符串
示例:
>>> symsum("(k, k, 1, 5)") # Σk from 1到5 → 15
>>> symsum("(1/k^2, k, 1, oo)") # π²/6
>>> symsum("[[k, k^2], [k^3, 1]]") # 矩阵各元素求和
"""
try:
# 安全解析输入字符串
expr = sp.sympify(input_str)
def get_summation(f, k=None, a=1, b=sp.oo):
"""
计算符号级数和
:param f: 被求和的表达式
:param k: 求和变量
:param a: 求和下限
:param b: 求和上限
:return: 符号级数和的结果
"""
if k is None:
if f.free_symbols:
k = tuple(f.free_symbols)[0]
else:
k = sp.symbols('x')
return sp.summation(f, (k, a, b))
if isinstance(expr, tuple):
# 处理元组输入
f = expr[0]
k = expr[1] if len(expr) > 1 else None
a = expr[2] if len(expr) > 2 else 1
b = expr[3] if len(expr) > 3 else sp.oo
if isinstance(f, list):
F = sp.Matrix(f)
# 处理矩阵输入
return sp.Matrix(F.shape[0], F.shape[1], lambda i, j: get_summation(F[i, j], k, a, b))
else:
return get_summation(f, k, a, b)
elif isinstance(expr, list):
F = sp.Matrix(expr)
# 处理矩阵输入，默认从1到无穷大
return sp.Matrix(F.shape[0], F.shape[1], lambda i, j: get_summation(F[i, j]))
else:
# 处理普通表达式输入，默认从1到无穷大
return get_summation(expr)
except SyntaxError:
return "语法错误: 输入格式不正确"
except ValueError as ve:
return f"数值错误: {str(ve)}"
except Exception as e:
return f"未处理的错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试1: 有限求和
print("Σk from 1到5:", symsum("(k, k, 1, 5)"))
# 15
# 测试2: 无限级数
print("Σ1/k² from 1到∞:", symsum("(1/k**2, k, 1, oo)"))
# pi**2/6
# 测试3: 自动变量选择
print("Σm from 1到∞:", symsum("m"))
# oo
# 测试4: 矩阵输入
print("矩阵求和结果:\n", symsum("[[k, k**2], [k**3, 1/k]]"))
# Matrix([[oo, oo],
# [oo, oo]])