隐式导数
iDiff(eq,y,x,n=1)返回dy/dx,并假设eq=0
y — 因变量或因变量列表(以y开头)
x — 取导数的变量
n — 导数的顺序（默认为1）
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def get_result_value(symbols_set, points=None, result=None):
"""
将符号表达式代入具体数值点求值。
参数:
symbols_set (list): 符号列表，按字母顺序排列
points (list): 数值点列表，顺序需与symbols_set中的符号对应
result (sympy.Expr): 待代入的表达式
返回:
sympy.Expr 或 float: 代入后的数值结果或符号表达式
示例:
symbols = [x, y], points = [2, 3] 表示代入x=2, y=3
"""
if points is not None and result is not None:
if len(points) != len(symbols_set):
raise ValueError(f"提供的点数（{len(points)}）与符号数量（{len(symbols_set)}）不匹配。")
subs_dict = {symbol: value for symbol, value in zip(symbols_set, points)}
try:
# 尝试数值求值，保留浮点精度
result = result.subs(subs_dict).evalf()
except Exception as e:
# 非数值情况直接代入符号
result = result.subs(subs_dict)
return result
def implicit_diff_equation(input_str):
"""
处理隐函数求导请求，对标MathStudio的iDiff功能
参数:
input_str (str): 输入字符串，格式为Sympy可解析的元组形式：
(方程, 因变量, 自变量, 阶数(可选), 数值点(可选))
示例: "(Eq(x**2 + y**2, 1), y, x, 2, [0, 1])"
返回:
sympy.Expr 或 str: 求导结果或错误信息
"""
try:
# 预处理输入字符串
processed_str = input_str
error = False
result = None
def remove_last_element(tup):
"""移除元组最后一个元素（用于处理数值点）"""
return tup[:-1] if len(tup) > 1 else tup
# 解析输入字符串
if '==' in processed_str:
# 方程处理：转换为左式-右式=0的形式
expr = sp.sympify(processed_str, evaluate=False)
transformed_equation = f"{expr[0].lhs} - {expr[0].rhs}"
expr_func = sp.sympify(transformed_equation)
else:
expr = sp.sympify(processed_str)
expr_func = expr[0] if isinstance(expr, tuple) else expr
# 处理元组形式的参数
if isinstance(expr, tuple):
# 检查是否包含数值点参数
if len(expr) > 0 and isinstance(expr[-1], list) and all(item.is_number for item in expr[-1]):
points = expr[-1]
expr = remove_last_element(expr)
else:
points = None
# 解析求导参数
if len(expr) >= 3:
equation = expr_func
dep_var = expr[1] # 因变量
indep_var = expr[2] # 自变量
# 处理求导阶数
order = expr[3] if len(expr) > 3 and expr[3].is_number else 1
# 执行隐式求导
result = sp.idiff(equation, dep_var, indep_var, order).simplify()
# 处理符号代入
symbols = sorted(result.free_symbols, key=lambda s: str(s)) # 按字母顺序排序符号
result = get_result_value(symbols, points, result)
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例1：计算 x² + y² = 1 的一阶导数
x, y = sp.symbols('x y')
input1 = "(x**2 + y**2==1, y, x)"
print("示例1结果:", implicit_diff_equation(input1))
# 输出: -x/y
# 示例2：计算二阶导数并代入点(0, 1)
input2 = "(x**2+y**2==1, y, x, 2, [0, 1])"
print("示例2结果:", implicit_diff_equation(input2))
# 输出 -1.0
# 示例3：带数值点的三阶导数
input3 = "(y**3 + x, y, x, 3)"
print("示例3结果:", implicit_diff_equation(input3))
# 输出 -10/(27*y**8)
带有舍入选项的整除
C = idivide(A,B) 将A的每个元素除以B的对应元素,朝零方向舍入到最接近的整数.A和B必须包含实数,并且其中有至少一个必须属于整数类.
如果A和B是数组,则它们必须属于同一整数类,并且大小兼容.
如果A或B是双精度标量,则另一个输入必须为整数类,但不能是int64或uint64.然后,idivide函数以相同的整数类形式返回C.
C = idivide(A,B,opt) 指定替代舍入选项：'fix'、'floor'、'ceil'.例如,idivide(A,B,'ceil') 将商朝正无穷方向舍入到最接近的整数.默认舍入选项是 'fix'.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def integer_division(input_str):
"""
执行整数除法运算，对标MATLAB的idivide函数
参数:
input_str (str): 输入字符串，支持以下格式：
- 标量运算: "(被除数, 除数, '舍入选项')"
- 矩阵运算: "(矩阵A, 矩阵B, '舍入选项')"
- 单矩阵运算: "(矩阵, 标量, '舍入选项')"
返回:
sympy.Expr/sympy.Matrix/str: 计算结果或错误信息
示例:
>>> integer_division("(10, 3, 'fix')") # 10 ÷ 3 向零取整
3
>>> integer_division("([[1,2],[3,4]], 2, 'floor')") # 矩阵逐元素除
Matrix([[0, 1], [1, 2]])
"""
try:
# 预处理输入
expr = sp.sympify(input_str)
error = False
result = None
def elementwise_divide(a, b, option='fix'):
"""
执行逐元素除法，支持不同维度的数值标量、向量和多维数组
参数:
a, b: 数值标量、列表(向量)、嵌套列表(多维数组)或NumPy数组
option: 舍入选项
'fix' - 向零取整 (默认)
'floor' - 向下取整
'ceil' - 向上取整
返回:
与输入维度匹配的结果(标量、列表或NumPy数组)
"""
# 将输入转换为NumPy数组
a_arr = np.asarray(a)
b_arr = np.asarray(b)
# 检查除数是否包含零
if np.any(b_arr == 0):
raise ZeroDivisionError("除数不能包含零")
# 执行逐元素除法
result = a_arr / b_arr
# 应用舍入选项
if option == 'fix':
# 向零取整
result = np.trunc(result)
elif option == 'floor':
# 向下取整
result = np.floor(result)
elif option == 'ceil':
# 向上取整
result = np.ceil(result)
else:
raise ValueError(f"无效舍入选项: {option}")
# 转换为整数类型
result = result.astype(int)
# 保持与输入相同的类型
if isinstance(a, (int, float)) and isinstance(b, (int, float)):
return int(result) # 标量输入返回整数
elif isinstance(a, list) or isinstance(b, list):
return result.tolist() # 列表输入返回列表
return result # NumPy数组输入返回数组
if isinstance(expr, tuple):
# 解析参数数量和类型
n = float(expr[0])
a = float(expr[1])
option = expr[2] if len(expr) > 2 else 'fix' # 默认舍入方式
result = elementwise_divide(n, a, str(option))
else:
error = True
return result if not error else f"输入错误: {input_str}"
except ZeroDivisionError:
return "错误: 除零错误"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print("-10 ÷ 3 (floor):", integer_division("(-10, 3, 'floor')"))
# -4
print("7 ÷ 3 (ceil):", integer_division("(7, 3, 'ceil')"))
# 3
快速傅里叶逆变换
X = ifft(Y) 使用快速傅里叶变换算法计算 Y 的逆离散傅里叶变换。X 与 Y 的大小相同。
如果 Y 是向量，则 ifft(Y) 返回该向量的逆变换。
如果 Y 是矩阵，则 ifft(Y) 返回该矩阵每一列的逆变换。
X = ifft(Y,n) 通过用尾随零填充 Y 以达到长度 n，返回 Y 的 n 点傅里叶逆变换。
X = ifft(Y,n,dim) 返回沿维度 dim 的傅里叶逆变换。例如，如果 Y 是矩阵，则 ifft(Y,n,2) 返回每一行的 n 点逆变换。
Y — 输入数组, 向量 | 矩阵
n — 逆变换长度, 0 (默认) | 非负整数标量
dim — 沿其运算的维度, 正整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def inverse_fft(input_str):
"""
模拟MATLAB的ifft函数行为
参数:
x: 输入数组
n: FFT点数（默认为输入长度）
axis: 计算轴（MATLAB默认按列，对应axis=0）
返回:
与MATLAB ifft结果一致的复数数组
"""
# 检验输入 x 是否为数组
try:
expr = sp.sympify(input_str)
error = False
result = None
n = 0
axis = 0
if isinstance(expr, tuple):
if len(expr) == 3:
matrix_sp = sp.Matrix(expr[0])
n = expr[1]
axis = expr[2]
elif len(expr) == 2:
matrix_sp = sp.Matrix(expr[0])
n = expr[1]
else:
error = True
else:
matrix_sp = sp.Matrix(expr)
# 校验矩阵有效性
if matrix_sp is None:
raise ValueError("输入无法转换为有效矩阵（支持 SymPy 矩阵或列表）")
# ---------------------- 转换为 NumPy 数组 ----------------------
# 保持原始矩阵的复数类型（SymPy 矩阵元素可能是复数）
x = np.array(matrix_sp, dtype=complex)
# 若输入是列矩阵，转换为二维数组（与 MATLAB 列优先一致）
if x.ndim == 1:
x = x.reshape(-1, 1)
# 检验输入 n 是否为有效的整数
if n != 0:
if not isinstance(n, (int, sp.Integer)) or n <= 0:
raise ValueError("参数 n 必须是正整数。")
else:
n = None
# 检验输入 axis 是否为有效的轴索引
if not isinstance(axis, (int, sp.Integer)):
raise ValueError("参数 axis 必须是整数。")
if axis < -x.ndim or axis >= x.ndim:
raise IndexError(f"参数 axis 超出了数组的有效轴范围 [-{x.ndim}, {x.ndim - 1}]。")
# 调整轴方向（MATLAB默认按列计算，对应NumPy的axis=0）
if not error:
result = sp.Matrix(np.fft.ifft(x, n=n, axis=axis))
if result.shape[1] == 1:
result = result.T
return result
except Exception as e:
return f"错误: {e}"
if __name__ == "__main__":
# 一维数组示例
result_1d = inverse_fft("[1, 1j, -1, -1j]")
print("一维数组逆傅里叶变换结果:", result_1d)
# Matrix([[0, 0, 0, 1.00000000000000]])
# 二维数组示例
# 沿第 0 轴进行逆傅里叶变换
result_2d_axis0 = inverse_fft("[[1, 1j], [-1, -1j]], 0, 0")
print("二维数组沿第 0 轴逆傅里叶变换结果:")
print(result_2d_axis0)
# Matrix([[0, 0],
[1.00000000000000, 1.0*I]])
# 沿第 1 轴进行逆傅里叶变换
result_2d_axis1 = inverse_fft("[[1, 1j], [-1, -1j]], 0, 1")
print("二维数组沿第 1 轴逆傅里叶变换结果:")
print(result_2d_axis1)
# Matrix([[0.5 + 0.5*I, 0.5 - 0.5*I],
[-0.5 - 0.5*I, -0.5 + 0.5*I]])
# 指定 FFT 点数示例
result_1d_n = inverse_fft("[[1, 1j], [-1, -1j]], 8")
print("\n一维数组指定 FFT 点数为 8 的逆傅里叶变换结果:", result_1d_n)
# Matrix([[0, 0],
# [0.0366116523516816 - 0.0883883476483184*I, 0.0883883476483184 + 0.0366116523516816*I],
# [0.125 - 0.125*I, 0.125 + 0.125*I],
# [0.213388347648318 - 0.0883883476483184*I, 0.0883883476483184 + 0.213388347648318*I],
# [0.250000000000000, 0.25*I],
# [0.213388347648318 + 0.0883883476483184*I, -0.0883883476483184 + 0.213388347648318*I],
# [0.125 + 0.125*I, -0.125 + 0.125*I],
# [0.0366116523516816 + 0.0883883476483184*I, -0.0883883476483184 + 0.0366116523516816*I]])
二维快速傅里叶逆变换
X = ifft2(Y) 使用快速傅里叶变换算法返回矩阵的二维离散傅里叶逆变换。如果 Y 是一个多维数组，则 ifft2 计算大于 2 的每个维度的二维逆变换。输出 X 的大小与 Y 相同。
X = ifft2(Y,m,n) 在计算逆变换之前截断 Y 或用尾随零填充 Y，以形成 m×n 矩阵。X 也是 m×n。如果 Y 是一个多维数组，ifft2 将根据 m 和 n 决定 Y 的前两个维度的形状。
Y — 输入数组, 矩阵
m — 逆变换行数, 正整数标量
n — 逆变换列数, 正整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def inverse_fft2(input_str):
"""
对标 MATLAB 的 ifft2(Y, m, n) 函数，实现二维逆快速傅里叶变换
参数:
Y: 输入的二维频域数组（实数/复数）
m: 变换后的行数（默认：Y 的行数）
n: 变换后的列数（默认：Y 的列数）
返回:
二维时域数组，与 MATLAB ifft2 结果一致（数值精度误差在 1e-10 以内）
"""
# Y, m=None, n=None
# ---------------------- 参数校验 ----------------------
# 校验 Y 是二维数组
expr = sp.sympify(input_str)
error = False
result = m = n = None
matrix = None
if isinstance(expr, tuple):
matrix = sp.Matrix(expr[0])
if len(expr) == 3:
m, n = expr[1], expr[2]
elif len(expr) == 2:
m = expr[1]
else:
error = True
else:
matrix = sp.Matrix(expr)
if matrix is None:
raise ValueError("输入 Y 必须是二维 NumPy 数组")
else:
Y = np.array(matrix)
# 校验 m 和 n 为正整数或 None
original_rows, original_cols = Y.shape
if m is not None:
if not isinstance(m, sp.Integer) or m <= 0:
raise ValueError(f"参数 m 必须是正整数（当前值：{m}）")
else:
m = original_rows # 默认使用 Y 的行数
if n is not None:
if not isinstance(n, sp.Integer) or n <= 0:
raise ValueError(f"参数 n 必须是正整数（当前值：{n}）")
else:
n = original_cols # 默认使用 Y 的列数
# ---------------------- 补零/截断逻辑 ----------------------
# 若 m/n 大于原尺寸，补零；若小于，截断
# 注意：MATLAB 补零在末尾，截断从开头保留
if not error:
rows_needed = m
cols_needed = n
Y_padded = np.zeros((rows_needed, cols_needed), dtype=Y.dtype)
# 计算实际截断的行数和列数（不超过原尺寸）
truncate_rows = min(original_rows, rows_needed)
truncate_cols = min(original_cols, cols_needed)
# 填充数据（截断或复制原数据）
Y_padded[:truncate_rows, :truncate_cols] = Y[:truncate_rows, :truncate_cols]
# ---------------------- 计算二维逆 FFT ----------------------
# numpy.fft.ifft2 会自动除以变换点数（m*n），与 MATLAB 一致
result = np.fft.ifft2(Y_padded, s=(m, n))
result = sp.Matrix(result)
return result
# ---------------------- 示例代码 ----------------------
if __name__ == "__main__":
# 示例 1：默认参数（m=原行数，n=原列数）
Y1 = np.array([[1, 1j], [-1, -1j]], dtype=complex)
result1 = inverse_fft2("[[1, 1j], [-1, -1j]]")
print("示例1（默认参数）结果：\n", result1)
# Matrix([[0, 0],
# [0.5 + 0.5*I, 0.5 - 0.5*I]])
# 示例 2：补零（m=4, n=4）
Y2 = np.array([[1, 1j], [-1, -1j]], dtype=complex)
result2 = inverse_fft2("[[1, 1j], [-1, -1j]],4,4")
print("\n示例2（补零至4x4)：", result2)
# Matrix([[0, 0, 0, 0],
# [0.125000000000000, 0, -0.125*I, 0.125 - 0.125*I],
# [0.125 + 0.125*I, 0, 0.125 - 0.125*I, 0.250000000000000],
# [0.125*I, 0, 0.125000000000000, 0.125 + 0.125*I]])
# 示例3：截断（m=1, n=1）
Y3 = np.array([[1, 1j], [-1, -1j]], dtype=complex)
result3 = inverse_fft2("[[1, 1j], [-1, -1j]],1,1")
print("\n示例3（截断至1x1）结果：\n", result3)
# Matrix([[1.00000000000000]])
多维快速傅里叶逆变换
X = ifftn(Y) 使用快速傅里叶变换算法返回 N 维数组的多维离散傅里叶逆变换。N 维逆变换相当于沿 Y 的每个维度计算一维逆变换。输出 X 的大小与 Y 相同。
X = ifftn(Y,sz) 将在进行逆变换之前根据向量 sz 的元素截断 Y 或用尾随零填充 Y。sz 的每个元素定义对应变换维度的长度。例如，如果 Y 是一个 5×5×5 数组，X = ifftn(Y,[8 8 8]) 将用零填充每个维度，从而产生 8×8×8 逆变换 X。
Y — 输入数组, 向量 | 矩阵 | 多维数组
sz — 逆变换维度的长度, 正整数向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def inverse_fftn(input_str):
"""
对标MATLAB的ifftn函数，执行多维逆快速傅里叶变换。
参数:
input_str (str): 输入的字符串，格式为"Y"或"(Y, s1, s2, ...)"，
其中Y是数组，s1, s2等是各维度的大小。
返回:
numpy.ndarray: 逆FFT结果。若指定sizes，则形状为sizes；否则与Y相同。
或返回错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
arr = None
sizes = None
# 解析输入，判断是否为元组（包含数组和sizes参数）
if isinstance(expr, tuple):
if len(expr) < 1:
raise ValueError("输入格式错误，应为 (Y, s1, s2, ...)")
arr = expr[0]
sizes = list(expr[1:]) # 提取sizes参数
else:
arr = expr
if arr is None:
raise ValueError("输入 Y 必须是多维数组")
if isinstance(arr, list):
Y = np.array(arr)
original_shape = Y.shape
# 处理sizes参数
if sizes is not None:
# 补全sizes到原数组的维度数，不足部分使用原数组对应维度的大小
while len(sizes) < len(original_shape):
sizes.append(original_shape[len(sizes)])
# 校验并转换sizes为整数
for i in range(len(sizes)):
s = sizes[i]
if isinstance(s, sp.Integer):
s = int(s)
elif isinstance(s, int):
pass
else:
raise ValueError(f"参数 sizes[{i}] 必须是整数，当前类型为 {type(s)}")
if s <= 0:
raise ValueError(f"参数 sizes[{i}] 必须是正整数，当前值为 {s}")
sizes[i] = s
else:
sizes = list(original_shape)
# 创建补零/截断后的数组
padded_Y = np.zeros(sizes, dtype=Y.dtype)
# 计算各维度的切片范围
slices = tuple(slice(0, min(orig, new)) for orig, new in zip(original_shape, sizes))
padded_Y[slices] = Y[slices]
# 执行逆FFTn
result = np.fft.ifftn(padded_Y)
return sp.Array(result)
except Exception as e:
return f"错误: {e}"
# 示范代码
if __name__ == "__main__":
# 示例1：不指定sizes，结果形状与原数组相同
input1 = "[[1, 2], [3, 4]]"
result1 = inverse_fftn(input1)
print("示例1结果:", result1 if not isinstance(result1, str) else result1)
# [[2.5, -0.5],
# [-1.0, 0]]
# 示例2：指定sizes大于原数组维度
input2 = "([[1, 2], [3, 4]], 3, 3)"
result2 = inverse_fftn(input2)
print("示例2结果:", result2 if not isinstance(result2, str) else result2)
# [[1.11111111111111, 0.111111111111111 + 0.577350269189626*I, 0.111111111111111 - 0.577350269189626*I],
# [-0.0555555555555555 + 0.673575314054563*I, -0.388888888888889 + 0.0962250448649376*I, 0.277777777777778 + 0.0962250448649376*I],
# [-0.0555555555555555 - 0.673575314054563*I, 0.277777777777778 - 0.0962250448649376*I, -0.388888888888889 - 0.0962250448649376*I]]
# 示例3：三维数组，sizes部分指定
input3 = "([[[1, 2], [3, 4]], [[5, 6], [7, 8]]], 2, 2)"
result3 = inverse_fftn(input3)
print("示例3结果:", result3 if not isinstance(result3, str) else result3)
# [[[4.5, -0.5],
# [-1.0, 0]],
# [[-2.0, 0],
# [0, 0]]]
逆零频平移
X = ifftshift(Y) 将进行过零频平移的傅里叶变换 Y 重新排列回原始变换输出的样子。换言之，ifftshift 就是撤消 fftshift 的结果。
如果 Y 是向量，则 ifftshift 会将 Y 的左右两半部分进行交换。
如果 Y 是矩阵，则 ifftshift 会将 Y 的第一象限与第三象限交换，将第二象限与第四象限交换。
如果 Y 是多维数组，则 ifftshift 会沿每个维度交换 Y 的半空间。
X = ifftshift(Y,dim) 沿 Y 的维度 dim 执行运算。例如，如果 Y 是矩阵，其行表示多个一维变换，则 ifftshift(Y,2) 会将 Y 的每一行的左右两半部分进行交换。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def inverse_fft_shift(input_str):
"""
实现类似 MATLAB 的 ifftshift 函数，支持指定维度进行操作。
该函数将数组的零频率分量移到频谱中心。
:param input_str: 输入的字符串，格式为数组或 (数组, 维度)，维度为MATLAB风格的1-based索引
:return: 零频率分量移到中心后的数组（SymPy数组格式）
"""
expr = sp.sympify(input_str)
arr = dim = None
error = False
# 解析输入，判断是否为元组（包含数组和维度参数）
if isinstance(expr, tuple):
if len(expr) < 1:
raise ValueError("输入格式错误，应为 (数组, 维度)")
arr_expr = expr[0]
dim_expr = expr[1] if len(expr) >= 2 else None
arr = arr_expr
# 解析维度参数
if dim_expr is not None:
if isinstance(dim_expr, (sp.Integer, int)):
dim = int(dim_expr)
elif isinstance(dim_expr, (list, tuple, sp.Array)):
dim = [int(d) for d in dim_expr]
else:
raise ValueError("维度参数格式错误")
else:
arr = expr
dim = None
if arr is not None:
# 转换为numpy数组进行处理
if isinstance(arr, list):
x = np.array(arr)
axes = None
# 处理维度参数，将MATLAB的1-based转换为Python的0-based
if dim is None:
axes = tuple(range(x.ndim))
else:
if isinstance(dim, int):
if dim < 1:
raise ValueError("维度必须大于等于1")
axes = (dim - 1,)
elif isinstance(dim, (list, tuple)):
axes = []
for d in dim:
if d < 1:
raise ValueError("维度必须大于等于1")
axes.append(d - 1)
axes = tuple(axes)
else:
raise ValueError("维度参数必须是整数、整数列表或整数元组")
# 检查轴是否有效
for ax in axes:
if ax < 0 or ax >= x.ndim:
raise ValueError(f"轴{ax + 1}超出数组的维度范围（1-based）")
# 计算每个轴的移动量（取半）
shift = [x.shape[ax] // 2 for ax in axes]
# 执行循环位移
result = np.roll(x, shift=shift, axis=axes)
else:
error = True
if not error:
return sp.Array(result.tolist()) # 转换回SymPy数组
else:
return f"输入错误: {input_str}"
if __name__ == "__main__":
# 示例 1: 一维数组，不指定维度
Y1 = [0, 1, 2, 3, 4]
print("示例 1 - 原始一维数组:")
result1 = inverse_fft_shift(f"{Y1}")
print("示例 1 - 不指定维度经过ifftshift处理后的一维数组:")
print(np.array(result1.tolist()))
# [3 4 0 1 2]
# 示例 2: 二维数组，沿第1个维度（行方向）执行ifftshift
Y2 = [[0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11]]
print("\n示例 2 - 原始二维数组:")
result2 = inverse_fft_shift(f"{Y2}, 1")
print("示例 2 - 沿第1个维度经过ifftshift处理后的二维数组:")
print(np.array(result2.tolist()))
# [[8 9 10 11]
# [0 1 2 3]
# [4 5 6 7]]
# 示例 3: 二维数组，沿第2个维度（列方向）执行ifftshift
print("\n示例 3 - 原始二维数组（同示例 2）:")
result3 = inverse_fft_shift(f"{Y2}, 2")
print("示例 3 - 沿第2个维度经过ifftshift处理后的二维数组:")
print(np.array(result3.tolist()))
# [[2 3 0 1]
# [6 7 4 5]
# [10 11 8 9]]
# 示例 4: 二维数组，同时沿第1和第2个维度执行ifftshift
print("\n示例 4 - 原始二维数组（同示例 2）:")
result4 = inverse_fft_shift(f"{Y2}, [1, 2]")
print("示例 4 - 同时沿第1和第2个维度经过ifftshift处理后的二维数组:")
print(np.array(result4.tolist()))
# [[10 11 8 9]
# [2 3 0 1]
# [6 7 4 5]]
# 示例 5: 三维数组，沿第3个维度执行ifftshift
Y3 = [
[[0, 1], [2, 3]],
[[4, 5], [6, 7]]
]
print("\n示例 5 - 原始三维数组:")
result5 = inverse_fft_shift(f"{Y3}, 3")
print("示例 5 - 沿第3个维度经过ifftshift处理后的三维数组:")
print(np.array(result5.tolist()))
# [[[1 0]
# [3 2]]
# [[5 4]
# [7 6]]]
逆傅里叶变换
ifourier(F) 返回 F 的逆傅里叶变换. 默认情况下独立变量为w, 变换变量为x. 如果F不包含w, ifourier将使用函数symvar.
ifourier(F,transVar) 使用变换变量 transVar 而不是 x.
ifourier(F,var,transVar) 分别使用独立变量 var 和变换变量 transVar 而不是 w 和 x.
F — 输入,符号表达式,符号函数,符号向量,符号矩阵
var — 自变量,x（默认）| 符号变量
自变量，指定为符号变量。此变量通常称为“频率变量”。如果您未指定变量，则 ifourier 使用 w。如果 F 不包含 w，则 ifourier 使用函数 symvar 来确定自变量。
transVar — 转换变量, x(默认), t | 符号变量 | 符号表达式 | 符号向量 | 符号矩阵
变换变量，指定为符号变量、表达式、向量或矩阵。它通常被称为“时间变量”或“空间变量”。默认情况下，ifourier 使用 x。如果 x 是 F 的独立变量，则 ifourier 使用 t。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def inverse_fourier_transform(input_str):
"""
执行逆傅里叶变换，对标MATLAB的ifourier函数
参数:
input_str (str): 输入字符串，支持以下格式：
- 单参数: "F"
- 双参数: "(F, transVar)"
- 三参数: "(F, var, transVar)"
返回:
sympy.Expr/sympy.Matrix/str: 计算结果或错误信息
示例:
>>> inverse_fourier_transform("exp(-w**2)") # 默认w→x变换
sqrt(pi)*exp(-x**2/4)/(2*sqrt(pi))
>>> inverse_fourier_transform("(exp(-s**2), t)") # 指定转换变量
sqrt(pi)*exp(-t**2/4)/(2*sqrt(pi))
"""
try:
# 预处理输入字符串
expr = sp.sympify(input_str)
error = False
result = None
def elementwise_ifourier(f_val, var_val=None, transVar_val=None):
# 确定变量
if var_val is None:
# 查找符号变量，优先使用w，否则取第一个
symbols_in_f = f_val.free_symbols
if not symbols_in_f:
raise ValueError("无符号变量")
var_val = sp.symbols('w') if sp.symbols('w') in symbols_in_f else next(iter(symbols_in_f))
if transVar_val is None:
# 转换变量默认为x，如果输入变量是x则用t
transVar_val = sp.symbols('x') if var_val != sp.symbols('x') else sp.symbols('t')
# 按照MATLAB的定义计算逆傅里叶变换: (1/(2π)) ∫ F(ω) e^{iωt} dω
integral_expr = f_val * sp.exp(sp.I * var_val * transVar_val)
result = (1 / (2 * sp.pi)) * sp.integrate(integral_expr, (var_val, -sp.oo, sp.oo))
simplified_result = result.simplify()
# 如果结果是分段函数，遍历分支找第一个非零项
if isinstance(simplified_result, sp.Piecewise):
for arg in simplified_result.args: # arg结构: (表达式, 条件)
expr, _ = arg
simplified_expr = expr.simplify()
# 检查表达式是否明确非零（考虑符号计算的不确定性）
if simplified_expr.is_zero is not True: # 允许None（无法判断时默认保留）
return simplified_expr
# 所有分支都明确为零时返回0
return sp.S.Zero
else:
return simplified_result
if isinstance(expr, tuple):
if len(expr) == 3:
n = expr[0]
a = expr[1]
x = expr[2]
elif len(expr) == 2:
n = expr[0]
a = None
x = expr[1]
else:
error = True
result = elementwise_ifourier(n, a, x)
else:
result = elementwise_ifourier(expr)
return result if not error else f"输入错误: {input_str}"
except sp.SympifyError:
return "错误: 输入表达式解析失败"
except IndexError:
return "错误: 未找到有效符号变量"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1：基本用法 (exp(-w^2) → 高斯函数)
print("示例1:")
print(inverse_fourier_transform("exp(-w**2/4)"))
# exp(-x**2)/sqrt(pi)
# 示例2：指定转换变量 (s→t)
print("\n示例2:")
print(inverse_fourier_transform("(exp(-s**2), t)"))
# exp(-t**2/4)/(2*sqrt(pi))
# 示例3：符号表达式保留
print("\n示例3:")
print(inverse_fourier_transform("F(w)"))
# Integral(F(w)*exp(I*w*x), (w, -oo, oo))/(2*pi)
逆希尔伯特变换
f = ihtrans(H) 返回符号函数H的逆希尔伯特变换. 默认情况下, 独立变量为x, 变换变量为t.
H — 输入, 符号表达式, 符号函数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
def inverse_hilbert_transform(input_str):
"""
计算给定符号表达式的逆希尔伯特变换，与MATLAB的ihtrans对齐
参数:
input_str (str): 输入的时间域符号表达式，变量为t
返回:
tuple: (时间数组, 逆希尔伯特变换结果) 或错误信息字符串
"""
try:
expr = sp.sympify(input_str)
t = sp.symbols('t')
# 检查变量合法性
if not expr.has(t):
return "错误: 表达式必须包含变量t"
if len(expr.free_symbols - {t}) > 0:
print("警告: 检测到额外符号，已自动替换为1")
# 符号替换与数值化
substituted = expr.subs({s: 1 for s in expr.free_symbols - {t}})
y_func = sp.lambdify(t, substituted, 'numpy')
# 生成优化时间轴 (避免端点奇点)
t_vals = np.linspace(-np.pi, 3 * np.pi, 3000)[500:2500] # 取中间稳定段
y = y_func(t_vals)
# 使用镜像延拓处理边界效应
y_padded = np.concatenate((-y[::-1], y, -y[::-1]))
analytic = signal.hilbert(y_padded)
H = -np.imag(analytic)[len(y):2 * len(y)] # 取中间有效段
return (t_vals, H)
except Exception as e:
return f"错误: {str(e)}"
# 测试案例1: H(sin(t)) = -cos(t) => iH(-cos(t))应得sin(t)
t, H = inverse_hilbert_transform("-cos(t)")
plt.figure(figsize=(10, 4))
plt.plot(t, H, label='数值解')
plt.plot(t, np.sin(t), 'r--', lw=1, label='理论解')
plt.title('案例1: iH(-cos(t)) vs sin(t)')
plt.legend()
plt.show()
# 测试案例2: H^{-1}(sin(t)) = cos(t) (验证逆变换特性)
t, H = inverse_hilbert_transform("sinc(t)")
plt.figure(figsize=(10, 4))
plt.plot(t, H, label='数值解')
plt.plot(t, np.cos(t), 'g--', lw=1, label='理论解')
plt.title('案例2: iH(sin(t)) vs cos(t)')
plt.legend()
plt.show()
# 精度验证 (检查最大相对误差)
rel_error = np.max(np.abs(H - np.cos(t))) / np.max(np.abs(np.cos(t)))
print(f"最大相对误差: {rel_error:.2e}")
拉普拉斯函数逆变换
f=ilaplace(F)返回f的拉普拉斯逆变换. 默认情况下,自变量为s,变换变量为t.如果f不包含s,ilaplace使用函数symvar.
f=ilaplace(F,transVar)使用转换变量transVar而不是t.
f=ilaplace(F,var,transVar)分别使用自变量var和变换变量transVar代替s和t.
F —— 输入,符号表达式,符号函数
var —— 自变量,s(默认), 符号变量. 自变量, 指定为符号变量,表达式,向量或矩阵. 此变量通常称为“复频率变量”. 如果您未指定变量, 则 ilaplace 使用s 如果 F 不包含 s,则 ilaplace 使用函数 symvar 来确定自变量.
transVar —— 变换变量,t(默认值),变换变量,指定为符号变量,表达式,向量或矩阵.它通常被称为“时间变量”或“空间变量”.默认情况下, ilaplace使用t. 如果t是F的独立变量, 则ilaplace使用x.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def inverse_laplace_transform(input_str):
"""
对标 MATLAB 的 ilaplace 函数，实现逆拉普拉斯变换
参数:
input_str: 输入参数字符串，支持以下格式：
"F" -> 自动识别变量
"F, s" -> 指定原域变量
"F, s, t" -> 指定原域变量和时域变量
"F, t" -> 直接使用时域变量
返回:
逆变换结果，形式与输入一致（标量、矩阵或符号表达式）。错误时返回错误信息
示例:
>>> inverse_laplace_transform("1/s")
Heaviside(t)
>>> inverse_laplace_transform("1/(s^2 + a^2), s, t")
sin(a*t)*Heaviside(t)/a
"""
try:
# 将输入字符串转换为 sympy 表达式
expr = sp.sympify(input_str)
# 初始化错误标志
error = False
# 初始化结果变量
result = None
def eval_inverse_laplace(f_val, var_val=None, transVar_val=None):
"""
计算逆拉普拉斯变换的辅助函数
参数:
f_val: 待变换的表达式
var_val: 原域变量，默认为 None
transVar_val: 时域变量，默认为 None
返回:
逆拉普拉斯变换的结果
"""
# 定义默认的原域变量 s
s = sp.symbols('s')
# 如果 var_val 未提供
if var_val is None:
# 若表达式中包含 s，则使用 s 作为原域变量，否则使用表达式中的第一个自由符号
var_val = s if f_val.has(s) else tuple(f_val.free_symbols)[0]
# 如果 transVar_val 未提供
if transVar_val is None:
# 如果原域变量不是 t，则使用 t 作为时域变量，否则使用 x
default_symbol = 't' if var_val != sp.symbols('t') else 'x'
# 创建正实数的时域变量
transVar_val = sp.symbols(default_symbol, positive=True)
else:
# 将输入的时域变量转换为正实数的符号
transVar_val = sp.Symbol(str(transVar_val), positive=True)
# 调用 sympy 的逆拉普拉斯变换函数进行计算
return sp.inverse_laplace_transform(f_val, var_val, transVar_val)
# 如果输入表达式是元组
if isinstance(expr, tuple):
if len(expr) == 3:
# 分别提取表达式、原域变量和时域变量
n = expr[0]
a = expr[1]
x = expr[2]
elif len(expr) == 2:
# 提取表达式和时域变量，原域变量设为 None
n = expr[0]
a = None
x = expr[1]
else:
# 输入格式错误，设置错误标志
error = True
# 逆拉普拉斯变换
result = eval_inverse_laplace(n, a, x)
# 如果输入表达式是数字或者包含自由符号
elif expr.is_number or expr.free_symbols:
# 直接进行逆拉普拉斯变换
result = eval_inverse_laplace(expr)
else:
# 输入格式错误，设置错误标志
error = True
# 如果没有错误，返回结果，否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print(inverse_laplace_transform("1/s"))
# 1
print(inverse_laplace_transform("exp(-a*s)/s"))
# InverseLaplaceTransform(exp(-a*s)/s, s, t, _None)
print(inverse_laplace_transform("1/(s-a)**2, x"))
# x*exp(a*x)
热图
ImagescPlot(C, coloroption)将阵列C中的数据显示为使用颜色图中的全部颜色范围的图像.C的每个元素指定图像的一个像素的颜色.所得到的图像是像素的m×n网格，其中m是C中的行数，n是列数.元素的行和列索引确定相应像素的中心.
C - 数值矩阵
coloroption:
内置默认序列：
'viridis'（默认）、'plasma'、'inferno'、'magma'、'cividis'（这些是基于科学可视化的优化配色）
'Greys'、'YlGnBu'、'Greens'、'YlOrRd'、'Bluered'、'RdBu' 等（经典渐变色）
分类配色：
单色系：'Greys'、'Reds'、'Blues'、'Oranges' 等（名称均为复数形式）
双色对比：'RdBu'（红 - 蓝）、'PuOr'（紫 - 橙）、'BrBG'（棕 - 蓝绿）等
温度分布图（气象学）:
红色区域为高温区（沙漠），黄色区域为低温区（沿海）
ImagescPlot([[22,24,26,28], # 沿海区域
[28,30,32,34], # 平原区域
[35,37,38,40]], # 内陆沙漠
YlOrRd)
基因表达热图（生物信息学）:
高表达基因呈黄色（如基因B在样本3），低表达呈紫色
ImagescPlot([[0.1,1.2,0.5], # 基因A
[2.8,0.3,4.1], # 基因B
[1.5,3.0,0.7]], # 基因C
viridis)
用户行为分析（网页点击热图）:
深蓝色块为高点击区域（主内容区左侧）
ImagescPlot([[5,20,3], # 顶部导航栏
[50,2,10], # 主内容区
[1,4,30]], # 侧边广告
Blues)
金融相关性矩阵（量化分析）:
红色表示正相关（A与B），蓝色表示负相关（A与C）
ImagescPlot([[1.0,0.7,-0.2], # 股票A
[0.7,1.0,0.1], # 股票B
[-0.2,0.1,1.0]], # 股票C
RdBu)
二维隐式函数图
ImplicitPlot(f(x,y)) 绘制由隐式方程（即未解出因变量的方程）定义的函数关系的图像, 隐式函数图 必须包含两个变量，因为描述二维平面中的曲线, 展示变量间的相互约束, 实际系统常涉及多变量交互.
天体轨道（开普勒定律）: 航天器轨道计算、深空探测路径规划
ImplicitPlot(x^2/4+y^2=1,x=[-15,15],y=[-15,15])
电磁场分布（麦克斯韦方程）: 两个同号点电荷的等势线分布, 高压输电线周围的电场强度梯度
ImplicitPlot(x^2-y^2=1,x=[-@pi,@pi],y=[-@pi,@pi])
医学成像（CT重建）: 描述人体器官的断层扫描轮廓, 肿瘤边缘检测（灰度突变边界）
ImplicitPlot(x^2+y^2-exp(-(x^2+y^2))=0,x=[-3,3],y=[-3,3])
化学反应平衡: 氨合成反应
ImplicitPlot((x^2)*y=0.5,x=[0,1], y=[0,2])
声波干涉: 噪声消除耳机中的反相声波叠加, 音乐会厅声学设计
ImplicitPlot(cos(x)+cos(y)=0,x=[-2@pi,2@pi], y=[-2@pi,2@pi])
神经网络决策边界: 猫狗分类示例
ImplicitPlot(1/(1+exp(-(0.8x-0.5y+0.1)))= 0.5,x=[-5,5], y=[-5,5])
流体力学（伯努利方程）: 机翼表面压力分布
ImplicitPlot(0.5*1.225*x^2+P=101350,x=[0,300], P=[80000,120000])
生态学种群竞争:
参数赋值（澳洲野猫(y) vs 袋鼠(x)）：
a=0.3（野猫捕食强度系数）
ImplicitPlot(x*(1-x-0.3y)=0,x=[0,2],y=[0,4])
三维隐式函数图
ImplicitPlot(f(x,y,z)) 表示的是一个由方程 F(x, y, z) = 0 定义的曲面（或点集）。
它不像显式函数 z = f(x, y) 那样直接解出 z，而是通过检查空间中每个点 (x, y, z) 是否满足方程来确定该点是否在曲面上。
这种表示方法非常强大，可以描述极其复杂的拓扑结构和光滑曲面，尤其是在几何建模、科学计算可视化（如等值面）、物理（如势能面）和数学研究中应用广泛。
球体 (最基础): 最基本的等距曲面。
表示天体（理想化）、粒子模型、约束条件（如分子动力学中的范德华半径约束的简化）、各向同性场（如点电荷的等势面）。
ImplicitPlot3D(x^2+y^2+z^2-1=0,x=[-1.5,1.5],y=[-1.5,1.5],z=[-1.5,1.5])
圆柱面: 表示无限长圆柱体的一部分。
应用：管道、杆件、轴、激光束的理想化模型、旋转对称结构。
ImplicitPlot3D(x^2+y^2-0.64=0,x=[-2,2],y=[-2,2],z=[-2,2])
环面 (甜甜圈): 拓扑学基本曲面（亏格1）。
应用：轮胎、环形线圈（电感器、托卡马克核聚变装置）、环形分子结构、DNA环、某些建筑结构。
ImplicitPlot3D((x^2+y^2+z^2+0.91)^2-4(x^2+y^2),x=[-1.5,1.5],y=[-1.5,1.5],z=[-0.5,0.5])
圆锥面: 表示锥形结构。
应用：喇叭口、光线锥（几何光学）、某些结晶习性、应力集中区域（断裂力学）、声波或冲击波传播的理想化模型（特征线）。
ImplicitPlot3D(x^1+y^2-z^2,x=[-1,1],y=[-1,1],z=[-1,1])
双曲面 (单叶): 具有负高斯曲率的直纹曲面。
应用：冷却塔（经典的双曲面结构，强度高且自然通风效果好）、核反应堆安全壳、某些现代建筑、齿轮几何、天体轨道（双曲线轨道是平面上的双曲线，但双曲面是其空间推广）。
ImplicitPlot3D(x^1+y^2-z^2-0.25=0,x=[-2,2],y=[-2,2],z=[-2,2])
Goursat 曲面 (复杂拓扑): 一类高亏格的代数曲面（亏格>1）。
应用：复杂的拓扑结构研究、几何建模中的复杂形状原型、数学可视化（展示高亏格曲面）、某些多孔材料或泡沫的理想化模型。
ImplicitPlot3D(x^4+y^4+z^4-(x^2*y^2+y^2*z^2+z^2*x^2)+0.1=0,x=[-1.5,1.5],y=[-1.5,1.5],z=[-1.5,1.5])
球形变形: 高次曲面
ImplicitPlot3D(x^4+y^4+z^4-x^2-y^2-z^2+0.3=0,x=[-2,2],y=[-2,2],z=[-2,2])
脉冲响应图
impulseplot(sys)描述了一个线性时不变系统（LTI）在输入为单位脉冲信号（Dirac Delta 函数）时的输出响应。
它直观展示了系统的动态特性（如稳定性、振荡频率、衰减速度等）。数学上，它是传递函数 H(s)H(s) 的拉普拉斯逆变换。
电路系统（RLC 电路）: 二阶低通滤波器
衰减振荡（因复数极点 -1±2j）。
电容电压在瞬时电流脉冲下的响应，振荡频率 2 rad/s，衰减率由实部 -1 决定。
impulseplot(1,s^2+2s+5)
机械系统（悬架减震）: 车辆悬架模型
初始快速上升（分子零点影响），后衰减振荡。
车轮受瞬时冲击（如过路障）时车身的位移响应，振荡频率反映系统固有特性。
impulseplot(0.5s+1,s^2+1.2s+4)
控制系统（无人机姿态调整）: 俯仰角控制器
快速上升至峰值，后衰减振荡至 0。
无人机受瞬时阵风扰动后的俯仰角恢复过程，超调量和调整时间反映控制器性能。
impulseplot(10*(s+0.5),s^2+3s+10)
地震工程（建筑结构抗震）: 3层建筑简化模型
初期剧烈振荡 → 结构受冲击时的晃动
后期指数衰减 → 阻尼材料消耗能量
模拟地震瞬间（脉冲输入）后各楼层的位移响应，衰减速度反映结构耗能能力。
impulseplot(0.02s+0.1,s^2+0.8s+4)
金融系统（市场冲击分析）: 股票指数对突发事件的响应
初始骤降 → 恐慌性抛售（如黑天鹅事件）
振荡恢复 → 市场自我修正过程
稳态归零 → 冲击被完全消化
2020年疫情爆发时美股熔断的脉冲响应拟合。
impulseplot(0.3s-2,s^2+1.5s+0.8)
机器人学（机械臂关节控制）: 谐波减速器关节
初始超调（20%）→ 关节柔性导致的过冲
高频衰减振荡（8Hz）→ 结构共振频率
2秒内收敛 → 控制系统带宽足够
通过调整响应曲线降低超调，避免精密操作时抖动。
impulseplot(8s+16,s^2+4s+64)
神经科学（神经元电信号传导）
1ms内产生尖峰 → 动作电位激发
负向波动 → 钾离子外流导致的超极化
长尾衰减（10ms）→ 不应期恢复过程
解释癫痫患者为何对连续刺激更敏感（脉冲响应叠加效应）。
impulseplot(100,s^2+20s+100)
低通滤波器（平滑噪声）
响应呈5个等高的矩形脉冲（长度=5），表明系统对输入进行短期平均，快速衰减高频噪声。
温度传感器数据去噪，消除瞬时干扰。
impulseplot([1],[0.2,0.2,0.2,0.2,0.2])
弹簧-阻尼系统（机械振动）
响应呈衰减振荡（如 [1.0, 1.2, 0.64, -0.128, ...]），振幅逐渐减小。
分析悬架对路面冲击的减震效果，振荡频率和衰减速率决定舒适性。
impulseplot([1],[1,-1.2,0.8])
经济预测模型（时间序列分析）
指数衰减序列（如 [0.6, 0.48, 0.384, ...]），反映冲击的持久影响。
评估经济政策冲击的持续时间（如加息对通胀的影响）。
impulseplot([0.6],[1,-0.8])
数字通信系统（码间干扰分析）
主峰两侧有对称拖尾（如 [0.4, 1.0, 0.4]），主峰后出现非零值。
设计均衡器消除码间干扰（ISI），提升数据传输可靠性。
impulseplot([1],[0.4,1.0,0.4])
生态模型（种群动态）
振荡持续但缓慢衰减（如初始响应 [1, 0.8, 0.46, ...]）。
预测物种数量对突发事件的恢复能力（如疾病或食物短缺）。
impulseplot([1,0.3],[1,-0.5,-0.2])
叠加脉冲响应图
impulseplot([sys1],[sys2]...[sysN])
叠加脉冲响应图是在同一坐标系中绘制多个系统的脉冲响应曲线，用于直观比较不同系统的动态特性（如响应速度、稳定性、振荡行为等）。
悬架系统设计对比
H1：单调上升后平缓衰减 → 基础减震效果
H2：初始快速响应（0.2s达峰）→ 主动控制优势
H2使豪华SUV过减速带的冲击加速度降低45%
impulseplot([2.5s+8,s^2+1.8s+16], # 系统1: 传统弹簧减震器
[3.2s^2+12s+40,s^3+4.5s^2+22s+60]) # 系统2: 磁流变智能悬架
医疗输液泵控制算法
H1超调量≈0%（阶跃响应验证）→ 适合化疗药物输送
H2上升时间快40% → 急救时快速补液
H2在休克抢救中使血压达标时间缩短至1.8分钟
impulseplot([12,s^2+4s+12], # 系统1: PID控制器
[8s+30,s^2+2.5s+25]) # 系统2: 模糊自适应控制器
无人机抗风性能
H2在6级风中姿态稳定性提升至H1的2.3倍
impulseplot([15*(s+2),s^2+5s+50], # 系统1: 标准PID控制器
[25*(s+3),s^2+8s+100]) # 系统2: 神经网络补偿控制器
声学降噪耳机
H1：单纯指数衰减 → 物理隔音
H2：初始负向响应 → 主动声波抵消
impulseplot([10,s+15], # 系统1: 被动降噪
[-6s-50,s^2+20s+400]) # 系统2: 主动降噪
定积分和不定积分
F=int(expr)计算expr的不定积分. int使用由symvar(expr,1)确定的默认积分变量.如果expr是一个常数,那么默认的积分变量是x.
F=int(expr,var)计算expr相对于符号标量变量var的不定积分.
F=int(expr,a,b）计算expr从a到b的定积分.int使用由symvar(expr,1)确定的默认积分变量.如果expr是一个常数,那么默认的积分变量是x.
F=int(expr,var,a,b)计算expr相对于符号标量变量var从a到b的定积分.
expr — 整数,符号表达式,符号函数,符号向量,符号矩阵,符号数
var — 积分变量,符号变量
a —— 下限,数字,符号数,符号变量,符号表达式,符号函数
b —— 上限,数字,符号数,符号变量,符号表达式,符号函数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def integrate_equation(input_str):
"""
对标MATLAB int函数，支持定积分、不定积分和矩阵积分
参数:
input_str: 字符串形式的数学表达式，支持以下格式:
- 不定积分: "x**2", "sin(y)"
- 定积分: "(x**2, x, 0, 1)", "(f, (x, a, b))"
- 矩阵积分: "[[x, x**2], [sin(x), cos(x)]]"
返回:
SymPy表达式或矩阵，若出错返回错误信息字符串
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def element_integrate(ele):
if ele.free_symbols:
sym_var = tuple(ele.free_symbols)[0] # 注意：可能因符号顺序导致问题
return sp.integrate(ele, sym_var).doit()
elif ele.is_number:
sym_var = sp.symbols('x')
return sp.integrate(ele, sym_var).doit()
else:
return sp.nan
# 处理四元组定积分形式 (被积函数, 积分变量, 下限, 上限)
if isinstance(expr, tuple) and len(expr) == 4:
if expr[0].free_symbols and expr[1].free_symbols: # 标量积分
result = sp.integrate(expr[0], (expr[1], expr[2], expr[3])).doit()
else:
error = True
# 处理三元组形式 (被积函数, 下限, 上限) 自动提取积分变量
elif isinstance(expr, tuple) and len(expr) == 3:
sym_var = tuple(expr[0].free_symbols)[0] # 提取第一个符号作为积分变量
if expr[0].free_symbols: # 标量积分
result = sp.integrate(expr[0], (sym_var, expr[1], expr[2])).doit()
else:
error = True
# 处理二元组形式 (被积函数, 积分变量) 不定积分
elif isinstance(expr, tuple) and len(expr) == 2:
if all(e.free_symbols for e in expr):
result = sp.integrate(*expr).doit()
else:
error = True
# 处理单个表达式不定积分（自动提取变量）
elif expr.free_symbols or expr.is_number:
result = element_integrate(expr)
else:
error = True
return sp.simplify(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示范代码
if __name__ == "__main__":
# 示例1: 不定积分
print("--- 不定积分示例 ---")
print("∫x² dx =", integrate_equation("x**2"))
# x**3/3
# 示例2: 定积分 (显式指定变量)
print("\n--- 定积分示例（四元组形式）---")
print("∫₀¹ x² dx =", integrate_equation("(x**2, x, 0, 1)"))
# 1/3
# 示例3: 定积分 (自动提取变量)
print("\n--- 定积分示例（三元组形式）---")
print("∫₀^π sin(x) dx =", integrate_equation("(sin(x), 0, pi)"))
# 2
数值积分
q = integral(fun,xmin,xmax) 使用全局自适应积分和默认误差容限在 xmin 至 xmax 间以数值形式为函数 fun 求积分.
fun — 被积函数,函数句柄
xmin — x 的下限,实数,复数
xmax — x 的上限,实数,复数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def integrate_numerical_equation(input_str):
"""
对标MATLAB的integral函数，进行数值积分计算
参数:
input_str: 输入的字符串，表示积分表达式和积分区间
返回:
积分结果 或 错误信息
"""
try:
# 将输入的字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 处理输入为元组且长度为3的情况，例如 (x**2, 0, 1)
if isinstance(expr, tuple) and len(expr) == 3:
# 找到被积函数中的自由符号作为积分变量
x = list(expr[0].free_symbols)[0]
# 构建积分表达式并计算数值积分
result = sp.Integral(expr[0], (x, expr[1], expr[2])).evalf()
# 处理输入为元组且长度大于3的情况，例如 (x**2, x, 0, 1)
elif isinstance(expr, tuple) and len(expr) > 3:
# 构建积分表达式并计算数值积分
result = sp.Integral(expr[0], (expr[1], expr[2], expr[3])).evalf()
else:
# 若输入格式不符合要求，标记为错误
error = True
# 如果没有错误，返回积分结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 若计算过程中出现错误，返回错误信息
return f"错误:{e}"
# 示范代码
# 示例1：计算 x**2 在 [0, 1] 上的积分
input_str1 = "(x**2, 0, 1)"
result1 = integrate_numerical_equation(input_str1)
print(f"示例1结果: {result1}")
# 0.333333333333333
对二重积分进行数值计算
q = integral2(fun,xmin,xmax,ymin,ymax) 在平面区域 xmin ≤ x ≤ xmax 和 ymin(x) ≤ y ≤ ymax(x) 上逼近函数 z = fun(x,y) 的积分
q = integral2(fun,xmin,xmax,ymin,ymax,Name,Value) 使用一个或多个 Name,Value 对组参量指定其他选项
fun — 被积函数,函数句柄
xmin — x 的下限,实数,复数
xmax — x 的上限,实数,复数
ymin — y 的下限,实数,复数
ymax — y 的上限,实数,复数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def integrate_multi_equation(input_str):
"""
对标MATLAB integral2函数的多重数值积分实现
参数:
input_str: 输入字符串，格式为：
"函数表达式, 变量1下限, 变量1上限, 变量2下限, 变量2上限,..."
示例："x*y, 0, 1, y, 0, 1-x"
**kwargs: 传递给scipy积分函数的参数 (如epsabs, epsrel等)
返回:
成功时返回积分结果和误差组成的元组，失败时返回错误信息
示例:
>>> integrate_multi_equation("x*y, 0, 1, 0, 1-x")
(0.04166666666666668, 4.623514334137052e-16)
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
# 初始化错误标志为False
error = False
# 初始化积分结果为None
result = None
def multi_integral(f, *args):
# 提取传递的所有积分范围，由于每个积分范围由下限和上限两个值组成，所以除以2
num_ranges = len(args) // 2
# 找到函数f中的符号变量，并按符号名称排序
symbols_in_f = sorted(f.free_symbols, key=lambda s: str(s))
# 检查符号数量和积分范围是否匹配，如果积分范围数量多于符号变量数量，抛出异常
if num_ranges > len(symbols_in_f):
raise ValueError("传入的积分范围比符号变量多。")
# 从外到内嵌套integrate，处理多重积分
for i in range(num_ranges):
# 依次取出符号变量
symbol = symbols_in_f[i]
# 取出当前符号变量的积分下限
lower_bound = args[2 * i]
# 取出当前符号变量的积分上限
upper_bound = args[2 * i + 1]
# 对函数f关于当前符号变量进行积分，积分区间为[下限, 上限]
f = sp.integrate(f, (symbol, lower_bound, upper_bound))
return f
# 检查expr是否为元组
if isinstance(expr, tuple):
# 检查元组的第一个元素是否包含符号变量
if expr[0].free_symbols:
# 取出元组的第一个元素作为被积函数
F = expr[0]
# 取出元组中除第一个元素外的其他元素作为积分范围
args = expr[1:]
# 调用multi_integral函数进行多重积分
result = multi_integral(F, *args)
else:
# 如果元组的第一个元素不包含符号变量，设置错误标志为True
error = True
else:
# 如果expr不是元组，设置错误标志为True
error = True
# 如果没有错误，返回积分结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 捕获并处理可能出现的异常，返回错误信息
return f"Error: {e}"
# 二重积分 ∫₀¹∫₀¹ xy dxdy
# 示例1: 简单二重积分 ∫₀¹∫₀¹ x*y dxdy
print("示例1:", integrate_multi_equation("x*y, 0, 1, 0, 1"))
# 1/4
# 示例2: 带变量依赖的积分 ∫₀¹∫₀^{1-x} x*y dydx
print("\n示例2:", integrate_multi_equation("x*y, 0, 1, 0, 1-x"))
# (1 - x)**2/4
# 示例3: 三重积分 ∫₀¹∫₀¹∫₀¹ x*y*z dxdydz
print("\n示例3:", integrate_multi_equation("x*y*z, 0, 1, 0, 1, 0, 1"))
# 1/8
对三重积分进行数值计算
q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax) 在平面区域 xmin ≤ x ≤ xmax 和 ymin(x) ≤ y ≤ ymax(x) 和 zmin(x,y) ≤ z ≤ zmax(x,y) 上逼近函数 z = fun(x,y,z) 的积分
q = integral3(fun,xmin,xmax,ymin,ymax,Name,Value) 使用一个或多个 Name,Value 对组参量指定其他选项
fun — 被积函数,函数句柄
xmin — x 的下限,实数,复数
xmax — x 的上限,实数,复数
ymin — y 的下限,实数,复数
ymax — y 的上限,实数,复数
zmin — z 的下限,实数,复数
zmax — z 的上限,实数,复数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def integrate_multi_equation(input_str):
"""
对标MATLAB integral2函数的多重数值积分实现
参数:
input_str: 输入字符串，格式为：
"函数表达式, 变量1下限, 变量1上限, 变量2下限, 变量2上限,..."
示例："x*y, 0, 1, y, 0, 1-x"
**kwargs: 传递给scipy积分函数的参数 (如epsabs, epsrel等)
返回:
成功时返回积分结果和误差组成的元组，失败时返回错误信息
示例:
>>> integrate_multi_equation("x*y, 0, 1, 0, 1-x")
(0.04166666666666668, 4.623514334137052e-16)
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
# 初始化错误标志为False
error = False
# 初始化积分结果为None
result = None
def multi_integral(f, *args):
# 提取传递的所有积分范围，由于每个积分范围由下限和上限两个值组成，所以除以2
num_ranges = len(args) // 2
# 找到函数f中的符号变量，并按符号名称排序
symbols_in_f = sorted(f.free_symbols, key=lambda s: str(s))
# 检查符号数量和积分范围是否匹配，如果积分范围数量多于符号变量数量，抛出异常
if num_ranges > len(symbols_in_f):
raise ValueError("传入的积分范围比符号变量多。")
# 从外到内嵌套integrate，处理多重积分
for i in range(num_ranges):
# 依次取出符号变量
symbol = symbols_in_f[i]
# 取出当前符号变量的积分下限
lower_bound = args[2 * i]
# 取出当前符号变量的积分上限
upper_bound = args[2 * i + 1]
# 对函数f关于当前符号变量进行积分，积分区间为[下限, 上限]
f = sp.integrate(f, (symbol, lower_bound, upper_bound))
return f
# 检查expr是否为元组
if isinstance(expr, tuple):
# 检查元组的第一个元素是否包含符号变量
if expr[0].free_symbols:
# 取出元组的第一个元素作为被积函数
F = expr[0]
# 取出元组中除第一个元素外的其他元素作为积分范围
args = expr[1:]
# 调用multi_integral函数进行多重积分
result = multi_integral(F, *args)
else:
# 如果元组的第一个元素不包含符号变量，设置错误标志为True
error = True
else:
# 如果expr不是元组，设置错误标志为True
error = True
# 如果没有错误，返回积分结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 捕获并处理可能出现的异常，返回错误信息
return f"Error: {e}"
# 二重积分 ∫₀¹∫₀¹ xy dxdy
# 示例1: 简单二重积分 ∫₀¹∫₀¹ x*y dxdy
print("示例1:", integrate_multi_equation("x*y, 0, 1, 0, 1"))
# 1/4
# 示例2: 带变量依赖的积分 ∫₀¹∫₀^{1-x} x*y dydx
print("\n示例2:", integrate_multi_equation("x*y, 0, 1, 0, 1-x"))
# (1 - x)**2/4
# 示例3: 三重积分 ∫₀¹∫₀¹∫₀¹ x*y*z dxdydz
print("\n示例3:", integrate_multi_equation("x*y*z, 0, 1, 0, 1, 0, 1"))
# 1/8
分部积分
G=integrateByParts(F，du)将分部积分应用于F中的积分,其中微分du被积分.
在F中指定积分时,您可以通过使用int函数并将“Hold”选项设置为true来返回未求值的积分形式.然后，您可以使用integrateByParts按部件显示集成的步骤.
F —— 包含积分的表达式
du —— 待整合的差异
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def integrate_by_parts(input_str):
"""
该函数对标 MATLAB 的 integrateByParts 函数，用于执行分部积分操作。
分部积分公式为 ∫u dv = u*v - ∫v du。
参数:
input_str (str): 输入的字符串，格式为 "被积函数,du表达式"，例如 "u*diff(v, x),diff(u)"。
返回:
sympy 表达式或错误信息字符串: 若计算成功，返回分部积分的结果；若出现错误，返回错误信息。
"""
# 定义符号变量和函数
try:
# 定义符号变量 x
x = sp.symbols('x')
# 定义函数 u(x) 和 v(x)
u = sp.Function('u')(x)
v = sp.Function('v')(x)
# 使用逗号分割输入字符串，获取被积函数和 du 表达式
params = input_str.split(',')
# 检查输入参数数量是否为 2
if len(params) != 2:
raise ValueError("输入格式错误，需要两个参数")
F_str, du_str = params
# 转换 Matlab 语法到 SymPy 表达式
try:
# 处理被积函数（替换 int 为 Integral，并解析）
F_expr = sp.sympify(
F_str,
locals={
'int': sp.Integral,
'diff': sp.diff,
'exp': sp.exp,
'sin': sp.sin,
'u': u,
'v': v,
'x': x,
},
)
# 提取被积函数，例如 u*Derivative(v, x)
integrand = F_expr.args[0]
# 解析微分部分（如 diff(u)）
du = sp.sympify(
du_str,
locals={
'diff': sp.diff,
'exp': sp.exp,
'u': u,
'x': x,
'Symbol': sp.Symbol
},
)
except Exception as e:
# 尝试直接解析输入字符串为元组
expr = sp.sympify(input_str)
if isinstance(expr, tuple) and len(expr) == 2:
F, du = expr
# 从 du 中推导出 u
u = sp.integrate(du, x) # u 是 du 的不定积分
# F = u * dv, 因此 dv = F / u
dv = F / u
# 计算 v
v = sp.integrate(dv, x)
# 返回积分分部公式的结果 u*v - ∫v*du
return u * v - sp.integrate(v * du, x)
# 计算 u 和 dv/dx
try:
# 从 du = du/dx 积分得到 u（假设积分常数为 0）
u_expr = sp.integrate(du, x)
# 计算 dv/dx = 被积函数 / u
dv_dx = integrand / u_expr
# 积分得到 v
v_expr = sp.integrate(dv_dx, x)
except Exception as e:
raise ValueError(f"积分计算错误: {e}")
# 应用分部积分公式
result = u_expr * v_expr - sp.Integral(v_expr * du, x)
return result
except Exception as e:
return f"错误: {str(e)}"
# 示范代码
if __name__ == "__main__":
# 示例输入
input_str = "int(u*diff(v)),diff(u)"
result = integrate_by_parts(input_str)
print("分部积分结果:", result)
# u(x)*v(x) - Integral(v(x)*Derivative(u(x), x), x)
一维数据插值
vq = interp1(x,v,xq) 使用线性插值返回一维函数在特定查询点的插入值。向量 x 包含样本点，v 包含对应值 v(x)。向量 xq 包含查询点的坐标。
如果您有多个在同一点坐标采样的数据集，则可以将 v 以数组的形式进行传递。数组 v 的每一列都包含一组不同的一维样本值。
vq = interp1(x,v,xq,method) 指定备选插值方法：'linear'、'nearest'、'next'、'previous'、'pchip'、'cubic'、'v5cubic'、'makima' 或 'spline'。默认方法为 'linear'。
vq = interp1(x,v,xq,method,extrapolation) 用于指定外插策略，来计算落在 x 域范围外的点。如果希望使用 method 算法进行外插，可将 extrapolation 设置为 'extrap'。您也可以指定一个标量值，这种情况下，interp1 将为所有落在 x 域范围外的点返回该标量值。
x - 样本点, 向量
v - 样本值, 向量 | 矩阵
xq - 查询点, 标量 | 向量 | 矩阵
method - 插值方法, 'linear' （默认） | 'nearest' | 'next' | 'previous' | 'pchip' | 'cubic' | 'v5cubic' | 'makima' | 'spline'
extrapolation - 外插策略, 'extrap' | 标量值
vq - 插入的值, 标量 | 向量 | 矩阵
pp - 分段多项式, 结构体
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.interpolate import interp1d, CubicSpline
def interpolation_one_dimensional(input_str):
"""
对标MATLAB的interp1函数，实现一维数据插值
参数:
input_str: 输入参数字符串，支持以下格式：
"y, xq" -> 自动生成x坐标
"x, y, xq" -> 基本插值
"x, y, xq, method" -> 指定插值方法
"x, y, xq, method, extrap" -> 指定方法和外推值
返回:
插值结果数组，错误时返回错误信息字符串
示例:
>>> interpolation_one_dimensional("[1,2,3], [4,5,6], 2.5")
array(5.5)
>>> interpolation_one_dimensional("[1,2,3], [4,5,6], [0,4], 'linear', 0")
array([4., 0.])
"""
try:
# 解析输入字符串为参数元组
expr = sp.sympify(input_str, evaluate=False)
params = list(expr) if isinstance(expr, tuple) else [expr]
# 有效插值方法集合
valid_methods = {
'linear', 'nearest', 'zero',
'slinear', 'quadratic', 'cubic',
'previous', 'next'
}
# 初始化默认参数
x, y, xq = None, None, None
method = 'linear'
extrap = np.nan
# 参数转换函数
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")
# 处理前三个数值参数
numeric_params = []
for i, param in enumerate(params[:3]):
try:
numeric_params.append(convert(param))
except:
if i < 2:
raise
numeric_params.append(np.array([float(param.evalf())]))
# 参数逻辑判断
if len(numeric_params) == 2: # y, xq 模式
y, xq = numeric_params
x = np.arange(1, len(y) + 1)
else: # x, y, xq 模式
x, y, xq = numeric_params[:3]
if len(x) != len(y):
raise ValueError("x和y的长度必须一致")
# 验证x的单调性
if not np.all(np.diff(x) > 0):
raise ValueError("x必须是严格单调递增的")
# 处理方法参数
if len(params) >= 4:
method_param = params[3]
method_candidate = str(method_param).lower()
if method_candidate in valid_methods:
method = method_candidate
else:
try: # 尝试解析为外推值
extrap = float(sp.sympify(method_param))
except:
pass
# 处理外推参数
if len(params) >= 5:
try:
extrap = float(params[4])
except:
raise ValueError("无效的外推值")
# 创建插值器
# 关键修改点：单独处理三次样条
if method in ['cubic', 'spline']:
# 使用CubicSpline并设置not-a-knot边界条件
if not np.all(np.diff(x) > 0):
raise ValueError("三次样条需要严格单调的x值")
# MATLAB的cubic对应SciPy的CubicSpline with not-a-knot
interpolator = CubicSpline(
x, y,
bc_type='not-a-knot',
extrapolate=extrap is not None
)
else:
# 其他方法保持原有interp1d调用
interpolator = interp1d(
x, y,
kind=method,
fill_value=extrap,
bounds_error=False,
assume_sorted=True
)
# 执行插值计算
return interpolator(xq)
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 基本线性插值
print(interpolation_one_dimensional("[1,2,3], [4,5,6], 2.5"))
# [5.5]
# 自动生成x坐标
print(interpolation_one_dimensional("[4,5,6], [1.5,2.5]"))
# [4.5, 5.5]
# 三次样条插值
print(interpolation_one_dimensional("[0,1,2], [0,1,4], 1.5, 'cubic'"))
# [2.25]
# 外推测试
print(interpolation_one_dimensional("[1,2,3], [4,5,6], 4, 'linear', 0"))
# [0]
meshgrid格式的二维网格数据的插值
Vq = interp2(X,Y,V,Xq,Yq) 使用线性插值返回双变量函数在特定查询点的插入值。结果始终穿过函数的原始采样。X 和 Y 包含样本点的坐标。V 包含各样本点处的对应函数值。Xq 和 Yq 包含查询点的坐标。
Vq = interp2(___,method) 指定备选插值方法：'linear'、'nearest'、'cubic'、'makima'。默认方法为 'linear'。
Vq = interp2(___,method,extrapval) 还指定标量值 extrapval，此参数会为处于样本点域范围外的所有查询点赋予该标量值。
X,Y — 样本网格点, 矩阵 | 向量
V — 样本值,矩阵
Xq,Yq — 查询点, 标量 | 向量 | 矩阵
method — 插值方法, 'linear' (默认) | 'nearest' | 'cubic' | 'makima'
extrapval — X 和 Y 域范围外的函数值, 标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import re
import sympy as sp
from scipy.interpolate import RegularGridInterpolator
def interpolation_two_dimensional(input_str):
"""
对标 MATLAB 的 interp2 函数，实现二维网格数据插值
参数:
input_str: 输入参数字符串，支持以下格式：
"X, Y, V, Xq, Yq" -> 基本插值（默认线性）
"X, Y, V, Xq, Yq, method" -> 指定插值方法
"X, Y, V, Xq, Yq, method, extrap" -> 指定方法和外推值
返回:
插值结果数组，错误时返回错误信息字符串
示例:
>>> interpolation_two_dimensional("[0,1], [0,1], [[1,2],[3,4]], 0.5, 0.5")
array(2.5)
>>> interpolation_two_dimensional("[0,1], [0,1], [[1,2],[3,4]], [0.5,2], [0.5,1], 'linear', 0)")
array([2.5, 0. ])
"""
try:
# ------------------------- 输入解析 -------------------------
# 解析输入字符串为参数元组
expr = sp.sympify(input_str, evaluate=False)
print(expr)
params = list(expr) if isinstance(expr, tuple) else [expr]
# 有效插值方法集合
valid_methods = {
'linear', 'nearest', 'cubic',
'makima' # 注意：SciPy 的 cubic 对应 MATLAB 的 spline
}
# 初始化默认参数
X, Y, V, Xq, Yq = None, None, None, None, None
extrap = np.nan
# 提取前5个必要参数
arrays = []
for i, param in enumerate(params[:5]):
try:
if isinstance(param, list):
arr = np.array(param, dtype=float)
elif param.is_number:
arr = np.array([float(param)])
else:
try:
arr = np.array([float(param.evalf())])
except:
raise ValueError(f"无法转换符号表达式: {param}")
arrays.append(arr.squeeze()) # 去除单一维度
except Exception as e:
if i < 3:
raise ValueError(f"参数 {i + 1} 转换失败: {str(e)}")
else:
# 查询点可能是单个数值
arrays.append(np.array([float(param.evalf())]))
# 分配参数
X, Y, V, Xq, Yq = arrays[:5]
# ------------------------- 参数验证 -------------------------
# 验证网格结构
if X.ndim != 1 or Y.ndim != 1:
raise ValueError("X 和 Y 必须是一维数组")
# 检查单调性
if not (np.all(np.diff(X) > 0) and np.all(np.diff(Y) > 0)):
raise ValueError("X 和 Y 必须是严格单调递增的")
# 验证 V 的维度
expected_shape = (len(Y), len(X)) # MATLAB 的 meshgrid 是 (Y,X) 顺序
if V.shape != expected_shape:
raise ValueError(f"V 的维度应为 {expected_shape}，实际为 {V.shape}")
# ------------------------- 处理方法参数 -------------------------
scipy_method = 'linear'
if len(params) >= 6:
scipy_method = str(params[5])
# ------------------------- 处理外推参数 -------------------------
if len(params) >= 7:
try:
extrap = float(params[6])
except:
raise ValueError("无效的外推值")
# ------------------------- 创建插值器 -------------------------
# 使用 RegularGridInterpolator (适用于规则网格)
interpolator = RegularGridInterpolator(
(Y, X), # 注意顺序：MATLAB 是 (X,Y)，但 SciPy 是 (Y,X)
V,
method=scipy_method,
bounds_error=False,
fill_value=extrap
)
# ------------------------- 准备查询点 -------------------------
# 生成网格点坐标 (兼容标量输入)
Xq_mesh, Yq_mesh = np.meshgrid(Xq, Yq)
query_points = np.column_stack((Yq_mesh.ravel(), Xq_mesh.ravel()))
# ------------------------- 执行插值 -------------------------
Vq = interpolator(query_points)
Vq = Vq.reshape(Xq_mesh.shape)
return Vq
except Exception as e:
return f"错误: {str(e)}"
# 示例使用 ---------------------------------------------------
if __name__ == "__main__":
# 测试案例1：外推值设置
input_str = "[0,1],[0,1],[[1,2],[3,4]],[0.5,2],[0.5,1],linear,0"
result = interpolation_two_dimensional(re.sub(r"\s+", "", input_str))
print("测试1结果:", result)
# [[2.5 0. ]
# [3.5 0. ]]
meshgrid格式的三维网格数据的插值
Vq = interp3(X,Y,Z,V,Xq,Yq,Zq) 使用线性插值返回三变量函数在特定查询点的插值。结果始终穿过函数的原始采样。X、Y 和 Z 包含样本点的坐标。V 包含各样本点处的对应函数值。Xq、Yq 和 Zq 包含查询点的坐标。
Vq = interp3(___,method) 指定备选插值方法：'linear'、'nearest'、'cubic'、'makima' 或 'spline'。默认方法为 'linear'。
Vq = interp3(___,method,extrapval) 还指定标量值 extrapval，此参数会为处于样本点域范围外的所有查询点赋予该标量值。
X,Y,Z — 样本网格点, 数组 | 向量
V — 样本值, 数组
Xq,Yq,Zq — 查询点, 标量 | 向量
k — 细化因子, 1 (默认) | 非负实整数标量
method — 插值方法, 'linear' (默认) | 'nearest'
extrapval — X、Y 和 Z 域范围外的函数值, 标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.interpolate import RegularGridInterpolator
def parse_3d_array(param):
"""递归将 SymPy 嵌套结构转换为三维 numpy 数组"""
if isinstance(param, (sp.Tuple, list, tuple)):
return np.array([parse_3d_array(e) for e in param], dtype=float)
elif param.is_Number:
return float(param.evalf())
else:
raise ValueError("无法解析的三维数组元素")
def interpolation_three_dimensional(input_str):
"""
对标 MATLAB 的 interp3 函数，实现三维网格数据插值
参数:
input_str: 输入参数字符串，支持以下格式：
"X, Y, Z, V, Xq, Yq, Zq" -> 基本插值（默认线性）
"X, Y, Z, V, Xq, Yq, Zq, method" -> 指定插值方法
"X, Y, Z, V, Xq, Yq, Zq, method, extrap" -> 指定方法和外推值
返回:
插值结果数组，错误时返回错误信息字符串
示例:
>>> interpolation_three_dimensional("[0,1], [0,1], [0,1], [[[1,2],[3,4]],[[5,6],[7,8]]], 0.5, 0.5, 0.5")
array(4.5)
"""
try:
import re
# ------------------------- 输入解析 -------------------------
# 移除注释内容（#之后的所有字符）
cleaned_str = re.sub(r'#.*', '', input_str)
expr = sp.sympify(cleaned_str, evaluate=False)
params = list(expr) if isinstance(expr, tuple) else [expr]
valid_methods = {'linear', 'nearest'}
# 初始化参数
X, Y, Z, V, Xq, Yq, Zq = None, None, None, None, None, None, None
extrap = np.nan
method = 'linear'
# 解析前7个参数
arrays = []
for i, param in enumerate(params[:7]):
if i == 3: # 处理三维数组 V
try:
if isinstance(param, (sp.Tuple, list, tuple)):
arr = parse_3d_array(param)
else:
raise ValueError("V 必须是三维数组")
arrays.append(arr)
except Exception as e:
raise ValueError(f"参数4 (V) 转换失败: {str(e)}")
else: # 处理其他参数
try:
if isinstance(param, list):
arr = np.array(param, dtype=float).squeeze()
elif param.is_Number:
arr = np.array([float(param.evalf())])
else:
arr = np.array([float(param.evalf())])
arrays.append(arr.squeeze())
except Exception as e:
raise ValueError(f"参数 {i + 1} 转换失败: {str(e)}")
X, Y, Z, V, Xq, Yq, Zq = arrays[:7]
# ------------------------- 参数验证 -------------------------
# 验证一维数组
for arr, name in zip([X, Y, Z], ['X', 'Y', 'Z']):
if arr.ndim != 1:
raise ValueError(f"{name} 必须是一维数组")
# 验证严格单调递增
for arr, name in zip([X, Y, Z], ['X', 'Y', 'Z']):
if not np.all(np.diff(arr) > 0):
raise ValueError(f"{name} 必须是严格单调递增的")
# 验证 V 的维度
expected_shape = (len(Z), len(Y), len(X))
if V.shape != expected_shape:
raise ValueError(f"V 的维度应为 {expected_shape}，实际为 {V.shape}")
# ------------------------- 处理可选参数 -------------------------
if len(params) >= 8:
method = str(params[7])
if method not in valid_methods:
raise ValueError(f"无效的插值方法: {method}，支持的方法: {valid_methods}")
if len(params) >= 9:
try:
extrap = float(params[8])
except:
raise ValueError("无效的外推值")
# ------------------------- 创建插值器 -------------------------
interpolator = RegularGridInterpolator(
(Z, Y, X), # 注意坐标顺序
V,
method=method,
bounds_error=False,
fill_value=extrap
)
# ------------------------- 生成查询点 -------------------------
# 创建三维网格坐标
Zq_mesh, Yq_mesh, Xq_mesh = np.meshgrid(Zq, Yq, Xq, indexing='ij')
query_points = np.column_stack((
Zq_mesh.ravel(),
Yq_mesh.ravel(),
Xq_mesh.ravel()
))
# ------------------------- 执行插值 -------------------------
Vq = interpolator(query_points)
Vq = Vq.reshape(Zq_mesh.shape)
return Vq
except Exception as e:
return f"错误: {str(e)}"
# 测试三维线性插值
input_str = "[0,1], [0,1], [0,1], [[[1,2],[3,4]],[[5,6],[7,8]]], 0.5, 0.5, 0.5"
result = interpolation_three_dimensional(input_str)
print(result)
# [[[4.5]]]
一维插值(FFT 方法)
y = interpft(X,n) 在 X 中内插函数值的傅里叶变换以生成 n 个等间距的点。interpft 对第一个大小不等于 1 的维度进行运算。
y = interpft(X,n,dim) 沿维度 dim 运算。例如，如果 X 是矩阵，interpft(X,n,2) 将在 X 行上进行运算。
X — 输入数组, 向量 | 矩阵
n — 点的数目, 正整数标量
dim — 沿其运算的维度, 正整数标量
y — 插入的点, 向量 | 矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.fft import fft, ifft
import sympy as sp
def interpolation_fft(input_str):
"""
实现类似于 MATLAB 的 interpft 函数，通过 FFT 方法进行一维插值。
参数:
x (array-like): 输入信号（实数或复数）。
n (int): 插值后的信号长度，必须大于原信号长度。
返回:
ndarray: 插值后的实数信号（虚部被丢弃）。
异常:
ValueError: 如果 n 小于原信号长度。
示例:
>>> import numpy as np
>>> x = np.array([1, -1, 1, -1])
>>> interpft(x, 6)
array([ 1. , -0.5, -0.5, 1. , -0.5, -0.5])
"""
try:
# 将输入的字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 标记是否出现错误
error = False
# 存储插值结果
result = None
def interpft_single(x, n):
"""
对一维数组进行插值的辅助函数
"""
# 获取输入信号的长度
m = len(x)
# 检查插值后的长度是否小于原信号长度
if n < m:
raise ValueError("插值后的长度 n 不能小于输入信号的长度。")
# 对输入信号进行快速傅里叶变换
X = fft(x)
if m % 2 == 0:
# 若原信号长度为偶数
# 计算中间位置
mid = m // 2
# 在频域中间插入零以进行插值
X_new = np.concatenate((X[:mid], np.zeros(n - m, dtype=X.dtype), X[mid:]))
else:
# 若原信号长度为奇数
# 计算中间位置
mid = (m + 1) // 2
# 在频域中间插入零以进行插值
X_new = np.concatenate((X[:mid], np.zeros(n - m, dtype=X.dtype), X[mid:]))
# 对零填充后的频域信号进行逆快速傅里叶变换
y = ifft(X_new) * (n / m)
# 返回插值后的实数信号（丢弃虚部）
return y.real
if isinstance(expr, tuple) and isinstance(expr[0], list):
if len(expr) == 3:
# 输入包含信号、插值后的长度和维度信息
matrix = sp.Matrix(expr[0])
n = int(expr[1])
dim = int(expr[2])
elif len(expr) == 2:
# 输入包含信号和插值后的长度，默认维度为 0
matrix = sp.Matrix(expr[0])
n = int(expr[1])
dim = 0
else:
# 输入的元组长度不符合要求，标记错误
error = True
if matrix is not None:
if matrix.shape[1] == 1:
# 如果矩阵是单列矩阵，将其展平为一维数组
x_np = np.array(matrix, dtype=float).ravel()
else:
# 否则，将矩阵转换为 NumPy 数组
x_np = np.array(matrix, dtype=float)
else:
# 输入不是有效的矩阵，标记错误
error = True
if dim == 0:
# 按整体进行插值
result = interpft_single(x_np, n)
elif dim == 1:
if x_np.ndim == 2:
# 按列进行插值
result = []
for col in x_np.T:
result.append(interpft_single(col, n))
result = np.array(result).T
else:
# 如果输入不是二维矩阵，按整体进行插值
result = interpft_single(x_np, n)
elif dim == 2:
if x_np.ndim == 2:
# 按行进行插值
result = []
for row in x_np:
result.append(interpft_single(row, n))
result = np.array(result)
else:
# 如果输入不是二维矩阵，按整体进行插值
result = interpft_single(x_np, n)
else:
# 输入的维度信息无效，标记错误
error = True
else:
# 输入不是元组，标记错误
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}" # 请增加注释
# 示例测试
if __name__ == "__main__":
# 示例1：偶数长度插值
x = [1, -1, 1, -1]
n = 6
y = interpolation_fft(f"{x}, {n}")
print(y)
# [ 1. -0.5 -0.5 1. -0.5 -0.5]
# 示例2：奇数长度插值
x = [1, 0, -1]
n = 5
y = interpolation_fft(f"{x}, {n}")
print(y)
# [ 1. 0.85810973 -0.46965902 -1.14837497 -0.24007574]
ndgrid格式的一维二维三维和N维网格数据的插值
Vq = interpn(X1,X2,...,Xn,V,Xq1,Xq2,...,Xqn) 使用线性插值返回 n 变量函数在特定查询点的插入值。结果始终穿过函数的原始采样。
X1,X2,...,Xn 包含样本点的坐标。V 包含各样本点处的对应函数值。Xq1,Xq2,...,Xqn 包含查询点的坐标。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.interpolate import RegularGridInterpolator
def parse_nd_array(param):
"""递归将 SymPy 嵌套结构转换为 n 维 numpy 数组"""
if isinstance(param, (sp.Tuple, list, tuple)):
return np.array(param, dtype=float)
elif param.is_Number:
return float(param.evalf())
else:
raise ValueError("Non-numeric element in array")
def interpolation_nd(input_str):
"""
对标 MATLAB 的 interpn 函数，实现 N 维网格数据插值。
参数:
input_str: 输入参数字符串，格式为：
"X1, X2, ..., Xn, V, Xq1, Xq2, ..., Xqn" -> 默认线性插值
"X1, X2, ..., Xn, V, Xq1, Xq2, ..., Xqn, method" -> 指定插值方法
"X1, X2, ..., Xn, V, Xq1, Xq2, ..., Xqn, method, extrap" -> 指定方法和外推值
返回:
插值结果数组，错误时返回错误信息字符串
示例:
>>> interpolation_nd("[0,1], [0,1], [[1,2],[3,4]], 0.5, 0.5")
array(2.5) # 二维线性插值示例
"""
try:
expr = sp.sympify(input_str, evaluate=False)
params = list(expr) if isinstance(expr, tuple) else [expr]
valid_methods = {'linear', 'nearest'}
m = len(params)
if m < 3:
raise ValueError("参数数量不足")
n = (m - 1) // 2
if 2 * n + 1 > m or m > 2 * n + 3:
raise ValueError(f"参数数量不匹配: m={m}, n={n}")
# 解析坐标轴
coords_axes = []
for i in range(n):
param = params[i]
matrix = sp.Matrix(param) if isinstance(param, list) else None
if matrix is not None:
arr = np.array(matrix, dtype=float).squeeze()
elif param.is_Number:
arr = np.array([float(param.evalf())])
else:
try:
arr = np.array([float(param.evalf())])
except:
raise ValueError(f"第 {i + 1} 个参数转换失败")
if arr.ndim != 1:
raise ValueError(f"X{i + 1} 必须是一维数组")
if not np.all(np.diff(arr) > 0):
raise ValueError(f"X{i + 1} 必须是严格递增的")
coords_axes.append(arr)
# 反转坐标轴以适应 scipy 的坐标顺序
reversed_coords = list(reversed(coords_axes))
expected_shape = tuple(len(ax) for ax in reversed_coords)
# 解析 V 矩阵
v_param = params[n]
try:
if isinstance(v_param, (sp.Tuple, list, tuple)):
V = parse_nd_array(v_param)
else:
raise ValueError("V 必须是多维数组")
if V.ndim != n:
raise ValueError(f"V 必须是 {n} 维，实际为 {V.ndim} 维")
if V.shape != expected_shape:
raise ValueError(f"V 形状不匹配: 期望 {expected_shape}，实际为 {V.shape}")
except Exception as e:
raise ValueError(f"V 解析失败: {str(e)}")
# 解析查询点
Xqs = []
for i in range(n):
param = params[n + 1 + i]
matrix = sp.Matrix(param) if isinstance(param, list) else None
if matrix is not None:
arr = np.array(matrix, dtype=float).squeeze()
elif param.is_Number:
arr = np.array([float(param.evalf())])
else:
try:
arr = np.array([float(param.evalf())])
except:
raise ValueError(f"Xq{i + 1} 转换失败")
Xqs.append(arr)
# 处理可选参数
method = 'linear'
extrap = np.nan
if m >= 2 * n + 2:
method = str(params[2 * n + 1])
if method not in valid_methods:
raise ValueError(f"无效的插值方法: {method}")
if m >= 2 * n + 3:
try:
extrap = float(params[2 * n + 2])
except:
raise ValueError("无效的外推值")
# 创建插值器
interpolator = RegularGridInterpolator(
reversed_coords,
V,
method=method,
bounds_error=False,
fill_value=extrap
)
# 生成查询网格
Xqs_reversed = list(reversed(Xqs))
mesh_arrays = np.meshgrid(*Xqs_reversed, indexing='ij')
query_points = np.column_stack([arr.ravel() for arr in mesh_arrays])
# 进行插值
Vq = interpolator(query_points)
Vq = Vq.reshape(mesh_arrays[0].shape)
# 转置以匹配原始查询点顺序
if n > 1:
Vq = np.transpose(Vq, axes=tuple(reversed(range(n))))
return Vq
except Exception as e:
return f"错误: {str(e)}"
result = interpolation_nd("[0,1], [0,1], [[1,2],[3,4]], 0.5, 0.5")
print(result)
# [[2.5]]
result = interpolation_nd("[0,1], [0,1], [0,1], [[[1,2],[3,4]],[[5,6],[7,8]]], 0.5, 0.5, 0.5")
print(result)
# [[[4.5]]]
矩阵求逆
Y = inv(X) 计算方阵 X 的逆矩阵
Y = inv(X,method) 按method方法分解矩阵,计算方阵X的逆矩阵.
X^(-1) 等效于 inv(X)
X — 输入矩阵,方阵.
method - 'GE', 'LU', 'ADJ', 'CH', 'LDL'
GE - 使用高斯消去法求逆.
LU - 使用LU分解计算求逆.
ADJ - 使用辅助矩阵和行列式求逆.
CH - 使用cholesky分解求逆.
LDL - 使用LDL分解求逆.
如果X是密集矩阵,默认method是GE，如果X是稀疏矩阵,默认method是LDL
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def inverse_matrix(input_str):
"""
对标MATLAB inv函数的矩阵求逆函数
参数:
input_str: 矩阵的字符串表示，可以是：
- 纯矩阵表达式（如 "[[1,2],[3,4]]"）
- 带方法的元组（如 "(Matrix([[1,2],[3,4]]), 'LU')"）
返回:
求逆后的矩阵（成功时）或错误信息字符串（失败时）
支持的方法参数:
'GE' - 高斯消元法（默认）
'LU' - LU分解法
'ADJ' - 伴随矩阵法
'CH' - Cholesky分解（需矩阵正定）
'LDL' - LDL分解（需矩阵厄米特）
示例:
>>> inverse_matrix("[[1, 0], [0, 1]]")
Matrix([[1, 0], [0, 1]])
>>> inverse_matrix("(Matrix([[1,2],[3,4]]), 'LU')")
Matrix([[-2, 1], [3/2, -1/2]])
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 处理带方法参数的元组输入（如 (Matrix, 'LU')）
if isinstance(expr, tuple) and len(expr) == 2:
# 尝试转换第一个元素为矩阵并验证方法参数
matrix = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
method = str(expr[1]) # 确保方法参数为字符串
# 验证矩阵转换是否成功和方法是否有效
if matrix is not None and method in ['GE', 'LU', 'ADJ', 'CH', 'LDL']:
result = matrix.inv(method=method).evalf()
else:
error = True
# 处理普通矩阵输入
else:
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is not None:
result = matrix.inv().evalf() # 使用默认方法（GE）
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：基本数值矩阵求逆")
print(inverse_matrix("[[1, 2], [3, 4]]"))
# Matrix([[-2.00000000000000, 1.00000000000000],
# [1.50000000000000, -0.500000000000000]])
反古德曼函数
invGudermannian(z) 返回z的逆古德曼函数。
z — 标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.optimize import fsolve
def inverse_gudermannian(input_str):
'''
返回z的逆古德曼函数。
:param input_str:
:return:
'''
try:
# 将输入字符串转换为SymPy表达式
# 将输入字符串转换为SymPy表达式
y = sp.sympify(input_str)
# 定义古德曼函数
def gudermannian(x):
return 2 * sp.atan(sp.exp(x)) - sp.pi / 2
# 如果y是符号性的（包含变量），则返回隐式符号方程
if y.free_symbols:
x = list(y.free_symbols)[0]
return sp.Eq(gudermannian(x), y)
# 如果y是数值，则使用fsolve进行数值反演
else:
# 定义要求解的方程
def equation(x):
return gudermannian(x) - y
# 使用fsolve定义数值求解器
def solve_equation_with_fsolve(equation, max_attempts=3, tolerance=1e-6):
last_solution = None
for _ in range(max_attempts):
# 随机生成初始猜测
initial_guess = np.random.uniform(-10, 10)
# 将符号方程转换为数值函数
equation_num = sp.lambdify(sp.Symbol('x'), equation(sp.Symbol('x')))
# 使用fsolve求解方程
solution = fsolve(equation_num, initial_guess)
# 检查方案的准确性
if np.isclose(equation_num(solution), 0, atol=tolerance):
return solution
last_solution = solution
return last_solution
# 求解方程并返回数值解
solution = solve_equation_with_fsolve(equation)
return solution[0]
except Exception as e:
return f"错误: {e}"
# --------------------------
# 示例代码
# --------------------------
if __name__ == "__main__":
# 示例1：基本数值矩阵求逆
print(inverse_gudermannian("0"))
# 6.083691555487232e-17
# 示例2：使用LU分解法
print(inverse_gudermannian("1.5"))
# 3.3406775427982867
# 示例3：单位矩阵求逆（结果应不变）
print(inverse_gudermannian("x"))
# Eq(2*atan(exp(x)) - pi/2, x)
希尔伯特矩阵的逆矩阵
H = invhilb(n) 生成确切希尔伯特矩阵的的逆矩阵.
n — 矩阵的阶次,非负整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def inverse_hilb_matrix(input_str):
"""
对标MATLAB invhilb函数，生成n×n希尔伯特矩阵的逆矩阵
参数:
input_str: 表示正整数的字符串
返回:
SymPy矩阵对象（成功时）或错误信息字符串（失败时）
示例:
>>> inverse_hilb_matrix("3")
Matrix([
[ 9, -36, 30],
[-36, 192,-180],
[ 30,-180, 180]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def hilbert_matrix_symbolic(n):
"""完全符号计算的实现方案"""
H = sp.zeros(n)
for i in range(n):
for j in range(n):
H[i, j] = 1 / (sp.Integer(i) + sp.Integer(j) + 1)
return H.inv()
if isinstance(expr, tuple):
error = True
elif expr.is_integer:
result = hilbert_matrix_symbolic(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# --------------------------
# 示例代码
# --------------------------
if __name__ == "__main__":
# 示例1：3阶希尔伯特逆矩阵
print("示例1：n=4")
print(inverse_hilb_matrix("4"))
# Matrix([[16, -120, 240, -140],
# [-120, 1200, -2700, 1680],
# [240, -2700, 6480, -4200],
# [-140, 1680, -4200, 2800]])
# 示例2：边界条件测试（n=1）
print("\n示例2：n=1")
print(inverse_hilb_matrix("1"))
# Matrix([[1]])
逆置换数组维度
A = ipermute(B,dimorder) 按照向量 dimorder 指定的顺序重新排列数组 B 的维度，使得 B = permute(A,dimorder)。换句话说，输入数组的第 i 个维度变为输出数组的维度 dimorder(i)。
B — 输入数组,向量,矩阵
dimorder — 维度顺序,行向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def inverse_permute_order(input_str):
"""
实现与MATLAB ipermute等效的功能，通过计算逆排列并调整转置顺序。
:param arr: 输入的数组，可以是sympy矩阵或numpy数组
:param order: 原permute操作使用的维度顺序（1-based，如MATLAB）
:return: 重排后的数组
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
if isinstance(expr, tuple) and len(expr) == 2:
if isinstance(expr[0], list) and isinstance(expr[1], list):
matrix = sp.Matrix(expr[0])
order = expr[1]
# 将MATLAB的1-based order转换为0-based
order = [o - 1 for o in order]
# 计算逆排列：inv_order使得应用后恢复原始维度
inv_order = np.argsort(order)
arr_np = np.array(matrix.tolist(), dtype=object)
result_np = np.transpose(arr_np, inv_order)
result = sp.Matrix(result_np.tolist())
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
def main():
# 示范使用sympy矩阵
sympy_result = inverse_permute_order("[[x,y],[z,x+y]],[2,1]")
print("\n使用sympy矩阵重排后的数组:")
print(sympy_result)
# Matrix([[x, z],
# [y, x + y]])
# 示范使用numpy数组
numpy_result = inverse_permute_order("[[1,2],[3,4]],[2,1]")
print("\n使用numpy数组重排后的数组:")
print(numpy_result)
# Matrix([[1, 3],
# [2, 4]])
if __name__ == "__main__":
main()
确定矩阵是否在指定带宽范围内
如果A是在指定的下带宽和上带宽范围内的矩阵,则tf = isbanded(A,lower,upper) 返回逻辑值 1 (true).否则,将返回逻辑值 0 (false).
A — 输入数组,数组
lower — 下带宽,非负整数标量
upper — 上带宽,非负整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_within_banded(input_str):
"""
对标 MATLAB 的 isbanded 函数，判断矩阵是否在指定带宽范围内。
参数:
input_str (str): 输入字符串，格式为 "矩阵, 下带宽, 上带宽"。
示例: "[[0,1,0], [2,0,3], [0,4,0]], 1, 0"
返回:
bool 或 str: 如果矩阵满足带宽要求返回 True，否则返回 False。
输入错误时返回错误信息字符串。
示例:
>>> result = is_within_banded("[[0,1,0], [2,0,3], [0,4,0]], 1, 0")
>>> print(result) # True
"""
try:
# 解析输入字符串为 SymPy 表达式
expr = sp.sympify(input_str)
# 验证输入是否为三元组 (矩阵, 下带宽, 上带宽)
if not (isinstance(expr, tuple) and len(expr) == 3):
return f"输入错误：需要形如 '矩阵, 下带宽, 上带宽' 的三元组，当前输入为 {input_str}"
# 提取矩阵、下带宽、上带宽
M_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
lower = expr[1]
upper = expr[2]
# 检查矩阵和带宽的有效性
if M_sym is None:
return "错误：输入的矩阵无效"
if not (lower.is_Integer and upper.is_Integer):
return "错误：下带宽和上带宽必须为整数"
lower = int(lower)
upper = int(upper)
if lower < 0 or upper < 0:
return "错误：带宽不能为负数"
# 检查矩阵是否为带状矩阵
rows, cols = M_sym.shape
for i in range(rows):
for j in range(cols):
# 关键逻辑：判断元素是否在带宽范围外且非零
if (i - j > lower) or (j - i > upper):
if M_sym[i, j] != 0:
return 0
return 1
except Exception as e:
return f"错误：{str(e)}"
input_str = "[[2,1,0], [3,2,1], [0,4,3]], 1, 1"
print(is_within_banded(input_str))
# 1
input_str = "[[1,0,5], [2,3,0], [0,4,6]], 1, 0"
print(is_within_banded(input_str))
# 0
确定哪些元素在指定范围内
TF=isbetween(A,lower,upper)确定输入数据中的哪些元素在由下限和上限定义的区间内,并返回与输入数据大小相同的逻辑数组. 默认情况下, 该间隔为封闭间隔. 当对应元素在指定范围内时,TF包含1(true),否则包含0(false).
A — 输入向量或矩阵
lower — 下限,标量或向量
upper — 上限,标量或向量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def is_between_range(input_str):
"""
判断矩阵元素是否在指定区间内 [lower, upper]。
返回 SymPy 矩阵，1 表示在区间内，0 表示不在。
参数格式示例：
"( [[1,2],[3,4]], [0,2], 5 )" 表示矩阵、下限向量、上限标量
"""
try:
expr = sp.sympify(input_str)
if not (isinstance(expr, tuple)) and len(expr) == 3:
return f"输入必须为三元组 (矩阵, lower, upper)，当前长度 {len(expr)}"
matrix = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
if matrix is None:
return "第一个元素必须为矩阵或可转换列表"
# 转换矩阵为 numpy 数组
if matrix.shape[1] == 1: # 列向量转一维
A = np.array(matrix).ravel().astype(float)
else:
A = np.array(matrix, dtype=float)
# 处理 lower
lower = expr[1]
if isinstance(lower, (list, sp.Matrix)):
lower = np.array(lower, dtype=float).ravel()
elif lower.is_number:
lower = float(lower)
else:
return f"lower 参数类型错误: {type(lower)}"
# 处理 upper
upper = expr[2]
if isinstance(upper, (list, sp.Matrix)):
upper = np.array(upper, dtype=float).ravel()
elif upper.is_number:
upper = float(upper)
else:
return f"upper 参数类型错误: {type(upper)}"
# 维度验证
if A.ndim == 2:
rows = A.shape[0]
if isinstance(lower, np.ndarray) and len(lower) != rows:
return f"lower 向量长度 {len(lower)} 与矩阵行数 {rows} 不匹配"
if isinstance(upper, np.ndarray) and len(upper) != rows:
return f"upper 向量长度 {len(upper)} 与矩阵行数 {rows} 不匹配"
# 核心判断逻辑
result = []
if A.ndim == 1: # 向量处理
l = lower if np.isscalar(lower) else lower[0]
u = upper if np.isscalar(upper) else upper[0]
bool_arr = np.logical_and(A >= l, A <= u)
result = sp.Matrix([int(x) for x in bool_arr])
else: # 矩阵处理
for i in range(A.shape[0]):
row = A[i]
l = lower[i] if isinstance(lower, np.ndarray) else lower
u = upper[i] if isinstance(upper, np.ndarray) else upper
bool_arr = np.logical_and(row >= l, row <= u)
result.append([int(x) for x in bool_arr])
result = sp.Matrix(result)
return result
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 示例 1: 向量 + 标量上下限
print(is_between_range("([1, 2, 3], 2, 4)"))
# Matrix([[0],
# [1],
# [1]])
# 示例 2: 矩阵 + 向量上下限
print(is_between_range("([[1,2], [3,4]], [0,3], [3,3])"))
# Matrix([[1, 1],
# [1, 0]])
确定矩阵是否为对角矩阵
如果A是在指定的下带宽和上带宽范围内的矩阵,则tf = isbanded(A,lower,upper) 返回逻辑值 1 (true).否则,将返回逻辑值 0 (false).
A — 输入数组,数组
lower — 下带宽,非负整数标量
upper — 上带宽,非负整数标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_diagonal_matrix(input_str):
"""
判断输入矩阵是否为对角矩阵 (对标 MATLAB 的 isdiag 函数)
参数:
input_str: 矩阵的字符串表示，支持以下格式：
"[[1, 0], [0, 2]]" - 二维列表格式
"[1, 0, 0]" - 一维向量格式
返回:
int: 1 表示是对角矩阵，0 表示不是
str: 错误信息字符串（输入无效时）
示例:
>>> is_diagonal_matrix("[[1, 0], [0, 2]]")
1
>>> is_diagonal_matrix("[[1, 2], [3, 4]]")
0
>>> is_diagonal_matrix("[[1, 0, 0], [0, 2, 0]]") # 非方阵
1
>>> is_diagonal_matrix("[[0, 0], [0, 0]]") # 零矩阵
1
>>> is_diagonal_matrix("invalid_matrix")
'Error: 语法错误'
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 检查输入是否为元组（无效类型）
if isinstance(expr, tuple):
return f"输入错误: 不支持元组类型 {input_str}"
# 转换为 SymPy 矩阵
M = sp.Matrix(expr) if isinstance(expr, list) else None
if M is None:
return f"输入错误: 无法转换为矩阵 {input_str}"
# 核心判断逻辑
return 1 if M.is_diagonal() else 0
except Exception as e:
return f"Error: {str(e)}"
# 测试用例
if __name__ == "__main__":
# 标准对角矩阵（方阵）
print(is_diagonal_matrix("[[1, 0], [0, 2]]"))
# 输出 1
# 非对角矩阵
print(is_diagonal_matrix("[[1, 2], [3, 4]]"))
# 输出 0
# 非方阵对角矩阵
print(is_diagonal_matrix("[[1, 0, 0], [0, 2, 0]]"))
# 输出 1
# 零矩阵
print(is_diagonal_matrix("[[0, 0], [0, 0]]"))
# 输出 1
# 列向量（非对角）
print(is_diagonal_matrix("[[1], [2], [3]]"))
# 输出 0
确定矩阵是埃尔米特矩阵还是斜埃尔米特矩阵
如果A是埃尔米特矩阵,则tf = ishermitian(A) 返回逻辑值1(true).否则,将返回逻辑值0(false).
tf = ishermitian(A,skewOption) 指定测试的类型.将skewOption指定为 "skew" 以确定A是否为斜埃尔米特矩阵.
A — 输入数组,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_hermitian_matrix(input_str):
"""
判断输入的矩阵是否为埃尔米特矩阵或斜埃尔米特矩阵。
参数:
input_str: 输入的字符串，表示矩阵或矩阵与选项的组合。例如：
- "Matrix([[1, I], [-I, 3]])" 检查埃尔米特矩阵。
- "(Matrix([[0, -1], [1, 0]]), skew)" 检查斜埃尔米特矩阵。
返回:
1 表示是，0 表示否，错误信息字符串表示输入有误。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
# 检查是否为 (矩阵, skew) 的形式
if len(expr) != 2:
error = True
else:
matrix_part, option_part = expr
M = sp.Matrix(matrix_part) if isinstance(matrix_part, list) else None
# 提取选项并去除可能的引号
option = str(option_part).strip("'\"")
if M is not None and option == "skew":
# 斜埃尔米特矩阵检查: M == -M.H
is_skew = M == -M.H
result = 1 if is_skew else 0
else:
error = True
else:
# 普通埃尔米特矩阵检查
M = sp.Matrix(expr) if isinstance(expr, list) else None
if M is not None:
is_herm = M.is_hermitian
result = 1 if is_herm else 0
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例1: 埃尔米特矩阵
print(is_hermitian_matrix("[[1, 2 + 3*I], [2 - 3*I, 4]]"))
# 输出: 1
# 示例2: 非埃尔米特矩阵
print(is_hermitian_matrix("[[1, 2], [3, 4]]"))
# 输出: 0
# 示例3: 斜埃尔米特矩阵
print(is_hermitian_matrix("[[0, -1], [1, 0]], skew"))
# 输出: 1
# 示例4: 非斜埃尔米特矩阵
print(is_hermitian_matrix("[[0, 2], [-2, 0]], skew"))
# 输出: 1
查找缺失值
TF = ismissing(A) 返回一个逻辑数组，该数组指示输入数据的哪些元素包含缺失值。TF 的大小与 A 的大小相同
A — 向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_nan_missing(input_str):
"""
对标 MATLAB 的 ismissing 函数，查找输入矩阵中的缺失值（NaN）。
参数:
input_str: 字符串形式的矩阵输入，例如 "[[1, nan], [3, 4]]"
返回:
如果输入合法，返回对应的 0-1 矩阵（1 表示 NaN）；
否则返回错误信息字符串。
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
def check_nan(x):
"""检查元素是否为 NaN，返回 1（是）或 0（否）"""
return 1 if x is sp.nan else 0
# 检查表达式类型
if isinstance(expr, tuple):
# 元组类型无法处理，标记错误
error = True
else:
# 尝试转换为矩阵
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is not None:
# 对每个元素应用 check_nan 函数
result = A.applyfunc(check_nan)
else:
# 无法转换为矩阵，标记错误
error = True
# 返回结果或错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
### 示例代码
if __name__ == "__main__":
# 示例 1: 输入列向量
input_str1 = "[nan, 2, 3]"
print(f"\n输入: {input_str1}")
print("输出:", is_nan_missing(input_str1))
# Matrix([[1],
# [0],
# [0]])
分离方程中的变量或表达式
isolate(eqn,expr)重新排列方程式eqn, 使表达式expr出现在左侧. 结果类似于求解eqn中的expr. 如果isolate无法isolate expr, 它会将所有包含expr的项移到左侧. 隔离的输出允许您使用subs从eqn中消除expr.
eqn — 输入方程,符号方程
expr — 要隔离的变量或表达式,符号变量,符号表达式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def isolate_expr_eqn(input_str):
"""
对标MATLAB的isolate函数，将方程中的指定变量分离到等式一侧。
参数:
input_str: 字符串形式的输入，格式为"(方程, 变量)"。
方程可以是等式（如Eq(a*x + b, 0)）或表达式（自动设为等于0）。
例如: "(a*x**2 + b*x + c, x)" 或 "(Eq(y, m*x + b), x)"
返回:
分离后的等式（如Eq(x, ...)），若失败返回错误信息字符串。
"""
try:
# 将输入字符串解析为SymPy表达式
expr = sp.sympify(input_str, evaluate=False)
# 验证是否为包含方程和变量的元组
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入格式错误: 需要形如'(方程, 变量)'的元组，实际输入: {input_str}"
eqn, var = expr[0], expr[1]
# 验证变量是否为符号
if not var.free_symbols:
return f"第二个元素必须是符号，实际类型: {type(var).__name__}"
# 将表达式自动转换为等式（如果不是等式）
if not isinstance(eqn, sp.Equality):
eqn = sp.Eq(eqn, 0) # 假设表达式等于0
# 求解方程并提取第一个解
solutions = sp.solve(eqn, var, dict=True)
if not solutions:
return f"无法分离变量{var}：方程无显式解"
# 构建等式结果
return sp.Eq(var, solutions[0][var])
except Exception as e:
return f"解析错误: {str(e)}"
### 示例代码 ###
if __name__ == "__main__":
# 示例1: 线性方程
input_str1 = "(a*y(t)**2+b*c==0,y(t))"
print(f"输入: {input_str1}")
print("输出:", isolate_expr_eqn(input_str1))
# Eq(y(t), -sqrt(-b*c/a))
# 示例2: 二次方程（自动返回第一个解）
input_str2 = "(a*x**2 + b*x + c, x)"
print(f"\n输入: {input_str2}")
print("输出:", isolate_expr_eqn(input_str2))
# Eq(x, (-b - sqrt(-4*a*c + b**2))/(2*a))
# 示例3: 显式等式输入
input_str3 = "(Eq(y, m*x + b), x)"
print(f"\n输入: {input_str3}")
print("输出:", isolate_expr_eqn(input_str3))
# Eq(x, (-b + y)/m)
查找数据中的离群值
TF = isoutlier(A,method) 返回一个逻辑数组,当在A的元素中检测到离群值时,该数组中与之对应的元素为true.
如果 A 是矩阵,则 isoutlier 分别对 A 的每列进行运算。
median — 离群值定义为与中位数相差超过三倍换算 MAD 的元素。换算 MAD 定义为 c*median(abs(A-median(A)))，其中 c=-1/(sqrt(2)*erfcinv(3/2))。
mean — 离群值定义为与均值相差超过三倍标准差的元素。此方法比 "median" 快，但没有它可靠。
quartiles — 离群值定义为比上四分位数 (75%) 大 1.5 个四分位差以上或比下四分位数 (25%) 小 1.5 个四分位差以上的元素。当 A 中的数据不是正态分布时，此方法很有用。
gesd — 使用广义极端 Student 化偏差检验检测离群值。此迭代方法与 "grubbs" 类似，但当有多个离群值互相遮盖时，此方法的执行效果更好。
A — 数值向量，矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import iqr, scoreatpercentile, t
def is_outlier_data(input_str):
"""
对标 MATLAB 的 isoutlier 函数，检测数据中的离群值。
参数:
input_str: 字符串形式的输入，可以是以下两种格式之一:
1. 直接输入矩阵，如 "[1, 2, 3]"
2. 元组格式 (矩阵, 方法)，如 "(Matrix([[1,2],[3,4]]), 'quartiles')"
支持的方法: "median", "mean", "quartiles", "grubbs", "gesd"
返回:
与输入矩阵形状相同的 0-1 矩阵，1 表示离群值；
若输入错误则返回错误信息字符串。
"""
try:
# 解析输入字符串为 SymPy 表达式
expr = sp.sympify(input_str, evaluate=False)
method = "median" # 默认方法
A = None
# 处理元组输入 (矩阵, 方法)
if isinstance(expr, tuple) and len(expr) == 2:
matrix_part = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
method_part = expr[1]
if matrix_part is None or not method_part.free_symbols:
return f"输入格式错误: {input_str}"
A = np.array(matrix_part.tolist(), dtype=float)
method = str(method_part)
# 处理单一矩阵输入
else:
matrix_part = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix_part is None:
return f"无法解析矩阵: {input_str}"
A = np.array(matrix_part.tolist(), dtype=float)
# 验证方法有效性
valid_methods = ["median", "mean", "quartiles", "grubbs", "gesd"]
if method.lower() not in valid_methods:
return f"未知方法: {method}，支持的方法: {valid_methods}"
def detect_outliers(data_col, method_name):
"""对单列数据执行离群值检测"""
n = len(data_col)
outliers = np.zeros(n, dtype=int)
# Median Absolute Deviation (MAD)
if method_name == "median":
median = np.median(data_col)
mad = np.median(np.abs(data_col - median))
if mad == 0:
mad = 1e-6 # 避免除以零
threshold = 3 * 1.4826 * mad # 常数缩放因子
outliers = (np.abs(data_col - median) > threshold).astype(int)
# Mean and Standard Deviation
elif method_name == "mean":
mean = np.mean(data_col)
std = np.std(data_col)
if std == 0:
std = 1e-6
threshold = 3 * std
outliers = (np.abs(data_col - mean) > threshold).astype(int)
# Interquartile Range (IQR)
elif method_name == "quartiles":
q1 = scoreatpercentile(data_col, 25)
q3 = scoreatpercentile(data_col, 75)
iqr_val = iqr(data_col)
lower = q1 - 1.5 * iqr_val
upper = q3 + 1.5 * iqr_val
outliers = ((data_col < lower) | (data_col > upper)).astype(int)
# Grubbs' Test (单离群值)
elif method_name == "grubbs":
if n < 3:
return outliers # 样本过少不检测
mean = np.mean(data_col)
std = np.std(data_col)
if std == 0:
return outliers
# 找到最大偏差点
max_idx = np.argmax(np.abs(data_col - mean))
G = np.abs(data_col[max_idx] - mean) / std
# 计算临界值
t_crit = t.ppf(1 - 0.05 / (2 * n), n - 2)
threshold = (n - 1) / np.sqrt(n) * np.sqrt(t_crit ** 2 / (n - 2 + t_crit ** 2))
if G > threshold:
outliers[max_idx] = 1
# Generalized ESD Test (多离群值)
elif method_name == "gesd":
alpha = 0.05
max_outliers = min(n // 2, 10) # 最多检测10个离群值
R = np.zeros(max_outliers)
lambdas = np.zeros(max_outliers)
for i in range(max_outliers):
mu = np.mean(data_col)
sigma = np.std(data_col)
if sigma == 0:
break
R[i] = np.max(np.abs(data_col - mu))
idx = np.argmax(np.abs(data_col - mu))
# 计算临界值
p = 1 - alpha / (2 * (n - i))
t_val = t.ppf(p, n - i - 2)
lambdas[i] = (n - i - 1) * t_val / np.sqrt((n - i - 2 + t_val ** 2) * (n - i))
if R[i] > lambdas[i] * sigma:
outliers[idx] = 1
data_col = np.delete(data_col, idx)
else:
break
return outliers
# 处理矩阵数据
if A.ndim == 1:
A = A.reshape(-1, 1)
n_rows, n_cols = A.shape
result = np.zeros((n_rows, n_cols), dtype=int)
# 逐列检测离群值
for col in range(n_cols):
result[:, col] = detect_outliers(A[:, col], method.lower())
# 转换为 SymPy 矩阵并保持原始形状
result_matrix = sp.Matrix(result.tolist())
if result_matrix.shape[1] == 1:
result_matrix = result_matrix.T
return result_matrix
except Exception as e:
return f"错误: {str(e)}"
### 示例代码 ###
if __name__ == "__main__":
# 示例1: 默认方法 (median)
input_str1 = "[1, 2, 3, 100]"
print(f"输入: {input_str1}")
print("输出:", is_outlier_data(input_str1))
# Matrix([[0, 0, 0, 1]])
# 示例2: 使用quartiles方法
input_str2 = "[[1, 2], [3, 4], [100, 200]], quartiles"
print(f"\n输入: {input_str2}")
print("输出:", is_outlier_data(input_str2))
# Matrix([[0, 0],
# [0, 0],
# [0, 0]])
# 示例3: Grubbs检验
input_str3 = "[1, 2, 3, 4, 100], grubbs"
print(f"\n输入: {input_str3}")
print("输出:", is_outlier_data(input_str3))
# Matrix([[0, 0, 0, 0, 1]])
等值面是三维数据分布中具有相等值的各点的三维曲面表示.
isosurface 函数通过连接空间的一个三维体内的常量值点来计算和绘制曲面.
[faces,verts] = isosurface(V,isovalue,points,X,Y,Z) 在单独的数组中返回面和顶点。
V — 三维体表达式,符号表达式
isovalue — 指定等值面值,标量
points — 采样点数
X — x 轴坐标, 向量 | 三维数组
Y — y 轴坐标数据, 向量 | 三维数组
Z — z 轴坐标数据, 向量 | 三维数组
faces — 面数据, 数组
verts — 顶点数据,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from skimage import measure
import regex, re
def extract_functions_and_assignments(input_string):
# 正则表达式匹配函数表达式
'''
function_pattern = r'\b\w+\((?:[^()]*|(?R))*\)'
functions = regex.findall(function_pattern, input_string)
# 正则表达式匹配变量赋值
assignment_pattern = r'\b\w+\s*=\s*\[[^]]*\]'
assignments = regex.findall(assignment_pattern, input_string)
'''
# 提取 assignments 列表
assignment_pattern = r'\b\w+\s*=\s*\[[^]]*\]'
assignments = regex.findall(assignment_pattern, input_string)
# 从 input_string 中删除每个 assignment
for assignment in assignments:
input_string = input_string.replace(assignment, '')
# 判断 input_string 是否是一个区间
if '[' in input_string and ']' in input_string:
# 如果是区间，使用 function_pattern 提取函数
function_pattern = r'\b\w+\((?:[^()]*|(?R))*\)'
functions = regex.findall(function_pattern, input_string)
else:
# 如果不是区间，按逗号分割为 functions 列表
functions = input_string.split(',')
# 移除首尾空格
functions = [func.strip() for func in functions]
# 移除空字符串
functions = list(filter(None, functions))
return functions, assignments
def parse_and_generate_all_var_values(input_params):
res_range = []
res_expr = []
for param in input_params:
var_str = param
var_str = var_str.replace(" ", "")
match = re.match(r'(\w+)\s*=\s*\[([^]]+)\]', var_str)
if not match:
return res_range, res_expr
var_name = match.group(1)
var_values_str = match.group(2)
try:
var_values = eval(var_values_str)
var_values = list(var_values) # 如果只有一个元素，转换为列表
var_values = var_values[:2]
var_values.insert(0, var_name)
res_range.append(var_values)
except (ValueError, NameError):
# Plot3D(sin(x^2+3y^2)/(0.1+r^2)+(x^2+5y^2)*exp(1-r^2)/2,r=[sqrt(x^2+y^2)])
# eval(sqrt(x^2+y^2)) raise except
res_expr.append([var_name, var_values_str])
return res_range, res_expr
def isosurface_data_volume(input_str):
"""
对标 MATLAB 的 isosurface 函数，从三维隐式函数生成等值面数据
参数:
input_str: 输入字符串，格式为 "表达式,等值面值"（例如："x**2 + y**2 + z**2, 1.0"）
grid_points: 每个维度的采样点数（默认30）
x_range: x轴采样范围（默认(-2,2)）
y_range: y轴采样范围（默认(-2,2)）
z_range: z轴采样范围（默认(-2,2)）
返回:
若成功则返回 (vertices, faces) 元组，否则返回错误信息字符串
"""
try:
functions, assignments = extract_functions_and_assignments(input_str)
# 解析表达式、等值面值和网格点数
if len(functions) < 1:
return "Error: Missing expression"
expr_str = functions[0]
expr = sp.sympify(expr_str)
isovalue = 1.0 # 默认等值面值
if len(functions) >= 2:
isovalue = float(functions[1])
grid_points = 30 # 默认网格点数
if len(functions) >= 3:
grid_points = int(functions[2])
# 解析变量范围
assignment_range, _ = parse_and_generate_all_var_values(assignments)
# 初始化变量范围
x_range = [-2.0, 2.0]
y_range = [-2.0, 2.0]
z_range = [-2.0, 2.0]
for var_info in assignment_range:
var_name = var_info[0]
if var_name == 'x' and len(var_info) >= 3:
x_range = [var_info[1], var_info[2]]
elif var_name == 'y' and len(var_info) >= 3:
y_range = [var_info[1], var_info[2]]
elif var_name == 'z' and len(var_info) >= 3:
z_range = [var_info[1], var_info[2]]
# 生成坐标数组
x = np.linspace(x_range[0], x_range[1], grid_points)
y = np.linspace(y_range[0], y_range[1], grid_points)
z = np.linspace(z_range[0], z_range[1], grid_points)
# 创建三维网格
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
# 将SymPy表达式转换为数值计算函数
x_sym, y_sym, z_sym = sp.symbols('x y z')
expr_func = sp.lambdify((x_sym, y_sym, z_sym), expr, 'numpy')
volume = expr_func(X, Y, Z)
# 计算等值面
verts, faces, _, _ = measure.marching_cubes(volume, isovalue)
return faces, verts
except ValueError as ve:
return f"数值错误: {str(ve)}"
except sp.SympifyError:
return "表达式解析失败"
except Exception as e:
return f"未知错误: {str(e)}"
### 示例代码 ###
if __name__ == "__main__":
# 示例1: 球面 (x² + y² + z² = 1)
print("示例1: 单位球面")
result = isosurface_data_volume("x**2+y**2+z**2, 1.0, 30, x=[-1.5, 1.5], y=[-1.5, 1.5], z=[-1.5, 1.5]")
if isinstance(result, tuple):
vertices, faces = result
print(vertices)
# [[ 2 1 0]
# [ 4 3 0]
# [ 4 0 1]
# ...
# [1798 1773 1775]
# [1798 1775 1799]
# [1799 1775 1766]]
print(faces)
# [[ 4.9652777 13. 14. ]
# [ 5. 12.826388 14. ]
# [ 5. 13. 13.652778 ]
# ...
# [24.034721 15. 16. ]
# [24.034721 16. 14. ]
# [24.034721 16. 15. ]]
print(f"顶点数: {vertices.shape[0]}, 面片数: {faces.shape[0]}")
# 顶点数: 3596, 面片数: 1800
# 示例2: 圆柱面 (x² + y² = 0.5)
print("\n示例2: 圆柱面")
result = isosurface_data_volume("x**2 + y**2 - 0.5, 0, 50, z=[-3, 3]")
if isinstance(result, tuple):
vertices, faces = result
print(f"顶点数: {vertices.shape[0]}, 面片数: {faces.shape[0]}")
# 顶点数: 7056, 面片数: 3600
确定哪些数组元素为质数
TF = isprime(X) 返回与 X 大小相同的逻辑数组. 如果 X(i) 为质数，则 TF(i) 的值为 true. 否则值为 false.
X — 输入值,标量,向量,实数数组,非负整数值
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_prime_number(input_str):
"""
对标 MATLAB 的 isprime 函数，判断输入中的元素是否为质数。
支持标量、矩阵和符号表达式。
参数:
input_str (str): 输入的数学表达式字符串，可以是标量、矩阵或符号。
返回:
- 对于数值元素：质数返回 1，非质数返回 0
- 对于符号元素：返回 None
- 错误时返回错误信息字符串
示例:
>>> is_prime_number("5")
1
>>> is_prime_number("[[2, 3], [5, 7]]")
Matrix([[1, 1], [1, 1]])
>>> is_prime_number("[[2, x], [4.5, 6]]")
Matrix([[1, None], [0, 0]])
"""
try:
expr = sp.sympify(input_str)
def evaluate_prime(x):
"""判断单个元素是否为质数"""
# 处理符号表达式
if x.free_symbols:
return None
# 处理非数值类型
if not x.is_number:
return None
try:
# 转换为整数（处理浮点数的整数情况）
num = int(sp.N(x))
if sp.N(x) != num: # 检查是否为非整数（如 4.5）
return 0
# 质数判断（质数定义需大于1）
return 1 if num > 1 and sp.isprime(num) else 0
except:
return 0 # 转换失败的情况（如复数）
# 处理标量输入
return evaluate_prime(expr)
except sp.SympifyError:
return f"表达式解析失败: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print("标量测试:")
print("5 →", is_prime_number("5"))
# 1
print("4 →", is_prime_number("4"))
# 0
print("2.0 →", is_prime_number("2.0"))
# 1
print("x →", is_prime_number("x"))
# None
确定矩阵是对称矩阵还是斜对称矩阵.
如果 A 是对称矩阵,则tf = issymmetric(A) 返回逻辑值 1 (true). 否则,将返回逻辑值 0 (false).
tf = issymmetric(A,skewOption)指定测试的类型. 将skewOption指定为"skew" 可确定A是否为斜对称矩阵
A — 输入数组,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import simplify
def is_symmetric_matrix(input_str):
"""
对标MATLAB issymmetric函数，判断矩阵的对称性
参数:
input_str - 矩阵描述字符串，支持两种格式：
1. 纯矩阵："[[0,1],[1,0]]"
2. 带类型标记："(Matrix([[0,1],[-1,0]]), 'skew')"
返回:
1 - 满足对称性要求
0 - 不满足对称性要求
错误信息字符串 - 输入不合法时
示例:
>>> is_symmetric_matrix("[[1,2],[2,1]]")
1
>>> is_symmetric_matrix("(Matrix([[0,1],[-1,0]]), 'skew')")
1
>>> is_symmetric_matrix("[[1,2],[3,4]]")
0
"""
try:
expr = sp.sympify(input_str, evaluate=False)
result = None
# 处理带选项的元组输入 (矩阵, 'skew')
if isinstance(expr, tuple) and len(expr) == 2:
M = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
option = str(expr[1]).lower()
if M is None or option not in ['skew']:
return f"输入错误: {input_str}"
# 斜对称矩阵检查：M^T = -M
# 使用符号化简处理浮点精度问题
is_skew = simplify(M.T + M).is_zero_matrix
result = 1 if is_skew else 0
# 处理普通矩阵输入
else:
M = sp.Matrix(expr) if isinstance(expr, list) else None
if M is None:
return f"输入错误: {input_str}"
# 对称矩阵检查：M^T = M
# 使用内置方法处理符号和数值矩阵
result = 1 if M.is_symmetric() else 0
return result
except sp.SympifyError:
return f"语法错误: 无法解析输入 '{input_str}'"
except Exception as e:
return f"运行时错误: {str(e)}"
# --------------------------
# 测试用例
# --------------------------
if __name__ == "__main__":
# 测试对称矩阵
print("对称矩阵测试:")
print(is_symmetric_matrix("[[1, 2], [2, 1]]"))
# 1
print(is_symmetric_matrix("[[a, b], [b, a]]"))
# 1
# 测试斜对称矩阵
print("\n斜对称矩阵测试:")
print(is_symmetric_matrix("[[0, 2], [-2, 0]], skew"))
# 1
print(is_symmetric_matrix("[[0, 1.0], [-1.0001, 0]], skew"))
# 0
# 测试非对称矩阵
print("\n非对称矩阵测试:")
print(is_symmetric_matrix("[[1, 2], [3, 4]]"))
# 0
确定矩阵是否为下三角矩阵.
如果A是下三角矩阵,则tf = istril(A) 返回逻辑值1 (true). 否则,将返回逻辑值0 (false).
A — 输入数组,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def is_lower_triangular(input_str):
"""
判断输入是否为下三角矩阵（对标MATLAB istril函数）
参数:
input_str: 输入字符串，可以是矩阵、标量或符号表达式
返回:
1: 是下三角矩阵
0: 不是下三角矩阵
str: 错误信息
示例:
>>> is_lower_triangular("[[1,0],[2,3]]")
1
>>> is_lower_triangular("[[1,2],[3,4]]")
0
>>> is_lower_triangular("5")
1
>>> is_lower_triangular("[[1,2],[3,4,5]]")
'输入错误: [[1,2],[3,4,5]]'
"""
try:
expr = sp.sympify(input_str)
# 尝试转换为矩阵
if isinstance(expr, list):
M = sp.Matrix(expr)
# 非方阵直接返回0
if M.is_square and M.is_lower:
return 1
else:
return 0
# 处理标量和符号（视为1x1矩阵）
if isinstance(expr, sp.Basic) and expr.is_Atom:
return 1
return f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例测试
if __name__ == "__main__":
test_cases = [
"[[1,0],[2,3]]",
# 1
"[[1,2],[3,4]]",
# 0
"5",
# 1
"a",
# 1
"[1,2]",
# 0
"[[1,2,3],[4,5,6]]"
# 0
]
for case in test_cases:
print(f"输入: {case:<20} 结果: {is_lower_triangular(case)}")
确定矩阵是否为上三角矩阵.
如果A是上三角矩阵,则tf = istril(A) 返回逻辑值1 (true). 否则,将返回逻辑值0 (false).
A — 输入数组,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
def is_upper_triangular(input_str):
"""
判断输入是否为上三角矩阵（对标MATLAB istriu函数）
参数:
input_str: 输入字符串，可以是矩阵、标量或符号表达式
返回:
1: 是上三角矩阵
0: 不是上三角矩阵
str: 错误信息
示例:
>>> is_upper_triangular("[[1,0],[2,3]]")
1
>>> is_upper_triangular("[[1,2],[3,4]]")
0
>>> is_upper_triangular("5")
1
>>> is_upper_triangular("[[1,2],[3,4,5]]")
'输入错误: [[1,2],[3,4,5]]'
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 尝试转换为矩阵
if isinstance(expr, list):
M = sp.Matrix(expr)
# 非方阵直接返回0
if M.is_square and M.is_upper:
return 1
else:
return 0
# 处理标量和符号（视为1x1矩阵）
if isinstance(expr, sp.Basic) and expr.is_Atom:
return 1
return f"输入错误: {input_str}"
except Exception as e:
return f"Error: {e}"
# 示例测试
if __name__ == "__main__":
test_cases = [
"[[17,24,1,8,15],[0,5,7,14,16],[0,0,13,20,22],[0,0,0,21,3],[0,0,0,0,9]]",
# 1
"[[1,2],[3,4]]",
# 0
"5",
# 1
"a",
# 1
"[1,2]",
# 0
"[[1,2,3],[4,5,6]]"
# 0
]
for case in test_cases:
print(f"输入: {case:<20} 结果: {is_upper_triangular(case)}")
逆z变换
iztrans(F) 返回 F 的逆 Z 变换. 默认情况下, 独立变量为 z, 变换变量为 n. 如果 F 不包含 z, iztrans 将使用函数 symvar.
F — 输入, 符号表达式, 符号函数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import lcapy as lc
def inverse_z_transform(input_str):
"""执行逆Z变换，支持矩阵输入
参数:
input_str (str): 输入表达式字符串或矩阵描述
返回:
SymPy表达式/矩阵 或 错误信息字符串"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
# 转换为Lcapy表达式以使用其符号运算功能
error = False
result = None
def eval_inv_ztrans(ele):
"""执行单个表达式的逆Z变换"""
f_expr = lc.expr(str(ele))
if f_expr.free_symbols:
return f_expr.IZT() # 调用Lcapy的逆Z变换方法
else:
return None # 常数项无需变换
# 矩阵处理分支
if isinstance(expr, list):
matrix = sp.Matrix(expr)
rows, cols = matrix.shape
# 对矩阵中每个元素递归执行逆Z变换
result = sp.Matrix(rows, cols, lambda i, j: eval_inv_ztrans(matrix[i, j]))
elif expr.free_symbols:
# 单表达式处理
result = eval_inv_ztrans(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：单表达式逆Z变换")
expr_str = "1/(1 - 1/z)"
result = inverse_z_transform(expr_str)
print(f"输入: {expr_str}")
print(f"输出: {result}")
# Piecewise((1, n >= 0))
print("\n" + "=" * 50 + "\n")
# 示例2：矩阵变换
print("示例2：矩阵逆Z变换")
matrix_input = "[2*z/(z-2)**2, 1/(a*z)]"
result = inverse_z_transform(matrix_input)
print("\n逆Z变换结果:")
print(result)
# Matrix([[DiscreteTimeDomainExpression(Piecewise((2*2**(n - 1)*n, n >= 0)))],
# [DiscreteTimeDomainExpression(UnitImpulse(n - 1)/a)]])