广义拉盖尔函数和拉盖尔多项式
如果n是非负整数,laguerreL(n,x)返回n次拉盖尔多项式.当n不是非负整数时.laguerreL返回拉盖尔函数.
如果n是非负整数,则laguerreL(n,a,x) 返回n次广义拉盖尔多项式.
n — 多项式的次数, 数,向量,矩阵,符号数,符号向量,符号矩阵,符号函数
x — 输入, 数,向量,矩阵,符号数,符号向量,符号矩阵,符号函数
a — 输入, 数,向量,矩阵,符号数,符号向量,符号矩阵,符号函数
示例1: 量子力学 - 氢原子波函数径向部分
氢原子波函数的径向部分包含广义拉盖尔多项式
n=2, l=1 (a = 2l+1 = 3) 的径向波函数在r=1处的值
result = laguerreL(1, 3, 1.0) # n=1, a=3, x=1.0
print(f" 输出：{result}")
#输出：3
#说明：这是氢原子2p轨道径向波函数的一部分
示例2: 量子谐振子 - 能量本征态
一维谐振子的波函数包含拉盖尔多项式
n=3 的谐振子波函数在x=0.5处的值
result = laguerreL(3, 0.5) # 普通拉盖尔多项式
print(f" 输出：{result}")
#输出：-0.145833333333333
#说明：用于计算量子谐振子的概率振幅
示例3: 数值分析 - 高斯-拉盖尔积分
计算积分 ∫₀^∞ e^{-x} f(x) dx 的节点和权重
3点高斯-拉盖尔求积的节点对应拉盖尔多项式的根
result = laguerreL(3, 0.5) # L₃(x)在x=0.5处的值
print(f" 输出：{result}")
#输出：-0.145833333333333
#说明：用于构造数值积分公式
示例4: 概率论 - 拉盖尔分布
在排队论和可靠性理论中的应用
计算拉盖尔多项式在特定点的值
result = laguerreL(2, 1, 0.8) # 广义拉盖尔多项式
print(f" 输出：{result}")
#输出：0.92
#说明：用于某些特殊概率分布的矩生成函数
示例5: 电磁学 - 激光腔模式
拉盖尔-高斯光束的横向模式分布
LG₀₁模式 (径向指数p=0,角向指数l=1)
result = laguerreL(0, 1, 0.5) # n=0, a=1, x=0.5
print(f" 输出：{result}")
#输出：1.0
#说明：描述激光束的横向强度分布
示例6: 热传导 - 柱坐标系中的温度分布
某些边界条件下的热传导方程解包含拉盖尔多项式
result = laguerreL(2, 0, 1.2) # 普通拉盖尔多项式
print(f" 输出：{result}")
#输出：-0.68
#说明：用于分析圆柱体内的温度分布
示例7: 量子化学 - Slater型轨道
在量子化学计算中，拉盖尔多项式用于构造原子轨道
result = laguerreL(1, 2, 0.3) # 广义拉盖尔多项式
print(f" 输出：{result}")
#输出：2.7
#说明：用于构建分子轨道的基函数
示例8: 信号处理 - 正交滤波器组
拉盖尔多项式用于设计某些类型的数字滤波器
result = laguerreL(4, 0.2) # 普通拉盖尔多项式
print(f" 输出：{result}")
#输出：0.314733333333333
#说明：在信号处理中构造正交基函数
示例9: 矩阵计算示例 - 多参数同时计算
计算不同阶数拉盖尔多项式在多个点的值
result = laguerreL([[0,1,2], [1,2,3]], [[0.1,0.2,0.3], [0.4,0.5,0.6]])
print(f" 输出：{result}")
#输出：[[1.0, 0.8, 0.445],
[0.6, 0.125, -0.296]]
#说明：批量计算提高效率
示例10: 符号计算应用 - 一般化表达式
使用符号变量的拉盖尔多项式
result = laguerreL(n, a, x)
print(f" 输出：{result}")
#输出：assoc_laguerre(n, a, x)
#说明：保持符号形式，用于解析推导
示例11: 机器学习 - 径向基函数网络
拉盖尔多项式作为径向基函数的扩展
result = laguerreL(2, 0, 1.5) # 用于函数逼近
print(f" 输出：{result}")
#输出：-0.875
#说明：在神经网络中作为激活函数
示例12: 量子信息 - 连续变量量子态
在量子光学中用于描述光场的量子态
result = laguerreL(3, 1, 0.7) # 用于Wigner函数的计算
print(f" 输出：{result}")
#输出：0.722833333333334
#说明：用于量子态层析和重构
示例13: 数学物理 - 特殊函数关系
展示拉盖尔多项式与其他特殊函数的关系
验证 L₀^a(x) = 1
result = laguerreL(0, 2, 1.5)
print(f" 输出：{result}")
#输出：1.0
#说明：零阶广义拉盖尔多项式恒为1
示例14: 渐近行为研究 - 大参数行为:
研究拉盖尔多项式在大n或大x时的行为
result = laguerreL(10, 0, 5.0) # 较高阶多项式
print(f" 输出：{result}")
#输出：1.75627617945334
#说明：用于研究特殊函数的渐近性质
示例15: 教育演示 - 正交多项式演示
展示拉盖尔多项式的正交性质
result = laguerreL([0,1,2], [0.5,1.0,1.5])
print(f" 输出：{result}")
#输出：[1.0, 0, -0.875]
#说明：演示加权正交性 ∫₀^∞ x^a e^{-x} L_m^a(x) L_n^a(x) dx = 0 (m≠n)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import eval_genlaguerre
def laguerre_polynomial(input_str):
"""
计算广义或普通 Laguerre 多项式，对标 MATLAB 的 laguerreL 函数
参数:
input_str: 字符串形式的输入，支持以下格式：
- 两个参数: "(n, x)" 对应普通 Laguerre 多项式
- 三个参数: "(n, a, x)" 对应广义 Laguerre 多项式
参数可以是标量或矩阵
返回:
SymPy 表达式/矩阵 或 错误信息字符串
示例:
>>> laguerre_polynomial("(2, 0.5)")
0.125
>>> laguerre_polynomial("(2, 1, 0.5)")
0.375
>>> laguerre_polynomial("([[1,2], [3,4]], [[0,1], [2,3]])")
Matrix([
[1, 0.5],
[0, -0.5]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 参数解析逻辑
if isinstance(expr, tuple):
# 处理广义 Laguerre 多项式 (n, a, x)
if len(expr) == 3:
if all(e.is_number for e in expr):
n = int(expr[0])
a = float(expr[1])
x = complex(expr[2])
result = eval_genlaguerre(n, a, x)
elif any(e.free_symbols for e in expr):
result = sp.assoc_laguerre(*expr).doit()
else:
error = True
elif len(expr) == 2:
if all(e.is_number for e in expr):
n = int(expr[0])
a = 0
x = complex(expr[1])
result = eval_genlaguerre(n, a, x)
elif any(e.free_symbols for e in expr):
result = sp.laguerre(*expr).doit()
else:
error = True
else:
error = True
else:
error = True
return result if not error else f"错误：输入必须是参数元组"
except Exception as e:
return f"错误：{str(e)}"
# 示例测试
if __name__ == "__main__":
# 普通 Laguerre 测试
print("普通 Laguerre 多项式测试:")
print(laguerre_polynomial("0, 5"))
# (1+0j)
print(laguerre_polynomial("1, 2"))
# (-1+0j)
print(laguerre_polynomial("x, 0.5"))
# laguerre(x, 0.5)
# 广义 Laguerre 测试
print("\n广义 Laguerre 多项式测试:")
print(laguerre_polynomial("(x, 0, 5)"))
# laguerre(x, 0.5)
print(laguerre_polynomial("(2, 1, 0.5)"))
# (1.625+0j)
朗伯W函数
lambertw(x)返回Lambert W函数的主分支.此语法等效于lambertw(0,x).
lambertw(k，x)是Lambert W函数的第k个分支.仅当k=0或k=-1时.此语法才返回实值.
x —— 输入,数字,矢量,矩阵,数组,符号数,符号变量,符号数组,符号函数,符号表达式
k —— Lambert W函数的分支,整数,整数的矢量或矩阵,符号整数,符号矢量或整数的矩阵
示例1: 数学 - 超越方程求解
print(" 解方程: x * e^x = 2")
print(" 解为: x = W(2)")
solution = 2
result = lambertw(solution)
print(f" 输出：{result}")
print(" 验证：x * e^x =", result * np.exp(result))
#解方程: x * e^x = 2
#解为: x = W(2)
#输出：0.8526055020137254
#验证：x * e^x = 1.9999999999999998
#说明：朗伯W函数是方程 x * e^x = a 的解
示例2: 物理学 - 理想玻色气体
理想玻色气体的临界温度计算
临界密度方程涉及朗伯W函数
bose_gas = 1.0 # 简化参数
result = lambertw(bose_gas)
print(f" 输出：{result}")
#输出：0.5671432904097838
#说明：在玻色-爱因斯坦凝聚理论中用于计算临界温度
示例3: 电子工程 - 二极管电流-电压关系
某些二极管模型中的显式解
diode_param = 0.5 # 归一化参数
result = lambertw(diode_param)
print(f" 输出：{result}")
#输出：0.35173371124919584
#说明：用于从隐式方程中解析求解二极管电流
示例4: 金融数学 - Black-76模型
某些期权定价模型的解析解
option_param = 0.3 # 隐含波动率相关参数
result = lambertw(option_param)
print(f" 输出：{result}")
#输出：0.23675531078855933
#说明：在金融衍生品定价中用于解析求解
示例5: 组合数学 - Cayley公式推广
带权树的枚举问题
tree_count = 1.5 # 参数化情况
result = lambertw(tree_count)
print(f" 输出：{result}")
#输出：0.7258613577662263
#说明：在组合计数中用于生成函数的反函数
示例6: 计算机科学 - 哈希表性能
分析哈希表冲突概率
hash_param = 0.8 # 负载因子相关
result = lambertw(hash_param)
print(f" 输出：{result}")
#输出：0.4900678588015799
#说明：在算法分析中用于求解某些递归关系
示例7: 控制理论 - 系统稳定性
分析含时滞的反馈系统稳定性
result = lambertw(-1, -0.5) # 使用分支k=-1
print(f" 输出：{result}")
#输出：-0.7940236323446893-0.7701117505103791j
#说明：时滞微分方程的特征值计算
示例8: 热力学 - 平均场理论
铁磁相变的临界行为
thermo_param = 0.6
result = lambertw(thermo_param)
print(f" 输出：{result}")
#输出：0.4015636367870726
#说明：在统计物理中用于求解自洽方程
示例9: 矩阵计算 - 多参数同时计算
计算不同参数下的朗伯W函数值
result = lambertw([[0.1, 0.5], [1.0, 2.0]])
print(f" 输出：{result}")
#输出：[[0.09127653, 0.35173371]
[0.56714329, 0.8526055]]
#说明：批量计算提高效率
示例10: 符号计算应用 - 一般化表达式
使用符号变量的朗伯W函数
symbolic_input = z
result = lambertw(symbolic_input)
print(f" 输出：{result}")
#输出：LambertW(z)
#说明：保持符号形式，用于解析推导
示例11: 通信理论 - 高斯信道容量
某些信道容量计算涉及朗伯W函数
channel_capacity = 2.0 # SNR相关参数
result = lambertw(channel_capacity)
print(f" 输出：{result}")
#输出：0.8526055020137254
#明：在信息论中用于信道容量分析
示例12: 生物学 - 种群动力学
某些非线性增长方程的解析解
population_param = 0.8
result = lambertw(population_param)
print(f" 输出：{result}")
#输出：0.4900678588015799
#说明：在生态学中用于种群动态建模
示例13: 化学动力学 - 复杂反应网络
某些化学反应速率方程的解析解
kinetics_param = 0.4
result = lambertw(kinetics_param)
print(f" 输出：{result}")
#输出：0.29716775067313855
#说明：在化学动力学中用于求解速率方程
示例14: 经济学 - 最优化问题
某些经济模型中的最优化条件
econ_param = 0.7
result = lambertw(econ_param)
print(f" 输出：{result}")
#输出：0.447470259269655
#说明：在微观经济学中用于求解最优化问题
示例15: 复变函数 - 多分支行为
展示朗伯W函数在不同分支的值
for k in range(-2, 3):
branch_input = (k, 0.5)
print(f" 分支 k={k}, 输入：{branch_input}")
result = lambertw(branch_input)
print(f" 输出：{result}")
#分支 k=-2, 输入：(-2, 0.5)
#输出：-3.1049770718920247-10.71348331130125j
#分支 k=-1, 输入：(-1, 0.5)
#输出：-2.259158898533606-4.220960969266197j
#分支 k=0, 输入：(0, 0.5)
#输出：0.35173371124919584
#分支 k=1, 输入：(1, 0.5)
#输出：-2.259158898533606+4.220960969266197j
#分支 k=2, 输入：(2, 0.5)
#输出：-3.1049770718920247+10.71348331130125j
#说明：朗伯W函数是多值函数，有无限多个分支
示例16: 工程应用 - 光伏系统
太阳能电池的单二极管模型")
pv_param = 1.2 # 与光照和温度相关的参数
result = lambertw(pv_param)
print(f" 输出：{result}")
#输出：0.6355640163648698
#说明：用于从隐式I-V关系中解析求解电流
示例17: 数值方法验证 - 验证恒等式
验证 W(z) * exp(W(z)) = z
test_values = [0.1, 0.5, 1.0, 2.0]
for val in test_values:
w_val = lambertw(val)
verification = w_val * np.exp(w_val)
print(f" z={val}: W(z)={w_val}, W(z)*exp(W(z))={verification}")
#z=0.1: W(z)=0.09127652716086226, W(z)*exp(W(z))=0.1
#z=0.5: W(z)=0.35173371124919584, W(z)*exp(W(z))=0.5
#z=1.0: W(z)=0.5671432904097838, W(z)*exp(W(z))=1
#z=2.0: W(z)=0.8526055020137254, W(z)*exp(W(z))=1.9999999999999998
#说明：验证朗伯W函数的基本性质
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import lambertw
def lambert_w_value(input_str):
"""
对标 MATLAB 的 lambertw 函数，计算 Lambert W 函数的值，支持标量和矩阵输入。
参数:
input_str: 输入的字符串，可以是数值、矩阵或元组（分支参数和变量参数）。
返回:
计算结果（数值、矩阵或符号表达式），或错误信息。
示例:
>>> lambert_w_value("0.5")
0.3517337112491958
>>> lambert_w_value("(-1, 0.5)")
-1.756431208626170
>>> lambert_w_value("[[1, 2], [3, 4]]")
Matrix([
[0.567143290409784, 0.852605502013725],
[1.04990889496404, 1.20216787319705]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理元组输入 (分支k, 变量x)
if isinstance(expr, tuple) and len(expr) == 2:
k_part = expr[0]
x_part = expr[1]
if k_part.is_integer and x_part.is_number:
k = int(k_part)
x = complex(x_part)
result = lambertw(x, k=k)
elif k_part.is_integer and x_part.free_symbols:
result = sp.LambertW(x_part, k_part)
else:
error = True
elif expr.is_number:
k = 0
x = complex(expr)
result = lambertw(x, k=k)
elif expr.free_symbols:
result = sp.LambertW(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: 计算主分支 W0(0.5)
print("示例1:", lambert_w_value("0.5"))
# (0.35173371124919584+0j)
# 示例2: 计算分支 k=-1 的 W_{-1}(0.5)
print("示例2:", lambert_w_value("(-1, 0.5)"))
# (-2.259158898533606-4.220960969266197j)
拉普拉斯变换
F = laplace(f) 返回 f 的拉普拉斯变换. 默认情况下, 独立变量为t, 变换变量为s.
F = laplace(f,transVar) 使用变换变量 transVar 而不是 s.
F = laplace(f,var,transVar) 分别使用独立变量 var 和变换变量 transVar 而不是 t 和 s.
f — 输入, 符号表达式,符号函数,符号向量,符号矩阵
var — 独立变量,t(默认),符号变量. 独立变量, 指定为符号变量. 此变量通常称为“时间变量”或“空间变量”. 如果您未指定变量, 则默认情况下, 拉普拉斯使用t. 如果f不包含t, 则拉普拉斯使用函数 symvar 来确定独立变量.
transVar — 变换变量, s（默认）| z , 符号变量 , 符号表达式 , 符号向量 , 符号矩阵 变换变量,指定为符号变量、表达式、向量或矩阵. 此变量通常称为“复频率变量”. 如果您未指定变量,则默认情况下,拉普拉斯使用s. 如果s是f的独立变量, 则拉普拉斯使用z.
示例1: 电路分析 - RC电路响应
RC电路的电压响应 v(t) = 1 - exp(-t/RC)
rc_circuit = 1 - exp(-t)
result = laplace(rc_circuit)
print(f" 输出：{result}")
#输出：-1/(s + 1) + 1/s
#说明：在s域分析电路瞬态响应
示例2: 控制系统 - 二阶系统传递函数
质量-弹簧-阻尼系统的微分方程
spring_mass = sin(2*t) # 简谐激励
result = laplace(spring_mass)
print(f" 输出：{result}")
#输出：2/(s**2 + 4)
#说明：将时域微分方程转换为s域代数方程
示例3: 信号处理 - 理想低通滤波器
矩形脉冲的拉普拉斯变换
rect_pulse = Heaviside(t) - Heaviside(t-1) # 单位矩形脉冲
result = laplace(rect_pulse)
print(f" 输出：{result}")
#输出：1/s - exp(-s)/s
#说明：用于分析滤波器的频率响应
示例4: 机械振动 - 单位脉冲响应
狄拉克δ函数的拉普拉斯变换
delta_func = DiracDelta(t)
result = laplace(delta_func)
#输出：1
#说明：冲击响应的s域表示恒为1
示例5: 热传导 - 一维热传导方程
瞬时点热源的响应
heat_source = exp(-a*t) # 衰减热源
result = laplace(heat_source)
print(f" 输出：{result}")
#输出：1/(a + s)
#说明：将偏微分方程转换为常微分方程
示例6: 通信系统 - 调制信号
载波信号的拉普拉斯变换
am_signal = cos(w*t)*exp(-a*t) # 衰减调幅波
result = laplace(am_signal)
print(f" 输出：{result}")
#输出：(a + s)/(w**2 + (a + s)**2)
#说明：分析调制信号的频谱特性
示例7: 电力系统 - 短路电流
电力系统暂态过程分析
short_circuit = t*exp(-t) # 短路电流近似
result = laplace(short_circuit)
print(f" 输出：{result}")
#输出：(s + 1)**(-2)
#说明：电力系统故障分析
示例8: 经济学 - 指数增长模型
资本积累的连续时间模型
growth_model = exp(r*t) # 指数增长
result = laplace(growth_model)
print(f" 输出：{result}")
#输出：1/(-r + s)
#说明：经济动力学系统的稳定性分析
示例9: 矩阵系统 - 状态空间模型
多输入多输出系统的变换
system_matrix = [[t, t**2], [exp(-t), sin(t)]]
result = laplace(system_matrix)
print(f" 输出：{result}")
#输出：[[s**(-2), 2/s**3],
[1/(s + 1), 1/(s**2 + 1)]]
#说明：多变量系统的s域分析
示例10: 符号计算 - 参数化系统
带参数的拉普拉斯变换
symbolic_system = exp(-a*t)*sin(w*t)
result = laplace(symbolic_system)
print(f" 输出：{result}")
#输出：w/(w**2 + (a + s)**2)
#说明：保持参数形式，用于解析研究
示例11: 工程综合 - PID控制器
PID控制器的s域表示
pid_components = [1, 1/t, t] # P, I, D项
for i, comp in enumerate(pid_components):
print(f" {['比例', '积分', '微分'][i]}项: {comp}")
result = laplace(comp)
print(f" 变换结果: {result}")
#比例项: 1
#变换结果: 1/s
#积分项: 1/t
#变换结果: LaplaceTransform(1/t, t, s)
#微分项: t
#变换结果: s**(-2)
#说明：控制器设计与稳定性分析
示例12: 验证拉普拉斯变换性质
线性性质: L{a*f(t) + b*g(t)} = a*L{f(t)} + b*L{g(t)}
f1, f2 = exp(-t), sin(t)
a, b = 2, 3
result = laplace(a*f1, b*f2)
print(f" 输出：{result}")
#输出：3/(s**2 + 1) + 2/(s + 1)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def laplace_transform(input_str):
"""
对标 MATLAB 的 laplace 函数，计算拉普拉斯变换
参数:
input_str: 输入字符串，可以是：
1. 单个表达式 (自动检测变量)
2. 元组 (f, transVar)
3. 元组 (f, var, transVar)
4. 矩阵
返回:
拉普拉斯变换结果 (数值、矩阵或符号表达式)，或错误信息
示例:
>>> laplace_transform("exp(-t)")
1/(s + 1)
>>> laplace_transform("(exp(-t), z)")
1/(z + 1)
>>> laplace_transform("(sin(a*x), x, s)")
a/(a**2 + s**2)
>>> laplace_transform("[[t, t**2], [t**3, t**4]]")
Matrix([
[ 1/s**2, 2/s**3],
[6/s**4, 24/s**5]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def compute_laplace(f, var=None, trans_var=None):
"""
核心计算函数
"""
# 自动检测变量
if var is None:
# 优先尝试常用变量 t
t = sp.symbols('t')
var = t if f.has(t) else next(iter(f.free_symbols), t)
# 自动生成变换变量
if trans_var is None:
trans_var = sp.symbols('s') if var == sp.symbols('t') else sp.symbols('z')
# 执行拉普拉斯变换并忽略收敛条件
return sp.laplace_transform(f, var, trans_var, noconds=True)
# 处理元组输入 (f, var, trans_var) 或 (f, trans_var)
if isinstance(expr, tuple):
params = list(expr)
n_params = len(params)
if n_params not in (2, 3):
return f"错误：需要 2 或 3 个参数，当前收到 {n_params} 个"
# 解析参数
f_part = params[0]
var_part = params[1] if n_params == 3 else None
trans_var_part = params[-1]
# 标量计算
result = compute_laplace(
f=f_part,
var=var_part,
trans_var=trans_var_part
)
# 处理标量输入
elif expr.free_symbols:
result = compute_laplace(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:", laplace_transform("exp(-t)"))
# 1/(s + 1)
# 测试2: 指定变换变量
print("测试2:", laplace_transform("(exp(-t), z)"))
# 1/(z + 1)
# 测试3: 指定变量和变换变量
print("测试3:", laplace_transform("(sin(a*x), x, s)"))
# a/(a**2 + s**2)
拉普拉斯算子
l = laplacian(f,v) 返回符号域 f 相对于笛卡尔坐标中的向量 v 的拉普拉斯算子. 如果 f 是一个数组, 则该函数计算 f 的每个元素的拉普拉斯算子并返回与 f 大小相同的输出 l.
l = laplacian(f) 返回符号域 f 相对于由 f 中的符号变量构造的默认向量的拉普拉斯算子.
f — 符号域,符号表达式,符号函数,符号矩阵变量,符号矩阵函数
v — 用于查找拉普拉斯算子的向量,符号标量变量的向量,符号函数,符号矩阵变量,符号矩阵函数
示例1：热传导方程 - 二维温度场
T(x,y) = x² + y² 的温度分布
result1 = laplacian("x**2 + y**2")
print(f" T(x,y) = x² + y² 的拉普拉斯: {result1}")
#T(x,y) = x² + y² 的拉普拉斯: 4
#物理意义: 均匀热源分布
示例2：静电场电势 - 点电荷电势
二维点电荷电势: φ(x,y) = ln(√(x²+y²))
result2 = laplacian(log(sqrt(x**2 + y**2)))
print(f" φ(x,y) = ln(r) 的拉普拉斯: {result2}")
#φ(x,y) = ln(r) 的拉普拉斯: (-2*x**2/(x**2 + y**2) + 1)/(x**2 + y**2) + (-2*y**2/(x**2 + y**2) + 1)/(x**2 + y**2)
#物理意义: 二维点电荷产生的电势场
示例3：量子力学波函数 - 谐振子基态
谐振子基态波函数: ψ(x) = exp(-x²/2)
result3 = laplacian(exp(-x**2/2))
print(f" ψ(x) = exp(-x²/2) 的拉普拉斯: {result3}")
#ψ(x) = exp(-x²/2) 的拉普拉斯: (x**2 - 1)*exp(-x**2/2)
#物理意义: 量子谐振子的动能项
示例4：流体力学速度势 - 源流场
速度势函数: φ(x,y) = x/(x²+y²)
result4 = laplacian(x/(x**2 + y**2))
print(f" φ(x,y) = x/(x²+y²) 的拉普拉斯: {result4}")
#φ(x,y) = x/(x²+y²) 的拉普拉斯: 2*x*(4*x**2/(x**2 + y**2) - 3)/(x**2 + y**2)**2 + 2*x*(4*y**2/(x**2 + y**2) - 1)/(x**2 + y**2)**2
#物理意义: 无旋不可压缩流体的速度势
示例5：矩阵场的拉普拉斯 - 向量值函数
result5 = laplacian([[x**2, x*y], [y**2, x**2 + y**2]], [x, y])
print(f" 矩阵函数 F(x,y) 的拉普拉斯:\n{result5}")
#矩阵函数 F(x,y) 的拉普拉斯:
#[[2, 0],
[2, 4]]
#物理意义: 向量场的拉普拉斯运算
示例6：三维波动方程 - 球面波
球面波: u(x,y,z) = sin(√(x²+y²+z²))
result6 = laplacian(sin(sqrt(x**2 + y**2 + z**2)))
print(f" u(r) = sin(r) 的拉普拉斯: {result6}")
#物理意义: 球对称波的传播
示例7：泊松方程验证
如果 ∇²φ = ρ，这里验证一个特解
phi = laplacian(x**4 + y**4)
print(f" φ = x⁴ + y⁴ 的拉普拉斯: {phi}")
#φ = x⁴ + y⁴ 的拉普拉斯: 12*x**2 + 12*y**2
#对应的源项密度: ρ = 12*x**2 + 12*y**2
示例8：调和函数验证
调和函数的拉普拉斯为0
harmonic_func = laplacian(x**3 - 3*x*y**2) # 实部解析函数
print(f" f = x³ - 3xy² (调和函数) 的拉普拉斯: {harmonic_func}")
#f = x³ - 3xy² (调和函数) 的拉普拉斯: 0
#物理意义: 无源场的验证
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def laplacian_equation(input_str):
"""
对标 MATLAB 的 laplacian 函数，计算拉普拉斯算子
参数:
input_str: 输入字符串，支持以下形式：
1. 标量表达式 (自动检测变量)
2. 元组 (表达式, 变量列表)
3. 矩阵 (逐元素计算)
返回:
拉普拉斯算子计算结果 (标量、矩阵或符号表达式)，或错误信息
示例:
>>> laplacian_equation("x**2 + y**2")
4
>>> laplacian_equation("([[x*y, x**2], [y**2, x+y]], [x, y])")
Matrix([
[0, 2],
[2, 0]])
"""
try:
expr = sp.sympify(input_str)
error_msg = None
result = None
def compute_laplacian(f, vars=None):
"""计算拉普拉斯算子的核心函数"""
# 自动检测变量（按字母顺序排序）
if vars is None:
vars = sorted(f.free_symbols, key=lambda s: s.name)
if not vars:
raise ValueError("表达式不含自由变量")
# 验证变量类型
if not all(isinstance(v, sp.Symbol) for v in vars):
raise TypeError("变量列表中包含非符号对象")
# 计算拉普拉斯算子：Σ ∂²f/∂v²
return sum(sp.diff(f, v, v) for v in vars)
# 处理元组输入 (表达式, 变量列表)
if isinstance(expr, tuple):
# 参数数量验证
if len(expr) != 2:
return f"错误：需要 2 个参数的元组 (表达式, 变量列表)，收到 {len(expr)} 个参数"
f_part, vars_part = expr[0], expr[1]
# 变量列表验证
if not isinstance(vars_part, (list, tuple)):
return "错误：第二个参数必须是变量列表"
# 转换为符号变量
try:
vars_list = [sp.sympify(v) for v in vars_part]
except Exception as e:
return f"变量列表解析错误: {str(e)}"
result = compute_laplacian(f_part, vars_list)
# 处理标量输入 (自动检测变量)
elif expr.free_symbols:
result = compute_laplacian(expr)
else:
error_msg = f"无法识别的输入类型: {type(expr)}"
return result if not error_msg else error_msg
except Exception as e:
return f"错误：{str(e)}"
# 测试用例
if __name__ == "__main__":
# 测试1: 标量函数
print("测试1:", laplacian_equation("f(x,y,z)"))
# Derivative(f(x, y, z), (x, 2)) + Derivative(f(x, y, z), (y, 2)) + Derivative(f(x, y, z), (z, 2))
print("测试2:", laplacian_equation("x**2 + y**2 + z**2"))
# 6
# 测试3: 带变量列表
print("测试3:", laplacian_equation("(sin(x)*cos(y), [x, y])"))
# -2*sin(x)*cos(y)
最小公倍数
L = lcm(A,B) 返回 A 和 B 元素的最小公倍数.
A,B — 输入值,标量,向量,实数数组或正整数值
L — 最小公倍数,实数,正整数值
示例1：分数运算中的通分
result1 = lcm(3, 4)
print(f" 分数 1/3 和 1/4 的公分母: {result1}")
print(f" 计算过程: 1/3 + 1/4 = {4}/{result1} + {3}/{result1} = 7/{result1}\n")
#分数 1/3 和 1/4 的公分母: 12
#计算过程: 1/3 + 1/4 = 4/12 + 3/12 = 7/12
示例2：周期性事件的同步
两个事件分别每6天和8天发生一次，求同时发生的周期
result2 = lcm(6, 8)
print(f" 事件A每6天发生，事件B每8天发生")
print(f" 两个事件同时发生的周期: 每{result2}天\n")
#事件A每6天发生，事件B每8天发生
#两个事件同时发生的周期: 每24天
示例3：齿轮传动系统
两个齿轮的齿数，求它们再次对齐的旋转次数
gear1, gear2 = 15, 20
result3 = lcm(gear1, gear2)
cycles1 = result3 // gear1
cycles2 = result3 // gear2
print(f" 齿轮A有{gear1}齿，齿轮B有{gear2}齿")
print(f" 两个齿轮再次对齐需要: 齿轮A转{cycles1}圈，齿轮B转{cycles2}圈\n")
#齿轮A有15齿，齿轮B有20齿
#两个齿轮再次对齐需要: 齿轮A转4圈，齿轮B转3圈
示例4：矩阵运算 - 图像处理
图像处理中的像素操作:
matrix_expr = ([[2, 3, 4], [5, 6, 7]], [[3, 4, 5], [6, 7, 8]])
result4 = lcm(matrix_expr)
print(f" 两个像素矩阵的LCM:\n{result4}")
#[[6,12,20]
[30,42,56]]
#应用: 图像混合时的周期对齐
示例5：音乐节奏的同步
beat1, beat2 = 4, 6 # 两个不同的节奏周期
result5 = lcm(beat1, beat2)
print(f" 节奏A每{beat1}拍循环，节奏B每{beat2}拍循环")
print(f" 两个节奏同时回到起点的拍数: {result5}拍\n")
#节奏A每4拍循环，节奏B每6拍循环
#两个节奏同时回到起点的拍数: 12拍
示例6：工程中的材料切割
两种规格的材料，求最小公倍数来优化切割方案
length1, length2 = 180, 240 # 厘米
result6 = lcm(length1, length2)
print(f" 材料A长{length1}cm，材料B长{length2}cm")
print(f" 最小公共长度单位: {result6}cm")
#材料A长180cm，材料B长240cm
#最小公共长度单位: 720cm
#可同时整除两种长度的最大公约切割方案
示例7：计算机科学中的调度算法
process1, process2 = 12, 18 # 两个进程的执行周期(毫秒)
result7 = lcm(process1, process2)
print(f" 进程A每{process1}ms执行，进程B每{process2}ms执行")
print(f" 两个进程同时调度的周期: 每{result7}ms\n")
#进程A每12ms执行，进程B每18ms执行
#两个进程同时调度的周期: 每36ms
示例8：化学反应的周期同步
reaction1, reaction2 = 9, 15 # 两个反应的周期(分钟)
result8 = lcm(reaction1, reaction2)
print(f" 反应A每{reaction1}分钟完成，反应B每{reaction2}分钟完成")
print(f" 两个反应同时完成的周期: 每{result8}分钟\n")
#反应A每9分钟完成，反应B每15分钟完成
#两个反应同时完成的周期: 每45分钟
示例9：三维数组运算
array1 = [[[2, 3], [4, 5]], [[6, 7], [8, 9]]]
array2 = [[[3, 4], [5, 6]], [[7, 8], [9, 10]]]
result9 = lcm_number(array1, array2)
print(f" 逐元素LCM结果:\n{result9}\n")
#逐元素LCM结果:
#[[[6,12]
[20,30]]
[[42,56]
[72,90]]]
示例10：符号计算
expr1 = 2*x
expr2 = 3*y
symbolic_lcm = lcm(expr1, expr2)
print(f" 符号表达式 {expr1} 和 {expr2} 的LCM: {symbolic_lcm}")
#符号表达式 2*x 和 3*y 的LCM: 6*x*y
示例11：大数据集的批量处理
large_array1 = [list(range(1, 6)), list(range(6, 11))]
large_array2 = [list(range(2, 7)), list(range(7, 12))]
result11 = lcm_number(large_array1, large_array2)
print(f" 大数据集A: {large_array1}")
print(f" 大数据集B: {large_array2}")
print(f" 批量LCM计算结果形状: {np.array(result11).shape}\n")
#大数据集A: [[1, 2, 3, 4, 5],
[6, 7, 8, 9, 10]]
#大数据集B: [[2, 3, 4, 5, 6],
[7, 8, 9, 10, 11]]
#批量LCM计算结果形状: (2, 5)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def lcm_number(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为包含两个元素的元组
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(int(e) for e in expr)
result = np.lcm(*params)
elif any(e.free_symbols for e in expr):
result = sp.lcm(*expr)
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：数值计算
print(lcm_number("(6, 8)"))
# 24
埃尔米特不定矩阵的分块LDL分解
[L,D,P] = ldl(A)将满矩阵A分解为置换的下三角矩阵L和满足A = L*D*L'的分块对角矩阵D, 以向量形式返回置换信息P.
A — 输入矩阵,实数或复数方阵
L — 下三角因子,方阵
D — 分块对角因子,方阵
P — 置换信息,矩阵,向量
示例1：结构力学刚度矩阵
stiffness_matrix = "[[2, -1, 0], [-1, 2, -1], [0, -1, 2]]"
result = ldl(stiffness_matrix)
L, D, perm = result
print("L 矩阵:")
print(L)
#L 矩阵:
#[[1.0, 0, 0],
[-0.5, 1.0, 0],
[0, -0.666666666666667, 1.0]]
print("D 矩阵:")
print(D)
#D 矩阵:
#[[2.0, 0, 0],
[0, 1.5, 0],
[0, 0, 1.33333333333333]]
print("置换矩阵:")
print(perm)
#置换矩阵:
#[[0],
[1],
[2]]
验证分解正确性：L * D * L^T 应该等于原矩阵
reconstructed = L * D * L.T
print("重建的矩阵:")
print(reconstructed)
#重建的矩阵:
#[[2.0, -1.0, 0],
[-1.0, 2.0, -1.0],
[0, -1.0, 2.0]]
示例2：金融协方差矩阵
covariance_matrix = [[0.04, 0.02, 0.01], [0.02, 0.09, 0.03], [0.01, 0.03, 0.16]]
result = ldl(covariance_matrix)
L, D, perm = result
print("L 矩阵:")
print(L)
#L 矩阵:
#[[1.0, 0, 0],
[0.5, 1.0, 0],
[0.25, 0.3125, 1.0]]
print("D 矩阵:")
print(D)
#D 矩阵:
#[[0.04, 0, 0],
[0, 0.08, 0],
[0, 0, 0.1496875]]
在风险管理中，D矩阵的对角线元素代表各因子的方差
print("因子方差:", [D[i, i] for i in range(D.shape[0])])
#因子方差: [0.04, 0.08, 0.1496875]
示例3：电路导纳矩阵（符号）
circuit_matrix = [[1/R1 + 1/R2, -1/R2], [-1/R2, 1/R2 + 1/R3]]
使用具体的电阻值
numeric_circuit = [[1/100 + 1/200, -1/200], [-1/200, 1/200 + 1/300]]
result = ldl(numeric_circuit)
L, D, perm = result
print("L 矩阵:")
print(L)
#L 矩阵:
#[[1.0, 0],
[-0.333333333333333, 1.0]]
print("D 矩阵:")
print(D)
#D 矩阵:
#[[0.015, 0],
[0, 0.00666666666666667]]
示例4：最小二乘正规方程
A_matrix = [[1, 1], [1, 2], [1, 3]]
A = eval(A_matrix)
ATA = A.T * A
print(f"正规方程矩阵 A^T A: {ATA}")
#正规方程矩阵 A^T A:
#[[3, 6],
[6, 14]]
示例5：优化问题Hessian矩阵
函数 f(x,y) = ax² + bxy + cy² 的Hessian矩阵
hessian_symbolic = [[2*a, b], [b, 2*c]]
result = ldl(hessian_symbolic)
L, D = result
print("L 矩阵:")
print(L)
#L 矩阵:
#[[1, 0],
[b/(2*a), 1]]
print("D 矩阵:")
print(D)
#D 矩阵:
#[[2*a, 0],
[0, 2*c - b**2/(2*a)]]
检查正定性：如果D的所有对角元素为正，则矩阵正定
print("D的对角元素:", [D[i, i] for i in range(D.shape[0])])
#D的对角元素: [2*a, 2*c - b**2/(2*a)]
示例6：图像结构张量
structure_tensor = [[125, 45], [45, 85]]
result = ldl(structure_tensor)
L, D, perm = result
print("L 矩阵:")
print(L)
#L 矩阵:
#[[1.0, 0],
[0.36, 1.0]]
print("D 矩阵:")
print(D)
#D 矩阵:
#[[125.0, 0],
[0, 68.8]]
特征值（D的对角线）可以用于角点检测
eigenvalues = [D[i, i] for i in range(D.shape[0])]
print("特征值（用于角点检测）:", eigenvalues)
#特征值（用于角点检测）: [125.0, 68.8]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.linalg import lu, ldl
def ldl_decomposition_matrix(input_str):
"""
执行 LDL 分解的主函数
参数:
input_str: 输入矩阵的字符串表达式(如 "[[1, 2], [3, 4]]")
返回:
- 符号矩阵: 返回 (L, D) 元组
- 数值矩阵: 返回 (L, D, perm) 元组
- 错误输入: 返回错误消息字符串
注意:
- 当矩阵包含符号时自动声明符号为实数类型
- 自动检查矩阵对称性（符号矩阵场景）
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is not None:
# 检测符号和虚数存在性
contains_symbols = A.is_symbolic()
has_imaginary = any(element.has(sp.I) for element in A)
if contains_symbols:
# 自动声明所有符号为实数
symbols = A.free_symbols
real_subs = {s: sp.Symbol(s.name, real=True) for s in symbols}
A = A.subs(real_subs)
# 显式检查矩阵对称性
if A != A.T:
return f"错误: 矩阵必须对称，当前矩阵不对称\n{A}"
# 执行 LDL 分解
try:
l, d = A.LDLdecomposition()
result = (l, d)
except ValueError as e:
return f"分解错误: {e}"
elif has_imaginary:
# 执行 LDL 分解
try:
l, d = A.LDLdecomposition()
result = (l, d)
except ValueError as e:
return f"分解错误: {e}"
else:
# 数值矩阵处理
try:
N_A = np.array(A, dtype=float)
lu_mat, d_mat, perm = ldl(N_A)
result = (
sp.Matrix(lu_mat),
sp.Matrix(d_mat),
sp.Matrix(perm)
)
except np.linalg.LinAlgError as e:
return f"数值分解错误: {e}"
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: 符号矩阵分解")
symbolic_matrix = "[[a, b], [b, c]]"
print(f"输入矩阵: {symbolic_matrix}")
L, D = ldl_decomposition_matrix(symbolic_matrix)
print("L 矩阵:")
print(L)
# Matrix([[1, 0],
# [b/a, 1]])
print("D 矩阵:")
print(D)
# Matrix([[a, 0],
# [0, c - b**2/a]])
print("\n" + "=" * 50 + "\n")
# 示例 2: 实数数值矩阵分解
print("示例 2: 数值矩阵分解")
numeric_matrix = "[[4, 12, -16], [12, 37, -43], [-16, -43, 98]]"
print(f"输入矩阵: {numeric_matrix}")
L, D, perm = ldl_decomposition_matrix(numeric_matrix)
print("L 矩阵:")
print(L)
# Matrix([[1.00000000000000, 0, 0],
# [3.00000000000000, 0.147058823529412, 1.00000000000000],
# [-4.00000000000000, 1.00000000000000, 0]])
print("D 矩阵:")
print(D)
# Matrix([[4.00000000000000, 0, 0],
# [0, 34.0000000000000, 0],
# [0, 0, 0.264705882352941]])
print("置换矩阵:")
print(perm)
# Matrix([[0],
# [2],
# [1]])
print("\n" + "=" * 50 + "\n")
勒奇超越函数
LenchPhi(z,s,a)对于某些特殊参数，LerchPhi会自动计算为精确值.
LerchPhi可以评估为任意的数值精度.
z,s,a -- 输入, 数字,符号变量,向量,矩阵
示例1：数论应用 - Riemann Zeta 和 Hurwitz Zeta 函数
LerchPhi 是更一般的函数，包含其他zeta函数作为特例
Riemann Zeta: ζ(s) = Φ(1, s, 1)
print("黎曼Zeta函数 ζ(2) =", LerchPhi(1, 2, 1))
#黎曼Zeta函数 ζ(2) = 1.64493406684823
Hurwitz Zeta: ζ(s,a) = Φ(1, s, a)
print("Hurwitz Zeta ζ(2, 0.5) =", LerchPhi(1, 2, 0.5))
#Hurwitz Zeta ζ(2, 0.5) = 4.93480220054468
Polylogarithm: Li_s(z) = zΦ(z, s, 1)
print("Polylogarithm Li_2(0.5) =", 0.5 * LerchPhi(0.5, 2, 1))
#Polylogarithm Li_2(0.5) = 0.582240526465013
示例2：统计物理 - 玻色-爱因斯坦分布
在玻色-爱因斯坦统计中，LerchPhi用于计算配分函数
对于理想玻色气体，某些量可以表示为LerchPhi函数
化学势相关计算
temperature = 0.1 # 归一化温度
chemical_potential = 0.05
粒子数密度计算（简化模型）
result = LerchPhi(exp(chemical_potential), 1.5, 1)
print(f"玻色气体粒子数密度（近似）: {result}")
#玻色气体粒子数密度（近似）: 2.41526400326651 - 0.754006708881948*I
不同温度下的行为
temperatures = [0.1, 0.5, 1.0, 2.0]
for T in temperatures:
z = np.exp(1 / T - 1) # 简化模型中的逸度
density = LerchPhi(z, 1.5, 1)
print(f"温度 {T}: 密度 ≈ {density}")
#温度 0.1: 密度 ≈ -0.00242966705008955 - 0.00131242909495758*I
#温度 0.5: 密度 ≈ 0.384146086327461 - 1.30409866434658*I
#温度 1.0: 密度 ≈ 2.61237534868549
#温度 2.0: 密度 ≈ 1.33627284845045
示例3：量子力学 - 库仑势能相关
在量子力学中，LerchPhi出现在库仑势的波函数计算中
特别是氢原子波函数的傅里叶变换
量子力学中的特殊值:
print("Φ(0.5, 1, 1) =", LerchPhi(0.5, 1, 1))
print("Φ(0.5, 2, 1) =", LerchPhi(0.5, 2, 1))
print("Φ(0.5, 2, 0.5) =", LerchPhi(0.5, 2, 0.5))
#Φ(0.5, 1, 1) = 1.38629436111989
#Φ(0.5, 2, 1) = 1.16448105293003
#Φ(0.5, 2, 0.5) = 4.27714550058095
参数扫描 (固定 a=1):
for s in [0.5, 1.0, 1.5, 2.0]:
for z in [0.1, 0.5, 0.9]:
value = LerchPhi(z, s, 1)
print(f"Φ({z}, {s}, 1) = {value:.6f}")
#Φ(0.1, 0.5, 1) = 1.077033
#Φ(0.5, 0.5, 1) = 1.612253
#Φ(0.9, 0.5, 1) = 4.468834
#Φ(0.1, 1.0, 1) = 1.053605
#Φ(0.5, 1.0, 1) = 1.386294
#Φ(0.9, 1.0, 1) = 2.558428
#Φ(0.1, 1.5, 1) = 1.037415
#Φ(0.5, 1.5, 1) = 1.249674
#Φ(0.9, 1.5, 1) = 1.793821
#Φ(0.1, 2.0, 1) = 1.026178
#Φ(0.5, 2.0, 1) = 1.164481
#Φ(0.9, 2.0, 1) = 1.444127
示例4：矩阵输入计算
创建参数矩阵
z_matrix = [[0.1, 0.2], [0.3, 0.4]]
s_value = 1.5
a_value = 1.0
矩阵形式的LerchPhi计算
matrix_input = (z_matrix, s_value, a_value)
result = LerchPhi(matrix_input)
print("LerchPhi 结果矩阵:")
print(result)
#LerchPhi 结果矩阵:
#[[1.03741452346169, 1.07957769907279],
[1.12770365181602, 1.18353042309162]]
示例5：级数求和分析
LerchPhi可以看作是一个广义的级数求和
def verify_lerch_series(z, s, a, terms=100):
"""通过级数部分和验证LerchPhi结果"""
partial_sum = 0
for n in range(terms):
term = (z ** n) / ((a + n) ** s)
partial_sum += term
return partial_sum
比较解析结果和级数部分和
z, s, a = 0.5, 1.2, 0.8
lerch_result = LerchPhi(z, s, a)
series_result = verify_lerch_series(float(z), float(s), float(a))
print(f"LerchPhi解析结果: {lerch_result}")
print(f"级数部分和(100项): {series_result}")
print(f"相对误差: {abs(float(lerch_result) - series_result) / abs(float(lerch_result)):.2e}")
#LerchPhi解析结果: 1.66792644603551
#级数部分和(100项): 1.6679264460355134
#相对误差: 1.33e-16
示例6：特殊函数关系验证
验证LerchPhi与其他特殊函数的关系
与digamma函数的关系: Φ(1,1,a) = -ψ(a)
print("Φ(1, 1, 0.5) =", LerchPhi(1, 1, 0.5))
#Φ(1, 1, 0.5) = oo
与指数积分的关系
z_val = 0.3
result = LerchPhi((z_val, 1, 1)
print(f"Φ({z_val}, 1, 1) = {result}")
#Φ(0.3, 1, 1) = 1.18891647979577
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from mpmath import lerchphi
def lench_phi_transcendent(input_str):
"""
对标 MATLAB 的 LerchPhi 函数的实现，支持标量和矩阵输入
Lerch Phi 函数定义:
Φ(z, s, a) = ∑_{n=0}^∞ z^n / (a + n)^s
分类归属:
属于特殊函数中的超越函数（Transcendental Function），
是广义的 Zeta 函数，包含 Riemann Zeta 和 Hurwitz Zeta 函数作为特例
参数:
input_str: 输入表达式字符串，格式为元组 (n, a, x) 或单个数值，
支持矩阵输入（需确保维度一致）
返回:
计算结果（标量/矩阵）或错误信息字符串
示例:
>>> lench_phi_transcendent("(0.5, 1, 0.2)")
1.64493406684823
>>> lench_phi_transcendent("([[0.5, 1], [2, 3]], 1, 0.2)")
Matrix([
[1.64493406684823, 2.16439397443791],
[3.30685281944005, 4.17469149400356]])
"""
try:
# 解析输入为 SymPy 表达式
expr = sp.sympify(input_str, evaluate=False)
error = False
# 参数校验
if not isinstance(expr, tuple) or len(expr) != 3:
return f"输入格式错误: 需要三元组 (n, a, x)，实际输入 {input_str}"
if all(e.is_number for e in expr):
params = tuple(complex(e.evalf()) for e in expr)
result = lerchphi(*params)
elif any(e.free_symbols for e in expr):
result = sp.lerchphi(*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(lench_phi_transcendent("(0.5, 1, 0.2)"))
# (5.59472668491108 + 0.0j)
print(lench_phi_transcendent("(2, 1, 1.5)"))
# (-0.376774759859769 - 1.11072073453959j)
第一类勒让德多项式
P = legendre(n,X) 为 X 中的每个元素计算阶数为n,级数为m = 0, 1, ..., n 时的第一类勒让德多项式.
n — 勒让德函数的阶,正整数
X — 输入值,标量,向量,矩阵
示例1：电磁学多极展开
在电磁学中，电势的多极展开使用Legendre多项式
def electric_potential_multipole(r, theta, charges):
"""
计算点电荷系统的电势（多极展开）
r: 观察点距离
theta: 观察点角度
charges: [(q_i, d_i)] 电荷量和位置列表
"""
potential = 0
for n in range(3): # 计算到四极矩
moment = 0
for q, d in charges:
moment += q * (d ** n)
P_n = legendreP(n, cos(theta))
term = moment * P_n / (r ** (n + 1))
potential += term
print(f"{n}极矩项: {term:.6f}")
return potential
#0极矩项: 0.000000
#1极矩项: 0.176777
#2极矩项: 0.000000
两个点电荷系统
charges = [(1.0, 0.5), (-1.0, -0.5)] # (电荷量, 位置)
r = 2.0
theta = sp.pi / 4
result = electric_potential_multipole(r, theta, charges)
print(f"总电势: {result:.6f}")
print("\n" + "=" * 50 + "\n")
#总电势: 0.176777
示例2：量子力学角动量波函数
def spherical_harmonic_theta_part(l, m, theta):
"""计算球谐函数的θ部分（关联Legendre多项式）"""
# 对于m=0的情况，就是普通的Legendre多项式
if m == 0:
return legendreP(l, cos(theta))
else:
# 对于m≠0，使用关联Legendre多项式
# 首先创建符号变量
x = sp.Symbol('x')
# 计算普通Legendre多项式的符号表达式
P_l = sp.legendre(l, x)
# 计算m阶导数
P_lm = (-1) ** m * (1 - x ** 2) ** (m / 2) * sp.diff(P_l, x, m)
# 代入x = cos(theta)
return P_lm.subs(x, sp.cos(theta))
计算几个量子态的角向波函数
quantum_states = [(0, 0), (1, 0), (1, 1), (2, 0)]
theta_vals = [0, sp.pi / 4, sp.pi / 2, 3 * sp.pi / 4, sp.pi]
球谐函数角向部分:
for l, m in quantum_states:
print(f"\nY_{l}^{m} 的θ部分:")
for theta_val in theta_vals:
value = spherical_harmonic_theta_part(l, m, theta_val)
# 如果是符号表达式，尝试数值化
if isinstance(value, sp.Expr):
try:
value_num = float(value)
print(f" θ={theta_val:.2f}: {value_num:.4f}")
except:
print(f" θ={theta_val:.2f}: {value}")
else:
print(f" θ={theta_val:.2f}: {float(value):.4f}")
#Y_0^0 的θ部分:
# θ=0.00: 1.0000
# θ=0.79: 1.0000
# θ=1.57: 1.0000
# θ=2.36: 1.0000
# θ=3.14: 1.0000
#Y_1^0 的θ部分:
# θ=0.00: 1.0000
# θ=0.79: 0.7071
# θ=1.57: 0.0000
# θ=2.36: -0.7071
# θ=3.14: -1.0000
#Y_1^1 的θ部分:
# θ=0.00: 0.0000
# θ=0.79: -0.7071
# θ=1.57: -1.0000
# θ=2.36: -0.7071
# θ=3.14: 0.0000
#Y_2^0 的θ部分:
# θ=0.00: 1.0000
# θ=0.79: 0.2500
# θ=1.57: -0.5000
# θ=2.36: 0.2500
# θ=3.14: 1.0000
示例3：高斯-勒让德求积法
def gauss_legendre_quadrature(n, func, a=-1, b=1):
"""
使用n点高斯-勒让德求积法计算积分
"""
# 获取求积节点和权重（这里简化为使用标准值）
# 实际应用中应该使用预先计算的节点和权重表
# 对于演示，我们计算几个低阶多项式的积分
if n == 2:
nodes = [-1 / sp.sqrt(3), 1 / sp.sqrt(3)]
weights = [1, 1]
elif n == 3:
nodes = [-sp.sqrt(3 / 5), 0, sp.sqrt(3 / 5)]
weights = [5 / 9, 8 / 9, 5 / 9]
else:
nodes = [0] # 简化处理
weights = [2]
# 变换积分区间
integral = 0
for i, node in enumerate(nodes):
x_transformed = (b - a) / 2 * node + (a + b) / 2
integral += weights[i] * func(x_transformed)
integral *= (b - a) / 2
return integral
测试积分 ∫_{-1}^{1} x² dx = 2/3
def test_func(x):
return x ** 2
for n_points in [2, 3]:
result = gauss_legendre_quadrature(n_points, test_func)
exact = 2 / 3
error = abs(result - exact)
print(f"{n_points}点高斯求积: {result}, 误差: {error}")
#2点高斯求积: 0.666666666666667, 误差: 0
#3点高斯求积: 0.666666666666667, 误差: 1.11022302462516E-16
示例4：信号处理中的正交基函数
def legendre_basis_expansion(signal, degree):
"""
使用Legendre多项式基展开信号
"""
coefficients = []
for n in range(degree + 1):
# 计算第n个基函数的系数（简化版本）
def basis_func(x):
return legendreP(n, x)
# 数值积分计算内积（实际应用中使用离散求和）
coeff = gauss_legendre_quadrature(5, lambda x: signal(x) * basis_func(x))
coefficients.append(coeff)
print(f"系数 a_{n}: {coeff:.6f}")
return coefficients
测试信号：f(x) = x³
def test_signal(x):
return x ** 3
print("信号 f(x) = x³ 的Legendre展开系数:")
coeffs = legendre_basis_expansion(test_signal, 3)
#信号 f(x) = x³ 的Legendre展开系数:
#系数 a_0: 0.000000
#系数 a_1: 0.000000
#系数 a_2: 0.000000
#系数 a_3: 0.000000
重建信号
def reconstructed_signal(x, coefficients):
result = 0
for n, coeff in enumerate(coefficients):
result += coeff * legendreP(n, x)
return result
验证重建精度
test_points = [-0.8, -0.4, 0, 0.4, 0.8]
print("\n重建信号验证:")
for x in test_points:
original = test_signal(x)
reconstructed = reconstructed_signal(x, coeffs)
error = abs(original - reconstructed)
print(f"x={x:5.1f}: 原值={original:7.4f}, 重建={reconstructed:7.4f}, 误差={error:.2e}")
#x= -0.8: 原值=-0.5120, 重建= 0.0000, 误差=5.12e-01
#x= -0.4: 原值=-0.0640, 重建= 0.0000, 误差=6.40e-02
#x= 0.0: 原值= 0.0000, 重建= 0.0000, 误差=0.00e+00
#x= 0.4: 原值= 0.0640, 重建= 0.0000, 误差=6.40e-02
#x= 0.8: 原值= 0.5120, 重建= 0.0000, 误差=5.12e-01
示例5：球坐标下的热传导
def spherical_heat_solution(r, theta, t, initial_temp=100):
"""
球对称热传导问题的解析解（使用Legendre多项式）
"""
# 简化模型：只考虑前几项展开
solution = 0
for n in range(3):
# 特征值 λ_n = n(n+1) 对于球对称问题
lambda_n = n * (n + 1)
# 时间衰减因子
time_factor = sp.exp(-lambda_n * t)
# 空间部分（Legendre多项式）
spatial_part = legendreP(n, cos(theta))
# 径向部分（球贝塞尔函数，这里简化）
radial_part = 1 / r if n == 0 else r ** n
term = time_factor * spatial_part * radial_part
solution += term
print(f"n={n}项: {float(term):.6f}")
return initial_temp * solution
计算不同位置和时间的温度分布
positions = [(1.0, 0), (1.0, sp.pi / 2), (2.0, 0)]
times = [0.1, 0.5, 1.0]
球坐标热传导温度分布:
for r, theta in positions:
print(f"\n位置 (r={r}, θ={theta:.2f}):")
for t in times:
temp = spherical_heat_solution(r, theta, t)
print(f" 时间 t={t}: 温度={float(temp):.2f}°C")
#位置 (r=1.0, θ=0.00):
#n=0项: 1.000000
#n=1项: 0.818731
#n=2项: 0.548812
# 时间 t=0.1: 温度=236.75°C
#n=0项: 1.000000
#n=1项: 0.367879
#n=2项: 0.049787
# 时间 t=0.5: 温度=141.77°C
#n=0项: 1.000000
#n=1项: 0.135335
#n=2项: 0.002479
# 时间 t=1.0: 温度=113.78°C
#位置 (r=1.0, θ=1.57):
#n=0项: 1.000000
#n=1项: 0.000000
#n=2项: -0.274406
# 时间 t=0.1: 温度=72.56°C
#n=0项: 1.000000
#n=1项: 0.000000
#n=2项: -0.024894
# 时间 t=0.5: 温度=97.51°C
#n=0项: 1.000000
#n=1项: 0.000000
#n=2项: -0.001239
# 时间 t=1.0: 温度=99.88°C
#位置 (r=2.0, θ=0.00):
#n=0项: 0.500000
#n=1项: 1.637462
#n=2项: 2.195247
# 时间 t=0.1: 温度=433.27°C
#n=0项: 0.500000
#n=1项: 0.735759
#n=2项: 0.199148
# 时间 t=0.5: 温度=143.49°C
#n=0项: 0.500000
#n=1项: 0.270671
#n=2项: 0.009915
# 时间 t=1.0: 温度=78.06°C
示例6：地球重力场建模
def earth_gravity_potential(r, theta, GM=3.986e14, R=6371e3):
"""
地球重力势的球谐展开（简化版本）
参数:
r: 距离地心的距离 (m)
theta: 余纬角 (弧度)
GM: 地球引力常数 (m³/s²)
R: 地球参考半径 (m)
"""
# 主要项：均匀球体势
potential = GM / r
# 添加非球形修正（使用Legendre多项式）
# J2项：地球扁率影响
J2 = 1.08263e-3
P2 = legendreP(2, cos(theta))
potential += (GM * R ** 2 / r ** 3) * J2 * P2
print(f"均匀球体势: {GM / r:.6e}")
print(f"J2修正项: {(GM * R ** 2 / r ** 3) * J2 * float(P2):.6e}")
return potential
计算不同纬度的重力势
GM = 3.986e14 # 地球引力常数 (m³/s²)
R = 6371e3 # 地球参考半径 (m)
latitudes = [0, 30, 60, 90] # 纬度（度）
height = 200e3 # 高度200km
地球重力势分布（高度200km）:
for lat in latitudes:
theta = (90 - lat) * sp.pi / 180 # 转换为余纬
r = R + height
potential = earth_gravity_potential(r, theta, GM, R)
print(f"纬度 {lat}°: 重力势 = {potential:.6e} m²/s²")
#均匀球体势: 6.066048e+07
#J2修正项: -3.086798e+04
#纬度 0°: 重力势 = 6.062961e+07 m²/s²
#均匀球体势: 6.066048e+07
#J2修正项: -7.716994e+03
#纬度 30°: 重力势 = 6.065276e+07 m²/s²
#均匀球体势: 6.066048e+07
#J2修正项: 3.858497e+04
#纬度 60°: 重力势 = 6.069906e+07 m²/s²
#均匀球体势: 6.066048e+07
#J2修正项: 6.173595e+04
#纬度 90°: 重力势 = 6.072221e+07 m²/s²
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.special import eval_legendre
def legendre_polynomial(input_str):
"""
对标 MATLAB 的 legendreP 函数，计算 Legendre 多项式
参数:
input_str: 输入字符串，可以是以下形式：
1. 元组 (n, x)：标量或矩阵参数
2. 单个数值/符号：默认视为 n=输入值，x=自动符号化
返回:
Legendre 多项式计算结果 (数值、矩阵或符号表达式)，或错误信息
示例:
>>> legendre_polynomial("(0, 0.5)")
1.00000000000000
>>> legendre_polynomial("([[1,2], [3,4]], [[0.1,0.2], [0.3,0.4]])")
Matrix([
[0.100000000000000, 0.200000000000000],
[ -0.365000000000000, -1.12000000000000]])
"""
try:
expr = sp.sympify(input_str)
error_msg = None
result = None
def legendre_p_sym(n, x):
"""核心计算函数，包含参数验证"""
# 数值型n的验证
if n.is_number:
# 检查是否为整数
if not n.is_integer:
raise ValueError("阶数n必须为整数")
# 转换为Python整数
n_int = int(n)
# 检查非负性
if n_int < 0:
raise ValueError("阶数n必须为非负整数")
# 计算Legendre多项式
return sp.legendre(n, x)
def legendre_p_sic(n, x):
"""
对标 MATLAB 的 legendreP 函数，计算Legendre多项式值（仅标量数值参数）
参数:
n (int): 多项式阶数（非负整数）
x (float): 自变量值
返回:
float: Legendre多项式在x处的值
异常:
ValueError: 如果n不是整数或为负数
"""
# 验证n是否为非负整数
if not isinstance(n, (int, np.integer)):
raise ValueError("阶数n必须是整数")
if n < 0:
raise ValueError("阶数n必须为非负整数")
# 计算Legendre多项式值
return eval_legendre(n, x)
# 处理元组输入 (n, x)
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
n = int(expr[0])
x = float(expr[1])
result = legendre_p_sic(n, x)
elif any(e.free_symbols for e in expr):
result = legendre_p_sym(*expr)
else:
error = True
else:
# 如果输入是纯数字，视为n=输入值，x=符号x
result = legendre_p_sym(expr, sp.symbols('x'))
return result
except Exception as e:
return f"错误：{str(e)}"
# 测试用例
if __name__ == "__main__":
# 测试1: 标量计算
print("测试1:", legendre_polynomial("(0, 0.5)"))
# 1.0
print("测试2:", legendre_polynomial("(2, 0.5)"))
# -0.125
print("测试4:", legendre_polynomial("(2, x)"))
# 3*x**2/2 - 1/2
print("测试5:", legendre_polynomial("3"))
# 5*x**3/2 - 3*x/2
print("测试6:", legendre_polynomial("(2.5, 0.1)"))
# -0.485
第二类勒让德多项式
P = legendreq(n,X) 为 X 中的每个元素计算阶数为n,级数为m = 0, 1, ..., n 时的第二类勒让德多项式.
n — 勒让德函数的阶,正整数
X — 输入值,标量,向量,矩阵
示例1：球坐标拉普拉斯方程的外部解
def spherical_laplace_exterior(r, theta, coefficients):
"""
球坐标下拉普拉斯方程的外部解（使用第二类Legendre函数）
适用于 r > R 的区域
"""
potential = 0
for n, A_n in enumerate(coefficients):
# 外部解使用 Q_n(cosθ)/r^(n+1) 形式
# 但实际中更常用 P_n(cosθ)/r^(n+1)，这里展示Q函数的应用
Q_n = legendreq(n, cos(theta))
term = A_n * Q_n / (r ** (n + 1))
potential += term
print(f"n={n}项: {float(term):.6e}")
return potential
计算外部电势分布
coefficients = [1.0, 0.5, 0.2] # 展开系数
r_values = [1.5, 2.0, 3.0]
theta_values = [0, pi / 4, pi / 2]
球坐标拉普拉斯方程外部解:
for r in r_values:
for theta in theta_values:
potential = spherical_laplace_exterior(r, theta, coefficients)
print(f"r={r}, θ={theta / sp.pi:.2f}π: 电势 = {float(potential):.6e}")
print()
#n=0项: 6.666667e+299
#n=1项: 2.222222e+299
#n=2项: 5.925926e+298
#r=1.5, θ=0.00π: 电势 = 9.481481e+299
#n=0项: 5.875824e-01
#n=1项: -8.372772e-02
#n=2项: -4.979655e-02
#r=1.5, θ=0.25π: 电势 = 4.540581e-01
#n=0项: 0.000000e+00
#n=1项: -2.222222e-01
#n=2项: -0.000000e+00
#r=1.5, θ=0.50π: 电势 = -2.222222e-01
#n=0项: 5.000000e+299
#n=1项: 1.250000e+299
#n=2项: 2.500000e+298
#r=2.0, θ=0.00π: 电势 = 6.500000e+299
#n=0项: 4.406868e-01
#n=1项: -4.709684e-02
#n=2项: -2.100792e-02
#r=2.0, θ=0.25π: 电势 = 3.725820e-01
#n=0项: 0.000000e+00
#n=1项: -1.250000e-01
#n=2项: -0.000000e+00
#r=2.0, θ=0.50π: 电势 = -1.250000e-01
#n=0项: 3.333333e+299
#n=1项: 5.555556e+298
#n=2项: 7.407407e+297
#r=3.0, θ=0.00π: 电势 = 3.962963e+299
#n=0项: 2.937912e-01
#n=1项: -2.093193e-02
#n=2项: -6.224569e-03
#r=3.0, θ=0.25π: 电势 = 2.666347e-01
#n=0项: 0.000000e+00
#n=1项: -5.555556e-02
#n=2项: -0.000000e+00
#r=3.0, θ=0.50π: 电势 = -5.555556e-02
示例2：量子力学散射问题
def quantum_scattering_phase_shift(k, l, V0=1.0, R=1.0):
"""
计算球对称势散射中的相移（使用第二类Legendre函数）
简化模型：方势阱散射
"""
# 对于高角动量量子数，第二类Legendre函数出现在散射幅度的计算中
# 这里展示其在渐近解中的应用
kr = k * R # 修正：使用已定义的变量kr
# 球贝塞尔函数与第二类Legendre函数的关系
# 在外部区域，波函数包含第二类球汉克尔函数（与Q函数相关）
if l == 0:
# s波散射 - 使用更准确的公式
# 对于方势阱，s波相移的近似公式
phase_shift = np.arctan(np.tan(kr) / kr - 1 / kr)
else:
# 高角动量散射（近似）
# 避免在奇点x=1处计算，使用接近1的值
x_val = 0.999 # 接近1但避免奇点
Q_l = legendreq(l, x_val)
phase_shift = float(Q_l) / (2 * l + 1) * V0
return phase_shift
计算不同角动量的散射相移
量子散射相移计算:
k_values = [0.5, 1.0, 2.0]
l_values = [0, 1, 2, 3]
for k in k_values:
print(f"\n波数 k = {k}:")
for l in l_values:
delta = quantum_scattering_phase_shift(k, l)
print(f" l={l}: 相移 δ_{l} = {delta:.4f} rad")
#波数 k = 0.5:
# l=0: 相移 δ_0 = -0.7369 rad
# l=1: 相移 δ_1 = 0.9321 rad
# l=2: 相移 δ_2 = 0.4581 rad
# l=3: 相移 δ_3 = 0.2784 rad
#波数 k = 1.0:
# l=0: 相移 δ_0 = 0.5085 rad
# l=1: 相移 δ_1 = 0.9321 rad
# l=2: 相移 δ_2 = 0.4581 rad
# l=3: 相移 δ_3 = 0.2784 rad
#波数 k = 2.0:
# l=0: 相移 δ_0 = -1.0101 rad
# l=1: 相移 δ_1 = 0.9321 rad
# l=2: 相移 δ_2 = 0.4581 rad
# l=3: 相移 δ_3 = 0.2784 rad
示例3：地球外部重力场的完整球谐展开
def earth_external_gravity(r, theta, phi, GM=3.986e14, R=6371e3, max_degree=2):
"""
地球外部重力场的完整球谐展开
包含第一类和第二类Legendre函数
"""
# 主要项
potential = GM / r
# 球谐展开（包含第二类函数项）
# 在实际应用中，第二类函数通常用于外部场展开
J2 = 1.08263e-3 # 地球扁率系数
for n in range(1, max_degree + 1):
if n == 2: # 主要非球形项
# 第一类Legendre多项式
P_n = sp.legendre(n, sp.cos(theta))
# 第二类Legendre函数（用于某些边界条件）
Q_n = legendreq(n, cos(theta))
# 组合贡献（简化模型）
term = (R / r) ** (n + 1) * (P_n + 0.1 * Q_n) * J2
potential += GM / R * term
print(f"n={n}: P_{n}项 = {float((R / r) ** (n + 1) * P_n * J2):.6e}")
print(f"n={n}: Q_{n}项 = {float(0.1 * (R / r) ** (n + 1) * Q_n * J2):.6e}")
return potential
地球外部重力场（包含第二类Legendre函数）:
latitudes = [0, 45, 90]
heights = [100e3, 500e3, 1000e3] # 高度
for lat in latitudes:
theta = (90 - lat) * sp.pi / 180 # 余纬
print(f"\n纬度 {lat}°:")
for h in heights:
r = 6371e3 + h
potential = earth_external_gravity(r, theta, 0)
print(f" 高度 {h / 1000:.0f}km: 重力势 = {potential:.6e} m²/s²")
#纬度 0°:
#n=2: P_2项 = -5.166051e-04
#n=2: Q_2项 = -0.000000e+00
# 高度 100km: 重力势 = 6.156558e+7 m²/s²
#n=2: P_2项 = -4.315320e-04
#n=2: Q_2项 = -0.000000e+00
# 高度 500km: 重力势 = 5.798494e+7 m²/s²
#n=2: P_2项 = -3.495374e-04
#n=2: Q_2项 = -0.000000e+00
# 高度 1000km: 重力势 = 5.405492e+7 m²/s²
#纬度 45°:
#n=2: P_2项 = 2.583025e-04
#n=2: Q_2项 = -8.682238e-05
# 高度 100km: 重力势 = 6.160863e+7 m²/s²
#n=2: P_2项 = 2.157660e-04
#n=2: Q_2项 = -7.252472e-05
# 高度 500km: 重力势 = 5.802090e+7 m²/s²
#n=2: P_2项 = 1.747687e-04
#n=2: Q_2项 = -5.874443e-05
# 高度 1000km: 重力势 = 5.408405e+7 m²/s²
#纬度 90°:
#n=2: P_2项 = 1.033210e-03
#n=2: Q_2项 = 1.033210e+296
# 高度 100km: 重力势 = 6.464253e+303 m²/s²
#n=2: P_2项 = 8.630641e-04
#n=2: Q_2项 = 8.630641e+295
# 高度 500km: 重力势 = 5.399738e+303 m²/s²
#n=2: P_2项 = 6.990748e-04
#n=2: Q_2项 = 6.990748e+295
# 高度 1000km: 重力势 = 4.373744e+303 m²/s²
示例4：球面声波辐射
def spherical_sound_radiation(r, theta, omega, k, a=1.0, p0=1.0):
"""
球面声波辐射问题的外部声场
使用第二类Legendre函数表示角向分布
"""
# 波数
k_val = k
# 声压场（外部解）
pressure = 0
for n in range(3):
# 球汉克尔函数（与第二类Legendre函数相关）
# 对于大r，球汉克尔函数 ~ exp(ikr)/(ikr)
# 角向分布（第二类Legendre函数）
Q_n = legendreq(n, cos(theta))
# 径向传播
radial_part = sp.exp(1j * k_val * r) / (1j * k_val * r)
# 球面辐射系数
A_n = p0 * (2 * n + 1) * 1j ** n / (k_val * a)
term = A_n * radial_part * Q_n
pressure += term
magnitude = abs(complex(term))
print(f"n={n}模式: 幅值 = {magnitude:.6f}")
return pressure
球面声波辐射场分布:
frequencies = [100, 500, 1000] # Hz
k_values = [0.1, 0.5, 1.0] # 波数
for k in k_values:
print(f"\n波数 k = {k}:")
r = 2.0
theta = sp.pi / 4
pressure = spherical_sound_radiation(r, theta, 2 * sp.pi * 100, k)
pressure_mag = abs(complex(pressure))
print(f" 声压幅值: {pressure_mag:.6f} Pa")
#波数 k = 0.1:
#n=0模式: 幅值 = 44.068679
#n=1模式: 幅值 = 56.516214
#n=2模式: 幅值 = 210.079194
# 声压幅值: 260.355956 Pa
#波数 k = 0.5:
#n=0模式: 幅值 = 1.762747
#n=1模式: 幅值 = 2.260649
#n=2模式: 幅值 = 8.403168
# 声压幅值: 10.414238 Pa
#波数 k = 1.0:
#n=0模式: 幅值 = 0.440687
#n=1模式: 幅值 = 0.565162
#n=2模式: 幅值 = 2.100792
# 声压幅值: 2.603560 Pa
示例5：复数参数和矩阵计算
复数点计算
complex_points = [
0.5 + 0.5j,
1.0 + 0.5j,
-0.5 + 0.5j,
0.5 - 0.5j
]
orders = [0, 1, 2]
第二类Legendre函数在复平面上的值:
for z in complex_points:
print(f"\nz = {z}:")
for n in orders:
Q_n = legendreq(n, z)
print(f" Q_{n}({z}) = {Q_n}")
#z = 0.5+0.5j:
# Q_0((0.5+0.5j)) = 0.4023594781085251+0.5535743588970452j
# Q_1((0.5+0.5j)) = -1.07560744039426+0.47796691850278517j
# Q_2((0.5+0.5j)) = -1.3663605082270465-0.7250175708671287j
#z = (1+0.5j):
# Q_0((1+0.5j)) = 0.708303336014054-0.6629088318340164j
# Q_1((1+0.5j)) = 0.0397577519310621-0.30875716382698926j
# Q_2((1+0.5j)) = -0.06294716724019189-0.10186301587517912j
#z = (-0.5+0.5j):
# Q_0((-0.5+0.5j)) = -0.4023594781085251+0.5535743588970452j
# Q_1((-0.5+0.5j)) = -1.07560744039426-0.47796691850278517j
# Q_2((-0.5+0.5j)) = 1.3663605082270465-0.7250175708671287j
#z = (0.5-0.5j):
# Q_0((0.5-0.5j)) = 0.4023594781085251-0.5535743588970452j
# Q_1((0.5-0.5j)) = -1.07560744039426-0.47796691850278517j
# Q_2((0.5-0.5j)) = -1.3663605082270465+0.7250175708671287j
矩阵输入示例:
n_matrix = [[0, 1], [2, 3]]
x_matrix = [[0.1, 0.2], [0.3, 0.4]]
逐元素矩阵计算:
for i in range(2):
for j in range(2):
n = n_matrix[i][j]
x = x_matrix[i][j]
Q_val = legendreq(n, x)
print(f" Q_{n}({x}) = {Q_val}")
#Q_0(0.1) = 0.10033534773107562
#Q_1(0.2) = -0.9594534891891836
#Q_2(0.3) = -0.5629746555341357
#Q_3(0.4) = 0.08026113738148181
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.special import lqmn
import sympy as sp
def legendre_polynomial_second_kind(input_str):
"""
对标 Wolfram 的 LegendreQ[n, 0, x]，计算第二类 Legendre 多项式（连带次数 m=0）
参数:
input_str: 输入字符串，解析为 (n, x) 的元组，支持标量或矩阵输入
返回:
SymPy 表达式/矩阵 或错误信息字符串
注意:
- 仅支持整数阶 n 和标量 x
- 当 n 为非整数时，符号计算可能不准确
- 矩阵输入需确保形状一致
示例:
>>> legendre_polynomial_second_kind("(2, 0.5)")
0.549306144334055
>>> legendre_polynomial_second_kind("(n, x)")
(log((x + 1)/(1 - x))*LegendreP(n, x)/2 - Sum(LegendreP(k, x)/(n - k), (k, 0, n - 1))/2)
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def legendre_q_scalar(n, x):
"""
计算第二类 Legendre 函数 Q_n(x)（m=0 的情况） - 标量输入版本
参数:
n: 阶数（非负整数）
x: 自变量（实数或复数）
返回:
Q_n(x) 的值
"""
# 确保 n 是整数
if not isinstance(n, (int, np.integer, sp.Integer)):
raise TypeError("阶数 n 必须是整数")
if n < 0:
raise ValueError("阶数 n 必须为非负整数")
# 计算 Q_n(x)
q, _ = lqmn(0, int(n), complex(x))
return q[0, -1]
# 输入必须为 (n, x) 的元组
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
result = legendre_q_scalar(*expr)
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. 基本用法
print("Q0(0.5) =", legendre_polynomial_second_kind("0, 0.5"))
# 0.5493061443340549
print("Q1(1.5) =", legendre_polynomial_second_kind("1, 1.5"))
# 0.20707843432557524
print("Q2(2.5) =", legendre_polynomial_second_kind("2, 2.5"))
# 0.009884255468216039
# 2. 复数输入
print("\n复数输入:")
x_complex = 0.5 + 0.5j
print(f"Q2({x_complex}) =", legendre_polynomial_second_kind(f"2, {x_complex}"))
# (-1.3663605082270465-0.7250175708671287j)
符号函数极限
limit(f,var,a)在var接近a时返回符号表达式f的双向极限.
limit(f,a)使用symvar找到的默认变量.
limit(f)返回的极限为0.
limit(f,var,a,left)在var接近a时返回f的左侧极限.
limit(f,var,a,right)在var接近a时返回f的右侧极限.
f —— 输入,符号表达式,符号函数,符号向量,符号矩阵
var — 自变量,符号变量
a —— 极限点, 数字
示例1：微积分基础与连续性分析
基本三角函数极限
print("sin(x)/x 在 x→0:", limit((sin(x)/x, x, 0))
print("(1-cos(x))/x² 在 x→0:", limit(((1-cos(x))/x**2, x, 0)))
#sin(x)/x 在 x→0: 1.0
#(1-cos(x))/x² 在 x→0: 0.5
指数和对数函数极限
print("(e^x - 1)/x 在 x→0:", limit(((exp(x)-1)/x, x, 0)))
print("ln(1+x)/x 在 x→0:", limit((ln(1+x)/x, x, 0)))
#(e^x - 1)/x 在 x→0: 1.0
#ln(1+x)/x 在 x→0: 1.0
重要极限
print("(1 + 1/x)^x 在 x→∞:", limit(((1 + 1/x)**x, x, oo)))
#(1 + 1/x)^x 在 x→∞: 2.71828182845905
不连续点分析:
print("(x²-1)/(x-1) 在 x→1:", limit(((x**2-1)/(x-1), x, 1)))
#(x²-1)/(x-1) 在 x→1: 2.0
跳跃间断点
print("1/x 在 x→0⁺:", limit((1/x, x, 0, right)))
print("1/x 在 x→0⁻:", limit((1/x, x, 0, left)))
#1/x 在 x→0⁺: oo
#1/x 在 x→0⁻: -oo
无穷间断点
print("1/x² 在 x→0:", limit_equation((1/x**2, x, 0)))
#1/x² 在 x→0: oo
示例2：物理中的极限应用
相对论动能公式在低速下的经典近似
relativistic_energy = m * c ** 2 * (1 / sqrt(1 - v ** 2 / c ** 2) - 1)
classical_approximation = limit(relativistic_energy, v, 0)
print(f"相对论动能: {relativistic_energy}")
print(f"低速近似: {classical_approximation}")
#相对论动能: c**2*m*(-1 + 1/sqrt(1 - v**2/c**2))
#低速近似: 0
简谐运动的小角度近似:
pendulum_equation = sin(theta)
small_angle_approx = limit(pendulum_equation/theta, theta, 0)
print(f"sin(θ) 在 θ→0: {small_angle_approx} (小角度近似: sinθ ≈ θ)")
#sin(θ) 在 θ→0: 1.0(小角度近似: sinθ ≈ θ)
黑体辐射的维恩位移定律推导:
普朗克定律在高温极限
planck_law = 8 * pi * h * c / (lambda_ ** 5 * (exp(h * c / (lambda_ * k * T)) - 1))
rayleigh_jeans = limit(planck_law, T, oo) # 高温极限→瑞利-金斯公式
print(f"高温极限下的辐射公式: {rayleigh_jeans}")
#高温极限下的辐射公式: oo*sign(k/lambda_**4)
示例3：工程控制系统稳定性分析
低通滤波器在截止频率处的增益
lowpass_gain = 1 / sqrt(1 + (f / f_c) ** 2)
gain_at_cutoff = limit(lowpass_gain, f, f_c)
print(f"低通滤波器在截止频率处增益: {gain_at_cutoff} ≈ -3dB")
#低通滤波器在截止频率处增益: 0.707106781186548 ≈ -3dB
高频衰减
high_freq_attenuation = limit(lowpass_gain, f, oo)
print(f"低通滤波器高频衰减: {high_freq_attenuation}")
#低通滤波器高频衰减: 0
示例4：经济学边际分析
弹性分析:
demand_function = 1000 * P ** (-1.5) # 弹性需求函数
价格弹性 = (dQ/dP) * (P/Q)
price_elasticity_def = ((demand_function.subs(P, P + h) - demand_function) / h) * (P / demand_function)
price_elasticity = limit(price_elasticity_def, h, 0)
print(f"需求函数: Q(P) = {demand_function}")
print(f"价格弹性: ε = {price_elasticity.simplify()}")
#需求函数: Q(P) = 1000/P**1.5
#价格弹性: ε = -1.5
示例5：矩阵极限与线性系统分析
矩阵指数函数的极限:
A = [[0, 1], [-1, 0]] # 旋转矩阵
矩阵指数在t→0的极限
I = eye(2)
matrix_exp_approx = I + A * t + (A * t) ** 2 / 2
计算矩阵极限
result = limit(matrix_exp_approx, t, 0)
print("矩阵指数在t→0的极限:")
print(result)
#矩阵指数在t→0的极限:
#[[1.0, 0],
[0, 1.0]]
状态转移矩阵的稳态分析:
马尔可夫链的转移矩阵
P = [[0.8, 0.2], [0.3, 0.7]]
计算P^n在n→∞的极限（稳态分布）
steady_state = limit(P ** n, n, oo)
print("马尔可夫链稳态分布:")
print(steady_state)
#马尔可夫链稳态分布:
#[[0.6, 0.4],
[0.6, 0.4]]
控制系统状态空间分析的极限:
A_sys = [[-1, 2], [0, -3]]
状态转移矩阵在t→∞的极限（稳定性分析）
state_transition = exp(A_sys * t)
steady_state_sys = limit(state_transition, t, oo)
print("系统状态转移矩阵在t→∞:")
print(steady_state_sys)
#系统状态转移矩阵在t→∞:
#[[0, 0],
[0, 0]]
示例6：复杂函数与特殊极限
狄利克雷函数类型的极限
oscillatory_function = sin(1 / x)
print(f"sin(1/x) 在 x→0: {limit((oscillatory_function, x, 0))} (极限不存在)")
#sin(1/x) 在 x→0: AccumBounds(-1, 1) (极限不存在)
但有界振荡
bounded_oscillation = x * sin(1 / x)
print(f"x·sin(1/x) 在 x→0: {limit((bounded_oscillation, x, 0))}")
#x·sin(1/x) 在 x→0: 0
分段函数的极限:
piecewise_func = Piecewise(
(x ** 2, x < 1),
(2 * x - 1, x >= 1)
)
left_limit = limit(piecewise_func, x, 1, left)
right_limit = limit(piecewise_func, x, 1, right)
print(f"分段函数在 x=1 的左极限: {left_limit}")
print(f"分段函数在 x=1 的右极限: {right_limit}")
#分段函数在 x=1 的左极限: 1.0
#分段函数在 x=1 的右极限: 1.0
参数曲线在参数趋近某值时的行为
x_param = t * cos(t)
y_param = t * sin(t)
当t→0时的极限
x_limit = limit(x_param, t, 0)
y_limit = limit(y_param, t, 0)
print(f"参数曲线在 t→0: ({x_limit}, {y_limit})")
#参数曲线在 t→0: (0, 0)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def limit_equation(input_str):
"""
对标 MATLAB 的 limit 函数实现，支持多元极限计算
功能特性:
1. 支持标量/矩阵输入
2. 支持左右极限计算
3. 自动符号变量检测
4. 数值与符号混合计算
参数格式:
"表达式, 变量->趋近点, 方向" 或元组形式 (表达式, 变量, 趋近点, 方向)
返回:
计算结果（标量/矩阵）或错误信息字符串
示例:
>>> limit_equation("(sin(x)/x, x, 0)")
1
>>> limit_equation("([[x,y], [1/x,sin(y)/y]], y, 0)")
Matrix([
[x, 0],
[1/x, 1]])
>>> limit_equation("(1/x, x, 0, 'right')")
∞
"""
try:
# 参数解析与标准化
def parse_parameters(expr):
"""将输入统一转换为四元组 (f, var, a, dir)"""
# 默认参数设置
default_var = None
default_a = 0
default_dir = '+-'
if isinstance(expr, tuple):
params = list(expr)
# 参数数量校验
if len(params) < 2 or len(params) > 4:
raise ValueError("参数数量错误，需要2-4个参数")
# 参数类型推断
f = params[0]
var = params[1] if len(params) >= 2 else default_var
a = params[2] if len(params) >= 3 else default_a
direction = params[3] if len(params) == 4 else default_dir
return f, var, a, direction
else:
# 处理非元组输入
return expr, default_var, default_a, default_dir
# 符号化输入并解析参数
expr = sp.sympify(input_str, evaluate=False)
f_expr, var_expr, a_expr, dir_expr = parse_parameters(expr)
# 主计算逻辑
def compute_limit(f, var, a, direction):
"""核心极限计算函数"""
# 自动检测变量
if var is None:
free_symbols = f.free_symbols
if len(free_symbols) != 1:
raise ValueError("自动变量检测失败：表达式包含多个自由符号")
var = next(iter(free_symbols))
# 方向参数标准化
dir_map = {'left': '-', 'right': '+', None: None}
direction = dir_map.get(str(direction).lower(), direction)
# 执行极限计算
try:
return sp.limit(f, var, a, dir=direction)
except (NotImplementedError, ValueError) as e:
return f"计算错误: {str(e)}"
# 分支处理逻辑
result = None
# 情况2：符号表达式输入
if isinstance(f_expr, sp.Expr):
result = compute_limit(f_expr, var_expr, a_expr, dir_expr)
else:
raise TypeError("不支持的输入类型")
return result.evalf() if isinstance(result, sp.Expr) else result
except Exception as e:
return f"错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 基础标量计算
print(limit_equation("(sin(x)/x, x, 0)"))
# 1.00000000000000
print(limit_equation("(exp(1/x), x, 0, 'right')"))
# oo
# 自动变量检测
print(limit_equation("(sin(z)/z, z, 0)"))
# 1.00000000000000
print(limit_equation("(sin(x)/x)"))
# 1.00000000000000
线性方程组求解
X = linsolve(A,B) 使用以下方法之一求解线性方程组 AX = B：
当A是方阵时,linsolve使用LU分解和部分主元消元法.
对于所有其他情况,linsolve使用QR分解和列主元消元法.
如果A为病态(对于方阵)或秩亏(对于矩形矩阵),则linsolve发出警告.
A - 系数矩阵, 符号矩阵
B - 输入数组, 向量,矩阵.
X - 输出,向量,矩阵.
示例1. 电路分析（基尔霍夫定律）
电路方程：根据基尔霍夫定律
节点电流方程: I1 = I2 + I3
回路1: 10 - 2*I1 - 3*I2 = 0
回路2: 3*I2 - 4*I3 - 5 = 0
test_input = (
[I1 - I2 - I3, # 节点方程
2 * I1 + 3 * I2 - 10, # 回路1
3 * I2 - 4 * I3 - 5], # 回路2
[I1, I2, I3]
)
result = linsolve(test_input)
print("电路电流分析:")
print(f"电流解: {result}")
#电路电流分析:
#电流解: {(55/26, 25/13, 5/26)}
示例2. 结构力学（桁架受力分析）
F1, F2, F3 # 各杆件受力
桁架节点平衡方程
水平平衡: F1*cos(45) + F2 = 0
垂直平衡: F1*sin(45) + F3 - 1000 = 0
力矩平衡: F3*2 - 1000*1 = 0
cos45 = sqrt(2) / 2
sin45 = sqrt(2) / 2
test_input = (
[F1 * cos45 + F2, # 水平平衡
F1 * sin45 + F3 - 1000, # 垂直平衡
2 * F3 - 1000], # 力矩平衡
[F1, F2, F3]
)
result = linsolve(test_input)
print("\n桁架受力分析:")
print(f"杆件受力: {result}")
#桁架受力分析:
#杆件受力: {(500*sqrt(2), -500, 500)}
示例3. 经济学（市场均衡模型）
P, Qd, Qs # 价格，需求量，供给量
需求函数: Qd = 100 - 2P
供给函数: Qs = 20 + 3P
市场均衡: Qd = Qs
test_input = (
[Qd - 100 + 2 * P, # 需求函数
Qs - 20 - 3 * P, # 供给函数
Qd - Qs], # 市场均衡条件
[P, Qd, Qs]
)
result = linsolve(test_input)
print("\n市场均衡分析:")
print(f"均衡价格和数量: {result}")
#市场均衡分析:
#均衡价格和数量: {(16, 68, 68)}
示例4. 化学计量学（化学反应配平）
配平反应: a C2H6 + b O2 → c CO2 + d H2O
a, b, c, d = sp.symbols('a b c d')
碳平衡: 2a = c
氢平衡: 6a = 2d
氧平衡: 2b = 2c + d
test_input = (
[2 * a - c, # 碳原子平衡
6 * a - 2 * d, # 氢原子平衡
2 * b - 2 * c - d], # 氧原子平衡
[a, b, c, d]
)
result = linsolve(test_input)
print("\n化学反应配平:")
print(f"化学计量系数: {result}")
#化学反应配平:
#化学计量系数: {(tau0/3, 7*tau0/6, 2*tau0/3, tau0)}
示例5. 网络流量分析
节点流量守恒方程
节点A: f1 = f2 + f3
节点B: f2 + f4 = 50
节点C: f3 = 30 + f4
外部输入: f1 = 100
test_input = (
[f1 - f2 - f3, # 节点A
f2 + f4 - 50, # 节点B
f3 - 30 - f4, # 节点C
f1 - 100], # 外部输入
[f1, f2, f3, f4]
)
result = linsolve(test_input)
print("\n网络流量分析:")
print(f"各支路流量: {result}")
#网络流量分析:
#各支路流量: EmptySet
示例6. 符号参数系统（通用性测试）
带参数的线性系统
test_input = (
[a * x + b * y - c,
b * x + c * y - a,
x + y - z],
[x, y, z]
)
result = linsolve(test_input)
print("\n符号参数系统:")
print(f"符号解: {result}")
#符号参数系统:
#符号解: {((-a*b + c**2)/(a*c - b**2), (a**2 - b*c)/(a*c - b**2),
(a**2 - a*b - b*c + c**2)/(a*c - b**2))}
示例7. 数值稳定性测试
测试病态矩阵（接近奇异的矩阵）
良态系统
result1 = linsolve([[1, 2], [2, 3]], [5, 8])
print(f"良态系统解: {result1}")
#良态系统解: {(1, 2)}
接近奇异的系统
epsilon = 1e-10
result2 = linsolve([[1, 1], [1, 1+epsilon]], [2, 2+epsilon])
print(f"接近奇异系统解: {result2}")
#接近奇异系统解: {(1.0, 1.0)}
示例8. 无解和无穷解情况
无解系统
result1 = linsolve([[1, 1], [1, 1]], [1, 2])
print(f"无解系统: {result1}")
#无解系统: EmptySet
无穷多解系统
result2 = linsolve("([[1, 1], [2, 2]], [2, 4])")
print(f"无穷多解系统: {result2}")
#无穷多解系统: {(2 - tau0, tau0)}
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Matrix, linsolve, linear_eq_to_matrix, Eq, sympify
def linear_equations_solve(input_str):
"""
改进后的线性方程组求解函数，支持符号输入及方程表达式直接输入
"""
try:
# 符号化输入并标准化处理
expr = sympify(input_str)
def parse_input(expr):
"""输入解析标准化为 (A, b) 形式"""
# 增广矩阵处理
M = sp.Matrix(expr) if isinstance(expr, list) else None
if M is not None:
if M.cols < 2:
raise ValueError("增广矩阵至少需要2列")
return M[:, :-1], M[:, -1]
# 元组输入处理
if isinstance(expr, tuple):
if len(expr) != 2:
raise ValueError("输入格式应为 (方程列表, 变量列表) 或 (系数矩阵, 常数项)")
equations_part, vars_part = expr
# 先判断是否为 (方程表达式, 变量列表) 形式
if all(isinstance(v, sp.Symbol) for v in vars_part):
# 检查 equations_part 是否真的是方程表达式列表
if all(isinstance(eq, (sp.Expr, sp.Eq)) for eq in equations_part):
# 将表达式列表转换为等式列表（假设表达式右边为0）
eq_list = [Eq(eq, 0) if not isinstance(eq, sp.Eq) else eq for eq in equations_part]
A, b = linear_eq_to_matrix(eq_list, vars_part)
return A, b
# 否则按 (系数矩阵, 常数项) 处理
A = Matrix(equations_part) if isinstance(equations_part, list) else equations_part
b = Matrix(vars_part) if isinstance(vars_part, list) else vars_part
if not isinstance(A, Matrix) or not isinstance(b, Matrix):
raise TypeError("系数矩阵和常数项需要为矩阵或可转换形式")
return A, b
raise TypeError("不支持的输入格式")
A, b = parse_input(expr)
# 维度校验
if A.rows != b.rows:
raise ValueError(f"维度不匹配: A({A.shape}) vs b({b.shape})")
# 符号解计算
return linsolve((A, b))
except Exception as e:
return f"求解错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
x, y, z = sp.symbols('x y z')
# 正确输入测试（方程表达式列表 + 变量列表）
test_input = (
[3 * x + 2 * y - z - 1,
2 * x - 2 * y + 4 * z + 2,
2 * x - y + 2 * z],
[x, y, z]
)
print(linear_equations_solve(str(test_input)))
# {(1, -2, -2)}
# 数值解示例（增广矩阵）
print(linear_equations_solve("[[3, 2, -1, 1], [2, -2, 4, -2], [2, -1, 2, 0]]"))
# {(1, -2, -2)}
# 符号解示例（系数矩阵 + 常数项）
print(linear_equations_solve("([[1, 1], [1, -1]], [5, 1])"))
# {(3, 2)}
点线图
ListPlot(x,(y)通过将抽象数字映射为几何图形，赋予数据空间意义.直观、快速地将离散数据点可视化.
趋势识别：销售数据监控
X = [1, 2, 3, 4, 5] (代表星期一至星期五)
Y = [1, 2, 3, 2, 1] (代表每日销量)
图的作用：一眼看出销量在周中（星期三）达到峰值，周初和周末较低。结论：可能需要加强周末促销或分析周中高峰原因。
异常点检测：传感器读数
X = [1, 2, 3, 4, 5] (连续5个时间点的采样)
Y = [10, 10.2, 10.1, 15.5, 10.3] (温度传感器读数)
图的作用：第4个点 (15.5) 明显脱离其他点形成的“平坦线”。结论：该时间点可能存在传感器瞬时故障、外部干扰或真实异常事件，需要检查。
算法迭代过程：优化算法收敛
X = [1, 2, 3, 4, 5] (迭代次数)
Y = [100, 50, 25, 12, 6] (目标函数值，如误差、成本)
图的作用：清晰显示误差随迭代次数快速下降并趋于稳定。结论：算法收敛良好；可能在第4/5次迭代后即可停止，节省计算资源。
简单系统响应：弹簧振子位移
X = [0, 1, 2, 3, 4] (时间 t)
Y = [0, 2, 0, -2, 0] (位移 s)
图的作用：点描绘出一个振荡衰减（或理想简谐运动）的轨迹。结论：直观验证了物理模型的预期行为（周期性运动）。
A/B测试结果：用户点击率对比
数据组A (旧方案)： Y_A = [5.1, 5.3, 5.0, 5.2, 5.1] (%)
数据组B (新方案)： Y_B = [5.5, 6.0, 5.8, 5.7, 6.2] (%)
图的作用：将两条线画在同一图上。结论：新方案（B线）稳定且显著高于旧方案（A线），初步证明新方案有效。
资源消耗监控：服务器内存占用
X = [9:00, 10:00, 11:00, 12:00, 13:00] (时间)
Y = [45, 48, 65, 62, 50] (内存使用率 %)
图的作用：显示在上午11点和12点出现明显峰值。结论：该时段负载较高，需关注是否常态或进行扩容/优化。
离散事件记录：网站每日故障次数
X = [1, 2, 3, 4, 5] (日期)
Y = [0, 0, 1, 0, 3] (故障次数)
图的作用：大部分日期无故障（Y=0），第3天有1次小故障，第5天发生3次故障（显著异常）。结论：需紧急排查第5天的系统问题。
三维点线图
ListPlot3D(f1(x),f2(x)) 将离散的三维数据点集转化为空间中的可视化图形，揭示复杂的高维关系.
对数螺旋线: 半径对数增长，避免大尺度发散
ListPlot3D(ln(a+1)*cos(a),ln(a+1)*sin(a),a=[0.1,50,500])
阻尼振荡弹簧: 振幅指数衰减的压缩螺旋
ListPlot3D(a*cos(a)*exp(-0.01a),a*sin(a)*exp(-0.01a),0.2a,a=[0.1,50,500])
双频率扭曲曲线: XY 平面螺旋 + Z 方向正弦振荡
ListPlot3D(a*cos(0.5a),a*sin(0.8a),10sin(0.1a),a=[0.1,50,500])
球面螺旋: 球坐标系下的缓慢展开螺旋
ListPlot3D(a*sin(0.02a)*cos(a),a*sin(0.02a)*sin(a),a*cos(0.02a),a=[0.1,50,500])
锥形莫比乌斯带: 经典莫比乌斯带叠加径向缩放
ListPlot3D((1+0.3cos(1.5a))*cos(a),(1+0.3cos(1.5a))*sin(a),0.5sin(1.5a)*log(a+1),a=[0.1,50,500])
指数增长螺旋线 (半径指数爆炸): 两个带电粒子在强电场和磁场耦合作用下的分离轨迹。
ListPlot3D(exp(0.1*a)*a,exp(0.1*a)*(-a),a,a=[0.1,50,500])
对数收敛空间曲线 (半径对数饱和): 演示非线性函数组合对空间形态的显著影响。
ListPlot3D(ln(a+10)*a,ln(a+10)*sqrt(a),0.5*a^2,a=[1,50,500])
幂律组合弹簧 (振幅幂律衰减): 更真实的阻尼模型（指数衰减是理想特例）。
ListPlot3D(a*exp(-0.02*(a^0.8)),(a^1.2)*exp(-0.02*a^0.8),0.3a,a=[0.5,50,500])
双曲螺旋线 (渐近线约束): 粒子在特定势场（如重力+斥力）中的逃逸轨迹。
ListPlot3D(sinh(0.1a),a*cosh(0.1a),0.2a,a=[0.1,50,500])
自然对数
Y = log(X)返回数组X中每个元素的自然对数ln(x).
Y = log(b, X)返回数组X中每个元素的以b为底数的对数ln(x).
log函数的域包含负数和复数,如果使用不当,可能会导致意外结果.对于负数和复数z=u + i*w,复数对数log(z)返回log(abs(z)) + 1i*angle(z)
X — 输入数组标量,向量,矩阵
示例1. 金融计算（复利和连续复利）
连续复利计算
principal = 1000 # 本金
annual_rate = 0.05 # 年利率5%
time = 3 # 3年
连续复利公式: A = P * e^(rt)
计算终值
future_value = principal * exp(annual_rate * time)
print(f"连续复利终值: {future_value:.2f}")
#连续复利终值: 1161.83
反过来计算利率
calculated_rate = log(future_value / principal) / time
print(f"计算出的年利率: {calculated_rate:.4f}")
#计算出的年利率: 0.0500
示例2. 物理和工程（放射性衰变）
半衰期计算
half_life = 5730 # 碳14半衰期（年）
decay_constant = log(2) / half_life
print(f"衰变常数: {decay_constant}")
#衰变常数: 0.000120968094338559
计算经过时间后的剩余比例
time_elapsed = 1000 # 经过1000年
remaining_ratio = exp(-decay_constant * time_elapsed)
print(f"1000年后剩余比例: {remaining_ratio:.4f}")
#1000年后剩余比例: 0.8861
示例3. 生物学（种群增长模型）
指数增长模型: N(t) = N0 * e^(rt)
initial_population = 100
growth_rate = 0.1 # 增长率10%
time_periods = [1, 5, 10, 20] # 不同时间点
for t in time_periods:
population = initial_population * exp(growth_rate * t)
print(f"经过 {t} 个时间单位后种群数量: {population:.1f}")
#经过 1 个时间单位后种群数量: 110.5
#经过 5 个时间单位后种群数量: 164.9
#经过 10 个时间单位后种群数量: 271.8
#经过 20 个时间单位后种群数量: 738.9
计算达到特定种群规模所需时间
target_population = 500
required_time = log(target_population / initial_population) / growth_rate
print(f"达到500个个体所需时间: {required_time:.2f}")
##达到500个个体所需时间: 16.09
示例4. 信息论（熵计算）
概率分布
probabilities = [0.25, 0.25, 0.25, 0.25] # 均匀分布
信息熵公式: H = -Σ p_i * ln(p_i)
entropy = -sum(p * log(p) for p in probabilities)
print(f"均匀分布的熵: {entropy:.4f}")
#均匀分布的熵: 1.3863
非均匀分布
prob_nonuniform = [0.5, 0.25, 0.125, 0.125]
entropy_nonuniform = -sum(p * log(p) for p in prob_nonuniform)
print(f"非均匀分布的熵: {entropy_nonuniform:.4f}")
#非均匀分布的熵: 1.2130
示例5. 化学（pH值和酸度计算）
H+离子浓度
H_plus_concentrations = [1e-7, 1e-3, 1e-9] # mol/L
for conc in H_plus_concentrations:
pH = -log(conc) / log(10)
print(f"[H+] = {conc:.1e} mol/L, pH = {pH:.2f}")
#[H+] = 1.0e-07 mol/L, pH = 7.00
#[H+] = 1.0e-03 mol/L, pH = 3.00
#[H+] = 1.0e-09 mol/L, pH = 9.00
计算pOH
OH_minus = 1e-5 # mol/L
pOH = -log(OH_minus) / log(10)
print(f"pOH值: {pOH:.2f}")
#pOH值: 5.00
示例6. 统计学（正态分布）
正态分布概率密度函数中的自然对数部分
x = 1.5 # 随机变量值
mu = 0 # 均值
sigma = 1 # 标准差
对数似然函数
log_likelihood = -0.5 * log(2) - 0.5 * log(2 * pi) \
- log(sigma) - 0.5 * ((x - mu) / sigma) ** 2
print(f"在x={x}处的对数似然: {log_likelihood:.4f}")
#在x=1.5处的对数似然: -2.3905
示例7. 矩阵运算（多元统计分析）
协方差矩阵的自然对数（在多元正态分布中用到）
covariance_matrix = [[4, 2], [2, 3]]
print("原始协方差矩阵:")
print(log(covariance_matrix))
#原始协方差矩阵:
#[[1.38629436111989, 0.693147180559945],
[0.693147180559945, 1.09861228866811]]
相关性矩阵
correlation_matrix = [[1, 0.5], [0.5, 1]]
print("相关性矩阵的对数:")
print(log(correlation_matrix))
#相关性矩阵的对数:
#[[0, -0.693147180559945],
[-0.693147180559945, 0]]
对角矩阵
diagonal_matrix = [[2, 0], [0, 3]]
print("对角矩阵的对数:")
print(log(diagonal_matrix))
#对角矩阵的对数:
#[[0.693147180559945, -oo],
[-oo, 1.09861228866811]]
示例8. 符号计算应用
符号表达式的自然对数
expr1 = x**2 + 1
result1 = log(expr1)
print(f"ln({expr1}) = {result1}")
#ln(x**2 + 1) = log(x**2 + 1)
复合函数
expr2 = exp(x) + sin(x)
result2 = log(expr2)
print(f"ln({expr2}) = {result2}")
#ln(exp(x) + sin(x)) = log(exp(x) + sin(x))
多变量函数
expr3 = x*y
result3 = log(expr3)
print(f"ln({expr3}) = {result3}")
#ln(x*y) = log(x*y)
示例9. 工程应用（信号处理）
分贝计算
power_ratio = 100 # 功率比
dB = 10 * log(power_ratio) / log(10)
print(f"功率比{power_ratio}对应的分贝值: {dB:.2f} dB")
#功率比100对应的分贝值: 20.00 dB
电压比的分贝计算
voltage_ratio = 10
dB_voltage = 20 * log(voltage_ratio) / log(10)
print(f"电压比{voltage_ratio}对应的分贝值: {dB_voltage:.2f} dB")
#电压比10对应的分贝值: 20.00 dB
示例10. 复杂数值处理
负数的自然对数（在复数域中）
negative_numbers = [-1, -2, -5]
for num in negative_numbers:
result = log(num)
print(f"ln({num}) = {result}")
#ln(-1) = 3.14159265358979*I
#ln(-2) = 0.693147180559945 + 3.14159265358979*I
#ln(-5) = 1.6094379124341 + 3.14159265358979*I
复数输入
complex_numbers = [1+1@i, 2-3@i, -1+2@i]
for num in complex_numbers:
result = log(num)
print(f"ln({num}) = {result}")
#ln(1+1i) = 0.346573590279973 + 0.785398163397448*I
#ln(2-3i) = 1.28247467873077 - 0.982793723247329*I
#ln(-1+2i) = 0.80471895621705 + 2.0344439357957*I
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def logarithmic_natural(input_str):
"""
对标MATLAB的log函数，计算自然对数
参数:
input_str: 输入的字符串表达式，可以是标量、矩阵或SymPy表达式
返回:
自然对数结果，矩阵或标量。若输入错误返回错误信息。
示例:
logarithmic_natural('5') -> ln(5)的数值
logarithmic_natural('[[1, 2], [3, 4]]') -> 各元素自然对数的矩阵
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
if expr.is_number:
z = complex(expr)
result = np.log(z)
elif expr.free_symbols:
# 处理标量或符号表达式
result = sp.log(expr, sp.E)
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(logarithmic_natural('2'))
# (0.6931471805599453+0j)
针对较小的X精确计算1+X的自然对数
Y = log1p(X) 计算数组X中每个元素的自然对数log(1+X),但不显式计算1+X.如果X<-1,则Y为复数.
此函数对于X中的小实数值更精确,因为它会补偿1+X中的舍入误差.
X — 输入数组标量,向量,矩阵
示例1. 数值稳定性测试（小值场景）
测试非常小的值
small_values = [1e-10, 1e-15, 1e-20]
for val in small_values:
# 直接计算 log(1+x)
direct_result = log1p(val)
# 对比传统方法
x = float(val)
traditional = log(1 + x)
print(f"log1p({val}) = {direct_result}")
print(f"log(1+{val}) = {traditional}")
print(f"相对误差: {abs((direct_result - traditional) / direct_result)}")
#log1p(1e-10) = 1.00000008269037E-10
#log(1+1e-10) = 1.000000082690371e-10
#相对误差: 0
#log1p(1e-15) = 1.11022302462516E-15
#log(1+1e-15) = 1.110223024625156e-15
#相对误差: 0
#log1p(1e-20) = 0
#log(1+1e-20) = 0.0
#相对误差: nan
示例2. 概率计算（对数概率和 softmax）
在机器学习中，经常需要计算 log(1 + exp(x)) 来避免数值溢出
这可以通过 log1p(exp(x)) 来实现
x_values = [-100, -10, -1, 0, 1, 10, 100]
for x in x_values:
# 计算 log(1 + exp(x))，数值稳定的方式
result = log1p(exp(x))
print(f"log(1 + exp({x})) = {result}")
#log(1 + exp(-100)) = 0
#log(1 + exp(-10)) = 0.0000453988992168705
#log(1 + exp(-1)) = 0.313261687518223
#log(1 + exp(0)) = 0.693147180559945
#log(1 + exp(1)) = 1.31326168751822
#log(1 + exp(10)) = 10.0000453988992
#log(1 + exp(100)) = 100.000000000000
示例3. 金融计算（连续复利的微小利率）
当利率非常小时，使用 log1p 可以提高精度
daily_rates = [0.0001, 0.00001, 0.000001] # 日利率
for rate_str in daily_rates:
rate = float(rate_str)
# 计算年化利率（假设365天）
annual_rate_log1p = 365 * log1p(rate_str)
# 传统方法
annual_rate_traditional = 365 * log(1 + rate)
print(f"日利率 {rate:.6f}:")
print(f" log1p方法年化: {annual_rate_log1p:.8f}")
print(f" 传统方法年化: {annual_rate_traditional:.8f}")
print(f" 差值: {abs(annual_rate_log1p - annual_rate_traditional)}")
#日利率 0.000100:
# log1p方法年化: 0.03649818
# 传统方法年化: 0.03649818
# 差值: 0
#日利率 0.000010:
# log1p方法年化: 0.00364998
# 传统方法年化: 0.00364998
# 差值: 0
#日利率 0.000001:
# log1p方法年化: 0.00036500
# 传统方法年化: 0.00036500
# 差值: 0
示例4. 统计学（对数似然函数）
在最大似然估计中，经常需要计算 log(1 + small_probability)
probabilities = [0.001, 0.0001, 0.00001]
for prob_str in probabilities:
prob = float(prob_str)
# 计算对数似然
log_likelihood = log1p(prob_str)
# 对比传统方法
traditional = log(1 + prob)
print(f"概率 {prob:.6f}:")
print(f" log1p结果: {log_likelihood:.10f}")
print(f" 传统结果: {traditional:.10f}")
if prob < 1e-4:
relative_error = abs((log_likelihood - traditional) / log_likelihood)
print(f" 相对误差: {relative_error}")
#概率 0.001000:
# log1p结果: 0.0009995003
# 传统结果: 0.0009995003
#概率 0.000100:
# log1p结果: 0.0000999950
# 传统结果: 0.0000999950
#概率 0.000010:
# log1p结果: 0.0000100000
# 传统结果: 0.0000100000
# 相对误差: 0
示例5. 物理计算（微小相对论效应）
相对论中的洛伦兹因子：γ = 1/sqrt(1 - v²/c²)
当 v << c 时，可以使用近似展开
beta_values = [0.001, 0.0001, 0.00001] # v/c
for beta_str in beta_values:
beta = float(beta_str)
beta_sq = beta ** 2
# 计算 log(γ) = -0.5 * log(1 - β²)
# 使用 log1p 处理接近 1 的值
log_gamma = -0.5 * log1p(-beta_sq)
print(f"β = {beta:.6f}: log(γ) = {log_gamma:.12f}")
#β = 0.001000: log(γ) = 0.0000005
#β = 0.000100: log(γ) = 0.000000005
#β = 0.000010: log(γ) = 0.00000000005
示例6. 矩阵计算（特征值处理）
接近单位矩阵的情况
small_perturbation = [[0.001, 0.0005], [0.0005, 0.001]]
print("小扰动矩阵的 log1p:")
result = log1p(small_perturbation)
#小扰动矩阵的 log1p:
#[[0.000999500333083423, 0.000499875041650993],
[0.000499875041650993, 0.000999500333083423]]
对角元素接近零的矩阵
near_zero_matrix = [[1e-10, 2e-10], [3e-10, 4e-10]]
print("\n近零矩阵的 log1p:")
result2 = log1p(near_zero_matrix)
#近零矩阵的 log1p:
#[[1.00000008269037e-10, 2.00000016528074e-10],
[3.00000024777111e-10, 4.00000033016148e-10]]
示例7. 复数运算
小复数
small_complex = [1e-10+2e-10@i, 1e-15+1e-15@i]
for comp_str in small_complex:
result = log1p(comp_str)
print(f"log1p({comp_str}) = {result}")
#log1p(1e-10+2e-10i) = 1.00000008269037e-10 + 1.9999999998e-10*I
#log1p(1e-15+1e-15i) = 1.11022302462516e-15 + 9.99999999999999e-16*I
对比传统方法
z = 1e-15 + 1e-15j
traditional = np.log(1 + z)
log1p_result = np.log1p(z)
print(f"\n传统方法: log(1+{z}) = {traditional}")
print(f"log1p方法: log1p({z}) = {log1p_result}")
#传统方法: log(1+(1e-15+1e-15j)) = (1.1102230246251563e-15+9.999999999999989e-16j)
#log1p方法: log1p((1e-15+1e-15j)) = (1.110223024625156e-15+9.999999999999989e-16j)
示例8. 符号计算和极限情况
计算 log1p(x) 在 x=0 附近的泰勒展开
expr = log1p(x)
taylor_expansion = series(expr, x, 0, 6)
print(f"log1p(x) 的泰勒展开: {taylor_expansion}")
#log1p(x) 的泰勒展开: x - x**2/2 + x**3/3 - x**4/4 + x**5/5 + O(x**6)
极限情况
limit_at_zero = limit(expr / x, x, 0)
print(f"lim(x→0) log1p(x)/x = {limit_at_zero}")
#lim(x→0) log1p(x)/x = 1
示例9. 实际工程应用（传感器数据处理）
假设传感器读数有微小偏差
sensor_offsets = [0.0001, 0.00005, 0.00001]
for offset_str in sensor_offsets:
offset = float(offset_str)
# 计算对数变换后的值（常用于方差稳定化）
log_transformed = log1p(offset_str)
print(f"传感器偏移 {offset:.6f}: 对数变换后 = {log_transformed:.10f}")
# 计算信噪比的改进
if offset > 0:
snr_improvement = log_transformed / offset
print(f" 信噪比改进因子: {snr_improvement:.6f}")
#传感器偏移 0.000100: 对数变换后 = 0.0000999950
# 信噪比改进因子: 0.999950
#传感器偏移 0.000050: 对数变换后 = 0.0000499988
# 信噪比改进因子: 0.999975
#传感器偏移 0.000010: 对数变换后 = 0.0000100000
# 信噪比改进因子: 0.999995
示例10. 误差分析和精度验证
测试各种边界情况
test_cases = [
0, # 零
1e-100, # 极小数
-0.5, # 负值
1e-10+1e-10@i, # 小复数
0.9999999999 # 接近1
]
for case in test_cases:
try:
result = log1p(case)
print(f"log1p({case}) = {result}")
except Exception as e:
print(f"log1p({case}) 错误: {e}")
#log1p(0) = 0
#log1p(1e-100) = 0
#log1p(-0.5) = -0.693147180559945
#log1p(1e-10+1e-10j) = 1.00000008269037e-10 + 9.999999999e-11*I
#log1p(0.9999999999) = 0.693147180509945
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.codegen.cfunctions import log1p
import numpy as np
def logarithmic_natural_1p(input_str):
"""
对标MATLAB的log1p函数，精确计算 log(1+X)
特别适用于X接近0时的小值场景
参数:
input_str: 输入的字符串表达式，可以是标量、矩阵或SymPy表达式
返回:
log(1+x) 计算结果，矩阵或标量。若输入错误返回错误信息
示例:
logarithmic_natural_1p('1e-20') -> 结果约等于1e-20
logarithmic_natural_1p('[[0, 0.5], [1, 2]]') -> 各元素计算log1p的矩阵
"""
try:
expr = sp.sympify(input_str) # 转换为SymPy表达式
error = False
result = None
# 处理标量或符号表达式
if expr.is_number:
z = complex(expr)
result = np.log1p(z)
elif expr.free_symbols:
result = log1p(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 小值输入:")
print(logarithmic_natural_1p('2+2j'))
# (1.2824746787307684+0.5880026035475675j)
# 示例2: 符号表达式输入
print("\n示例2 符号表达式:")
x = sp.Symbol('x')
print(logarithmic_natural_1p('x'))
# log1p(x)
对数积分
A = logint(x)计算对数积分函数(积分对数)
logint根据您使用的参数返回浮点或精确的符号结果.
计算这些数字的整数对数.因为这些数字不是符号对象,所以logint返回浮点结果.
x —— 输入,符号数,符号变量,符号表达式,符号函数,符号向量,符号矩阵
示例1. 数论应用（素数分布）
素数定理：π(x) ~ li(x)，其中π(x)是小于x的素数个数
x_values = [10, 100, 1000, 10000, 100000]
print("x\t\tli(x)\t\t实际素数个数\t相对误差")
print("-" * 60)
for x in x_values:
li_x = logint(x)
# 实际素数个数（近似值）
if x == 10:
actual_primes = 4
elif x == 100:
actual_primes = 25
elif x == 1000:
actual_primes = 168
elif x == 10000:
actual_primes = 1229
elif x == 100000:
actual_primes = 9592
error = abs(li_x - actual_primes) / actual_primes
print(f"{x}\t\t{float(li_x):.2f}\t\t{actual_primes}\t\t{error:.4f}")
#x li(x) 实际素数个数 相对误差
#------------------------------------------------------------
#10 6.17 4 0.5414
#100 30.13 25 0.2050
#1000 177.61 168 0.0572
#10000 1246.14 1229 0.0139
#100000 9629.81 9592 0.0039
示例2. 物理学应用（热传导问题）
在某些热传导问题中，温度分布涉及对数积分
positions = [0.1, 0.5, 1.0, 2.0, 5.0]
print("位置\t\t温度分布(相对值)")
print("-" * 40)
for pos in positions:
# 简化模型：T(r) ∝ li(exp(-r²/4αt))
r_squared = pos ** 2
argument = exp(-r_squared) # 假设 αt = 1
temp_dist = logint(argument)
print(f"{pos}\t\t{float(temp_dist):.6f}")
#位置 温度分布(相对值)
#----------------------------------------
#0.1 -4.037930
#0.5 -1.044283
#1.0 -0.219384
#2.0 -0.003779
#5.0 -0.000000
示例3. 电磁学应用（场强计算）
在特定电荷分布下，电势可能涉及对数积分
distances = [1.0, 2.0, 5.0, 10.0]
print("距离\t\t电势(相对值)")
print("-" * 35)
for d in distances:
# 简化的电势模型：V(r) ∝ li(1/r)
argument = 1/d
potential = logint(argument)
print(f"{d}\t\t{float(potential):.6f}")
#距离 电势(相对值)
#-----------------------------------
#1.0 -inf
#2.0 -0.378671
#5.0 -0.085126
#10.0 -0.032390
示例4. 统计学应用（极值分布）
在某些极值分布中，累积分布函数涉及对数积分
x_values = [1.5, 2.0, 3.0, 5.0, 10.0]
print("x\t\t累积分布函数值")
print("-" * 35)
for x in x_values:
# 简化的极值分布模型
cdf_value = logint(x) / logint(10)
print(f"{x}\t\t{float(cdf_value):.6f}")
#x 累积分布函数值
#-----------------------------------
#1.5 0.020284
#2.0 0.169515
#3.0 0.350913
#5.0 0.589495
#10.0 1.000000
示例5. 工程学应用（信号衰减）
# 在波导或传输线中，信号衰减可能涉及对数积分
frequencies = [1e6, 10e6, 100e6, 1e9] # Hz
print("频率(Hz)\t\t衰减系数(相对值)")
print("-" * 45)
for freq in frequencies:
# 简化的衰减模型：α(f) ∝ li(f/f₀)
f_normalized = freq / 1e6 # 归一化频率
argument = f_normalized
attenuation = logint(argument)
print(f"{freq:.1e}\t\t{float(attenuation):.6f}")
#频率(Hz) 衰减系数(相对值)
#---------------------------------------------
#1.0e+06 -inf
#1.0e+07 6.165600
#1.0e+08 30.126142
#1.0e+09 177.609658
示例6. 矩阵计算（多变量系统）
对角矩阵
diagonal_matrix = [[2, 0, 0], [0, 3, 0], [0, 0, 5]]
print("对角矩阵的对数积分:")
result1 = logint(diagonal_matrix)
print(result1)
#对角矩阵的对数积分:
#[[1.04516378011749, 0, 0],
[0, 2.16358859466719, 0],
[0, 0, 3.63458831003265]]
一般矩阵
general_matrix = [[1.5, 0.5], [0.5, 2.5]]
print("\n一般矩阵的对数积分:")
result2 = logint(general_matrix)
print(result2)
#[[0.125064986315296, -0.378671043061088],
[-0.378671043061088, 1.66729466750632]]
示例7. 复数域对数积分
复数参数
complex_values = [1+1@i, 2-1@i, 0.5+0.5@i]
for comp in complex_values:
result = logint(comp)
print(f"li({comp}) = {result}")
#li(1+1j) = 0.613911669221196 + 2.05958421419258*I
#li(2-1j) = 1.4112590420178 - 1.2247069384103*I
#li(0.5+0.5j) = -0.0139535417669452 + 2.62968411537039*I
示例8. 符号计算和解析性质
基本性质
expr1 = logint(x)
print(f"li(x) = {expr1}")
#li(x) = li(x)
导数
derivative = diff(expr1, x)
print(f"d(li(x))/dx = {derivative}")
#d(li(x))/dx = 0
级数展开
series_expansion = series(expr1, x, 1, 6)
print(f"li(x) 在 x=1 附近的级数展开: {series_expansion}")
#li(x) 在 x=1 附近的级数展开: li(x)
示例9. 特殊值验证
已知的特殊值
special_cases = [
1,
E,
2,
]
for input_val in special_cases:
result = logint(input_val)
print(f"li({input_val}) = {result}")
#li(1) = -oo
#li(E) = 1.89511781635594
#li(2) = 1.04516378011749
示例10. 与指数积分的关系
对数积分与指数积分的关系: li(x) = Ei(ln(x))
x_values = [1.5, 2.0, 3.0, 5.0]
for x in x_values:
li_direct = logint(x)
# 通过指数积分计算
ei_result = Ei(log(x))
print(f"x = {x}:")
print(f" 直接计算 li(x) = {li_direct}")
print(f" 通过 Ei(ln(x)) = {ei_result}")
print(f" 差值: {abs(li_direct - ei_result)}")
#x = 1.5:
# 直接计算 li(x) = 0.125064986315296
# 通过 Ei(ln(x)) = 0.125064986315296
# 差值: 0
#x = 2.0:
# 直接计算 li(x) = 1.04516378011749
# 通过 Ei(ln(x)) = 1.04516378011749
# 差值: 2.22044604925031E-16
#x = 3.0:
# 直接计算 li(x) = 2.16358859466719
# 通过 Ei(ln(x)) = 2.16358859466719
# 差值: 4.44089209850063E-16
#x = 5.0:
# 直接计算 li(x) = 3.63458831003265
# 通过 Ei(ln(x)) = 3.63458831003265
# 差值: 4.44089209850063E-16
示例11. 渐近行为分析
大x时的渐近行为: li(x) ~ x/ln(x)
large_values = [100, 1000, 10000, 100000]
print("x\t\tli(x)\t\tx/ln(x)\t\t相对误差")
print("-" * 60)
for x in large_values:
li_x = logint(x)
asymptotic = x / log(x)
error = abs(li_x - asymptotic) / li_x
print(f"{x}\t\t{float(li_x):.2f}\t\t{asymptotic:.2f}\t\t{error:.4f}")
#x li(x) x/ln(x) 相对误差
#------------------------------------------------------------
#100 30.13 21.71 0.2792
#1000 177.61 144.76 0.1849
#10000 1246.14 1085.74 0.1287
#100000 9629.81 8685.89 0.0980
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def logarithmic_integral(input_str):
"""
对标 MATLAB 的 logint 函数的实现，支持标量和矩阵输入
对数积分定义 (SymPy 实现):
li(x) = ∫₀^x 1/ln(t) dt (Cauchy主值积分)
分类归属:
属于特殊函数中的积分函数，与指数积分、三角函数积分同属超越积分类
参数:
input_str: 输入表达式字符串，可以是标量、矩阵或符号表达式
返回:
计算结果（标量/矩阵）或错误信息字符串
示例:
>>> logarithmic_integral("2")
1.04516378011749
>>> logarithmic_integral("[[1, 2], [3, 4]]")
Matrix([
[0, 1.04516378011749],
[2.16358859424961, 4.06403425769027]])
>>> logarithmic_integral("x")
li(x)
"""
try:
# 解析输入为 SymPy 表达式
expr = sp.sympify(input_str, evaluate=False)
return sp.li(expr).evalf()
except Exception as e:
return f"运行时错误: {str(e)}"
# ----------------------
# 示例测试代码
# ----------------------
if __name__ == "__main__":
# 标量计算示例
print(logarithmic_integral("2"))
# 1.04516378011749
print(logarithmic_integral("1.5"))
# 0.125064986315296
# 符号计算示例
x = sp.Symbol('x')
print(logarithmic_integral("x"))
# li(x)
矩阵对数
L = logm(A) 是 A 的主矩阵对数,即 expm(A) 的倒数.
输出L是每个特征值都具有严格位于π和π之间的虚部的唯一对数. 如果A是奇异矩阵或在负实轴上具有特征值,则未定义主对数.
A — 输入矩阵,方阵
示例1. 控制理论（线性系统分析）
连续时间线性系统: dx/dt = A*x
离散化后: x[k+1] = exp(A*Δt)*x[k]
矩阵对数用于从离散系统恢复连续系统矩阵
假设我们有一个离散状态转移矩阵
discrete_matrix = [[0.8187, 0.0861], [-0.1722, 0.7321]]
dt = 0.1 # 采样时间
print("离散状态转移矩阵:")
Phi_d = logm(discrete_matrix)
print(Phi_d)
#离散状态转移矩阵:
#[[-0.188273748658856, 0.110252305505317 + 1.38777878078145e-17*I],
[-0.220504611010633, -0.299166311455144]]
恢复连续系统矩阵: A = logm(Φ_d)/Δt
A_continuous = Phi_d / dt
print(f"\n恢复的连续系统矩阵 A = logm(Φ_d)/{dt}:")
print(A_continuous)
#恢复的连续系统矩阵 A = logm(Φ_d)/0.1:
#[[-1.88273748658856, 1.10252305505317 + 1.38777878078145e-16*I],
[-2.20504611010633, -2.99166311455144]]
示例2. 机器人学（姿态插值）
在机器人学中，矩阵对数用于 SO(3) 群上的姿态插值
旋转矩阵 R1 和 R2 之间的插值: R(t) = exp(t * logm(R2 * R1^T)) * R1
定义两个旋转矩阵
R1 = [[1, 0, 0], [0, 0.866, -0.5], [0, 0.5, 0.866]] # 绕x轴旋转30度
R2 = [[0.866, 0, 0.5], [0, 1, 0], [-0.5, 0, 0.866]] # 绕y轴旋转30度
print("初始旋转矩阵 R1:")
R1_mat = logm(R1)
print(R1_mat)
#初始旋转矩阵 R1:
#[[0, 0, 0],
[0, -2.20004840142047e-5, -0.523611477769969],
[0, 0.523611477769969, -2.20004840142047e-5]]
print("\n目标旋转矩阵 R2:")
R2_mat = logm(R2)
print(R2_mat)
#目标旋转矩阵 R2:
#[[-2.20004840142047e-5, 0, 0.523611477769969],
[0, 0, 0],
[-0.523611477769969, 0, -2.20004840142047e-5]]
计算相对旋转
R_rel = R2_mat * R1_mat.T
print(f"\n相对旋转矩阵 R2 * R1^T:")
print(R_rel)
#相对旋转矩阵 R2 * R1^T:
#[[0, -0.274168979652451, -1.15197059463323e-5],
[0, 0, 0],
[0, 1.15197059463323e-5, 4.84021296859278e-10]]
示例3. 量子力学（时间演化算子）
在量子力学中，时间演化算子 U = exp(-iHt/ℏ)
矩阵对数用于从时间演化算子恢复哈密顿量
假设的时间演化算子（简化情况）
U = [[0.8776, -0.4794], [0.4794, 0.8776]] # 近似旋转矩阵
print("时间演化算子 U:")
U_mat = logm(U)
print(U_mat)
#时间演化算子 U:
#[[3.05999063647575e-6, -0.499969227585355],
[0.499969227585355, 3.05999063647575e-6]]
恢复哈密顿量: H = iℏ * logm(U) / t
假设 ℏ = 1, t = 1
H_recovered = I * U_mat
print(f"\n恢复的哈密顿量 H = i * logm(U):")
print(H_recovered)
#恢复的哈密顿量 H = i * logm(U):
#[[3.05999063647575e-6*I, -0.499969227585355*I],
[0.499969227585355*I, 3.05999063647575e-6*I]]
示例4. 统计学（协方差矩阵处理）
在多变量正态分布中，协方差矩阵的几何均值涉及矩阵对数
协方差矩阵的黎曼均值: exp( (logm(Σ1) + logm(Σ2)) / 2 )
两个协方差矩阵
Sigma1 = [[4, 2], [2, 3]]
Sigma2 = [[3, 1], [1, 5]]
print("协方差矩阵 Σ1:")
log_S1 = logm(Sigma1)
print(log_S1)
#协方差矩阵 Σ1:
#[[1.20371282991029, 0.655968236281469],
[0.655968236281469, 0.875728711769551]]
print("\n协方差矩阵 Σ2:")
log_S2 = logm(Sigma2)
print(log_S2)
#协方差矩阵 Σ2:
#[[1.05825343611739, 0.261275228690240],
[0.261275228690240, 1.58080389349787]]
计算黎曼均值
log_mean = (log_S1 + log_S2) / 2
print(f"\n对数域的均值 (logm(Σ1) + logm(Σ2))/2:")
print(log_mean)
#对数域的均值 (logm(Σ1) + logm(Σ2))/2:
#[[1.13098313301384, 0.458621732485855],
[0.458621732485855, 1.22826630263371]]
示例5. 特殊矩阵类型
单位矩阵
identity = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
print("单位矩阵的对数:")
print(logm(identity))
#单位矩阵的对数:
#[[0, 0, 0],
[0, 0, 0],
[0, 0, 0]]
对角矩阵
diagonal = [[2, 0, 0], [0, 3, 0], [0, 0, 5]]
print("\n对角矩阵的对数:")
print(logm(diagonal))
#对角矩阵的对数:
#[[0.693147180559945, 0, 0],
[0, 1.09861228866811, 0],
[0, 0, 1.60943791243410]]
旋转矩阵
rotation = [[0.866, -0.5], [0.5, 0.866]] # 30度旋转
print("\n旋转矩阵的对数:")
print(logm(rotation))
#旋转矩阵的对数:
#[[-2.20004840142047e-5, -0.523611477769969],
[0.523611477769969, -2.20004840142047e-5]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def logarithmic_matrix(input_str):
"""
计算输入矩阵的矩阵对数，对标 MATLAB 的 logm 函数。
参数:
input_str: 输入的字符串，代表矩阵的表达式。
返回:
如果输入是有效的方阵，则返回矩阵的对数（计算结果为浮点数）；
否则返回错误信息。
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 尝试将输入转换为矩阵
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is not None:
# 检查矩阵是否为方阵
if A.is_square:
# 计算矩阵的对数并转换为浮点数
result = A.log().evalf()
else:
# 非方阵输入判定为错误
error = True
else:
# 无法转换为矩阵的输入判定为错误
error = True
# 根据是否有错误返回结果或错误信息
return result if not error else f"输入错误: {input_str}，输入必须是方阵。"
except Exception as e:
# 捕获其他异常并返回错误信息
return f"错误: {e}"
# 示范代码
# 示例 1: 有效的方阵输入
input_str1 = "[[1, 2], [3, 4]]"
result1 = logarithmic_matrix(input_str1)
print(f"输入: {input_str1}, 结果: {result1}")
# Matrix([[-0.350439813998554 + 2.39111795445172*I, 0.929351205704702 - 1.0937621702091*I],
# [1.39402680855705 - 1.64064325531366*I, 1.0435869945585 + 0.750474699138069*I]])
以2为底的对数和浮点数分解数
Y = log2(X)计算X的元素以2为底的对数，满足2^Y=X.
X — 输入数组标量,向量,矩阵
示例1. 计算机科学（数据结构和算法）
二叉树高度计算
nodes = [1, 3, 7, 15, 31, 63] # 完全二叉树的节点数
完全二叉树的高度计算:
for n in nodes:
height = logTwo(n + 1) # 高度 ≈ log₂(n+1)
print(f"节点数 {n}: 树高度 ≈ {float(height):.2f}")
#节点数 1: 树高度 ≈ 1.00
#节点数 3: 树高度 ≈ 2.00
#节点数 7: 树高度 ≈ 3.00
#节点数 15: 树高度 ≈ 4.00
#节点数 31: 树高度 ≈ 5.00
#节点数 63: 树高度 ≈ 6.00
算法复杂度分析
input_sizes = [8, 64, 512, 4096]
二分查找的比较次数:
for size in input_sizes:
comparisons = logTwo(size)
print(f"数组大小 {size}: 最多比较次数 = {float(comparisons):.1f}")
#数组大小 8: 最多比较次数 = 3.0
#数组大小 64: 最多比较次数 = 6.0
#数组大小 512: 最多比较次数 = 9.0
#数组大小 4096: 最多比较次数 = 12.0
示例2. 信息论（熵和编码）
计算信息熵（以比特为单位）
probabilities = [0.125, 0.25, 0.5, 0.75, 0.875]
print("单个事件的信息量 (I = -log₂(p)):")
for p in probabilities:
information = -logTwo(p)
print(f"p = {p}: I = {float(information):.3f} 比特")
#单个事件的信息量 (I = -log₂(p)):
#p = 0.125: I = 3.000 比特
#p = 0.25: I = 2.000 比特
#p = 0.5: I = 1.000 比特
#p = 0.75: I = 0.415 比特
#p = 0.875: I = 0.193 比特
均匀分布的信息熵
alphabet_sizes = [2, 4, 8, 16, 32] # 符号表大小
print("\n均匀分布的熵 (H = log₂(N)):")
for N in alphabet_sizes:
entropy = logTwo(N)
print(f"符号数 {N}: 熵 = {float(entropy):.3f} 比特/符号")
#均匀分布的熵 (H = log₂(N)):
#符号数 2: 熵 = 1.000 比特/符号
#符号数 4: 熵 = 2.000 比特/符号
#符号数 8: 熵 = 3.000 比特/符号
#符号数 16: 熵 = 4.000 比特/符号
#符号数 32: 熵 = 5.000 比特/符号
示例3. 信号处理（动态范围）
计算动态范围（分贝）
bit_depths = [8, 12, 16, 24, 32] # ADC/DAC 位深
数字系统的动态范围:
for bits in bit_depths:
# 动态范围 ≈ 6.02 × 位深 (dB)
dynamic_range = 6.02 * float(logTwo(2 ** bits))
print(f"{bits}-位系统: 动态范围 ≈ {dynamic_range:.1f} dB")
#8-位系统: 动态范围 ≈ 48.2 dB
#12-位系统: 动态范围 ≈ 72.2 dB
#16-位系统: 动态范围 ≈ 96.3 dB
#24-位系统: 动态范围 ≈ 144.5 dB
#32-位系统: 动态范围 ≈ 192.6 dB
采样率分析
nyquist_ratios = [2, 4, 8, 16]
过采样比与信噪比改善:
for ratio in nyquist_ratios:
snr_improvement = 10 * float(logTwo(ratio))
print(f"过采样比 {ratio}: SNR改善 ≈ {snr_improvement:.1f} dB")
#过采样比 2: SNR改善 ≈ 10.0 dB
#过采样比 4: SNR改善 ≈ 20.0 dB
#过采样比 8: SNR改善 ≈ 30.0 dB
#过采样比 16: SNR改善 ≈ 40.0 dB
示例4. 图像处理（像素深度）
不同颜色深度的颜色数量
color_depths = [1, 4, 8, 16, 24, 32]
颜色深度与可表示颜色数:
for depth in color_depths:
colors = 2 ** depth
print(f"{depth}-位: {colors:,} 种颜色")
#1-位: 2 种颜色
#4-位: 16 种颜色
#8-位: 256 种颜色
#16-位: 65,536 种颜色
#24-位: 16,777,216 种颜色
#32-位: 4,294,967,296 种颜色
图像文件大小估算
resolutions = [(640, 480), (1280, 720), (1920, 1080), (3840, 2160)]
bpp = 24 # 24位真彩色
不同分辨率未压缩图像大小:
for width, height in resolutions:
pixels = width * height
bits = pixels * bpp
bytes_size = bits / 8
megabytes = bytes_size / (1024 * 1024)
print(f"{width}×{height}: {megabytes:.1f} MB")
#640×480: 0.9 MB
#1280×720: 2.6 MB
#1920×1080: 5.9 MB
#3840×2160: 23.7 MB
示例5. 密码学（密钥空间）
不同密钥长度的安全性
key_lengths = [64, 128, 256, 512, 1024]
密钥长度与暴力破解难度:
for length in key_lengths:
combinations = 2 ** length
security_level = logTwo(combinations)
print(f"{length}-位密钥: 2^{length} = 10^{float(security_level) / logTwo('10'):.1f} 种可能")
#64-位密钥: 2^64 = 10^19.3+0.0j 种可能
#128-位密钥: 2^128 = 10^38.5+0.0j 种可能
#256-位密钥: 2^256 = 10^77.1+0.0j 种可能
#512-位密钥: 2^512 = 10^154.1+0.0j 种可能
#1024-位密钥: 2^1024 = 10^inf+nanj 种可能
生日攻击复杂度
hash_lengths = [64, 128, 160, 256, 512] # 哈希输出长度（位）
生日攻击复杂度 (≈ √(2^n)):
for n in hash_lengths:
complexity = 2 ** (n / 2)
print(f"{n}-位哈希: 生日攻击复杂度 ≈ 2^{n / 2} = {complexity:.1e}")
#64-位哈希: 生日攻击复杂度 ≈ 2^32.0 = 4.3e+09
#128-位哈希: 生日攻击复杂度 ≈ 2^64.0 = 1.8e+19
#160-位哈希: 生日攻击复杂度 ≈ 2^80.0 = 1.2e+24
#256-位哈希: 生日攻击复杂度 ≈ 2^128.0 = 3.4e+38
#512-位哈希: 生日攻击复杂度 ≈ 2^256.0 = 1.2e+77
示例6. 数值分析和浮点数
浮点数精度分析
float_types = [
("半精度", 16),
("单精度", 32),
("双精度", 64),
("四精度", 128)
]
浮点数格式的精度:
for name, bits in float_types:
# 尾数位数近似计算
if bits == 16:
mantissa_bits = 10
elif bits == 32:
mantissa_bits = 23
elif bits == 64:
mantissa_bits = 52
else:
mantissa_bits = 112
precision_decimal = mantissa_bits * float(logTwo(10)) / logTwo(2)
print(f"{name} ({bits}位): 约 {precision_decimal:.1f} 位十进制精度")
#半精度 (16位): 约 33.2+0.0j 位十进制精度
#单精度 (32位): 约 76.4+0.0j 位十进制精度
#双精度 (64位): 约 172.7+0.0j 位十进制精度
#四精度 (128位): 约 372.1+0.0j 位十进制精度
条件数计算
matrix_sizes = [2, 4, 8, 16]
矩阵条件数的对数尺度:
for size in matrix_sizes:
# 随机矩阵的条件数通常随维度增加
cond_estimate = size ** 2
log_cond = logTwo(cond_estimate)
print(f"{size}×{size} 矩阵: log₂(条件数) ≈ {float(log_cond):.1f}")
#2×2 矩阵: log₂(条件数) ≈ 2.0
#4×4 矩阵: log₂(条件数) ≈ 4.0
#8×8 矩阵: log₂(条件数) ≈ 6.0
#16×16 矩阵: log₂(条件数) ≈ 8.0
示例7. 矩阵运算（特征值分析）
正定矩阵
positive_definite = [[4, 1, 0], [1, 4, 1], [0, 1, 4]]
正定矩阵的 log₂:
result1 = logTwo(positive_definite)
print(result1)
#[[2, 0, -inf]
[0, 2, 0]
[-inf, 0, 2]]
对角矩阵（特别简单的情况）
diagonal = [[2, 0, 0], [0, 4, 0], [0, 0, 8]]
对角矩阵的 log₂:
result2 = logTwo(diagonal)
print(result2)
#[[1, -inf, -inf]
[-inf, 2, -inf]
[-inf, -inf, 3]]
验证：log₂(diag(2,4,8)) = diag(1,2,3)
expected = [[1, 0, 0], [0, 2, 0], [0, 0, 3]]
print(f"期望结果: {expected}")
#期望结果:
#[[1, 0, 0],
[0, 2, 0],
[0, 0, 3]]
示例8. 符号计算和数学性质
基本性质验证
expr1 = logTwo(x*y)
print(f"log₂(x*y) = {expr1}")
#log₂(x*y) = log(x*y)/log(2)
expr2 = logTwo(x/y)
#log₂(x/y) = log(x/y)/log(2)
expr3 = logTwo(x**y)
print(f"log₂(x^y) = {expr3}")
换底公式验证
print(f"\n换底公式: log₂(x) = {expr1}")
#换底公式: log₂(x) = log(x*y)/log(2)
示例9. 实际工程计算
数据存储计算
data_sizes_gb = [1, 10, 100, 1000]
数据存储的地址空间需求:
for size in data_sizes_gb:
bytes_needed = size * 1024 ** 3
address_bits = logTwo(bytes_needed)
print(f"{size} GB: 需要 {float(address_bits):.1f} 位地址总线")
#1 GB: 需要 30.0 位地址总线
#10 GB: 需要 33.3 位地址总线
#100 GB: 需要 36.6 位地址总线
#1000 GB: 需要 40.0 位地址总线
网络带宽需求
file_sizes_mb = [1, 10, 100, 500]
download_times = [1, 10, 60, 300] # 秒
下载时间与所需带宽:
for size in file_sizes_mb:
for time in download_times:
bandwidth_mbps = (size * 8) / time
print(f"{size} MB 在 {time} 秒内: 需要 {bandwidth_mbps:.1f} Mbps")
#1 MB 在 1 秒内: 需要 8.0 Mbps
#1 MB 在 10 秒内: 需要 0.8 Mbps
#1 MB 在 60 秒内: 需要 0.1 Mbps
#1 MB 在 300 秒内: 需要 0.0 Mbps
#10 MB 在 1 秒内: 需要 80.0 Mbps
#10 MB 在 10 秒内: 需要 8.0 Mbps
#10 MB 在 60 秒内: 需要 1.3 Mbps
#10 MB 在 300 秒内: 需要 0.3 Mbps
#100 MB 在 1 秒内: 需要 800.0 Mbps
#100 MB 在 10 秒内: 需要 80.0 Mbps
#100 MB 在 60 秒内: 需要 13.3 Mbps
#100 MB 在 300 秒内: 需要 2.7 Mbps
#500 MB 在 1 秒内: 需要 4000.0 Mbps
#500 MB 在 10 秒内: 需要 400.0 Mbps
#500 MB 在 60 秒内: 需要 66.7 Mbps
#500 MB 在 300 秒内: 需要 13.3 Mbps
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def logarithmic_natural_2(input_str):
"""
对标MATLAB的log2函数（基础对数部分）
计算以2为底的对数，支持标量/矩阵/符号输入
参数:
input_str: 输入的字符串表达式，可以是标量、矩阵或SymPy表达式
返回:
以2为底的对数值，支持以下类型返回：
- 标量数值：返回浮点数结果
- 矩阵输入：返回同型矩阵，逐元素计算结果
- 符号表达式：返回符号表达式
- 错误输入：返回错误描述
注意:
本函数仅实现对数计算部分，MATLAB的浮点数分解功能（[f,e] = log2(x)）需额外实现
示例:
logarithmic_natural_2('8') -> 3.0
logarithmic_natural_2('[[2,4],[8,16]]') -> [[1.0, 2.0], [3.0, 4.0]]
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理数值型标量
if expr.is_number:
z = complex(expr)
result = np.log2(z)
# 处理符号表达式
elif expr.free_symbols:
result = sp.log(expr, 2) # 保持符号表达式形式
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(logarithmic_natural_2('8'))
# (2.9999999999999996+0j)
print(logarithmic_natural_2('3'))
# (1.5849625007211563+0j)
# 示例3: 符号表达式
print("\n示例3 符号表达式:")
print(logarithmic_natural_2('x'))
# log(x)/log(2)
print(logarithmic_natural_2('2**y'))
# log(2**y)/log(2)
# 示例4: 边界值测试
print("\n示例4 边界值:")
print(logarithmic_natural_2('0.5'))
# (-1+0j)
print(logarithmic_natural_2('1'))
# 0j
对数正态参数估计
[PAt，pCI]=lognfit(x)也返回参数估计值的95%置信区间。
[PAt，pCI]=lognfit(x，alpha)指定置信区间的置信水平为100（1-alpha）%。
x —— 样本数据, 矢量
alpha——显著性水平, 0.05（默认）|范围（0,1）内的标量
示例1. 金融资产收益率分析
股票日收益率数据（假设为正，实际中需要处理负收益率）
stock_returns = [1.02, 1.05, 0.98, 1.08, 1.01, 1.03, 1.06, 0.99, 1.04, 1.07]
result = lognfit(stock_returns, 0.05)
params, ci = result
mu, sigma = params
收益率对数正态分布参数:
print(f"μ (位置参数) = {mu:.4f}")
print(f"σ (尺度参数) = {sigma:.4f}")
print(f"μ的95%置信区间: [{ci[0][0]:.4f}, {ci[0][1]:.4f}]")
print(f"σ的95%置信区间: [{ci[1][0]:.4f}, {ci[1][1]:.4f}]")
#μ (位置参数) = 0.0320
#σ (尺度参数) = 0.0307
#μ的95%置信区间: [0.0129, 0.0511]
#σ的95%置信区间: [0.0173, 0.0442]
计算几何均值（原始尺度下的均值）
geometric_mean = np.exp(mu + sigma ** 2 / 2)
print(f"几何均值: {geometric_mean:.4f}")
#几何均值: 1.0330
示例2. 环境科学（污染物浓度）
某地区PM2.5浓度数据 (μg/m³)
pm25_concentrations = [15.2, 22.8, 18.5, 35.1, 28.9, 42.3, 19.7, 31.5, 26.8, 38.2,
24.1, 33.7, 29.4, 45.6, 21.3]
result = lognfit(pm25_concentrations, 0.05)
params, ci = result
mu, sigma = params
PM2.5浓度对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 3.3169, σ = 0.3070
计算原始尺度的统计量
mean_original = np.exp(mu + sigma ** 2 / 2)
median_original = np.exp(mu)
mode_original = np.exp(mu - sigma ** 2)
print(f"原始尺度均值: {mean_original:.2f} μg/m³")
print(f"原始尺度中位数: {median_original:.2f} μg/m³")
print(f"原始尺度众数: {mode_original:.2f} μg/m³")
#原始尺度均值: 28.91 μg/m³
#原始尺度中位数: 27.58 μg/m³
#原始尺度众数: 25.10 μg/m³
示例3. 工程可靠性分析
机械部件故障时间（小时）
failure_times = [1250, 890, 2100, 1500, 3200, 1800, 950, 2700, 1400, 2300,
1600, 1900, 1100, 2500, 1700]
result = lognfit(failure_times, 0.1)
params, ci = result
mu, sigma = params
部件寿命对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 7.4258, σ = 0.3661
计算可靠性指标
t_1000 = 1000 # 1000小时可靠性
reliability_1000 = 1 - norm.cdf((np.log(t_1000) - mu) / sigma)
print(f"1000小时可靠性: {reliability_1000:.4f}")
#1000小时可靠性: 0.9215
计算中位寿命
median_life = np.exp(mu)
print(f"中位寿命: {median_life:.0f} 小时")
#中位寿命: 1679 小时
示例4. 医学研究（药物浓度）
患者血药浓度数据 (mg/L)
drug_concentrations = [2.1, 3.8, 1.5, 4.2, 2.9, 5.1, 3.2, 4.8, 2.7, 3.9,
4.1, 2.3, 3.5, 4.9, 3.1, 2.8, 4.5, 3.7, 2.6, 4.3]
result = lognfit(drug_concentrations, 0.05)
params, ci = result
mu, sigma = params
血药浓度对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 1.2085, σ = 0.3098
计算治疗窗口
therapeutic_min = 2.0 # 最低有效浓度
therapeutic_max = 5.0 # 最高安全浓度
prob_below_min = norm.cdf((np.log(therapeutic_min) - mu) / sigma)
prob_above_max = 1 - norm.cdf((np.log(therapeutic_max) - mu) / sigma)
prob_therapeutic = 1 - prob_below_min - prob_above_max
print(f"低于治疗窗口概率: {prob_below_min:.4f}")
print(f"在治疗窗口内概率: {prob_therapeutic:.4f}")
print(f"高于治疗窗口概率: {prob_above_max:.4f}")
#低于治疗窗口概率: 0.0481
#在治疗窗口内概率: 0.8541
#高于治疗窗口概率: 0.0978
示例5. 经济学（收入分布）
模拟收入数据（万元/年）
incomes = [8.5, 12.3, 15.8, 22.1, 18.9, 35.2, 28.7, 42.5, 55.1, 68.9,
25.3, 32.8, 48.6, 75.3, 92.1, 38.4, 52.7, 65.8, 82.3, 45.6]
result = lognfit(incomes, 0.05)
params, ci = result
mu, sigma = params
收入对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 3.5916, σ = 0.6384
计算收入分位数
percentiles = [0.1, 0.25, 0.5, 0.75, 0.9]
for p in percentiles:
income_p = np.exp(mu + sigma * norm.ppf(p))
print(f"{int(p * 100)}%分位点收入: {income_p:.1f} 万元")
#10%分位点收入: 16.0 万元
#25%分位点收入: 23.6 万元
#50%分位点收入: 36.3 万元
#75%分位点收入: 55.8 万元
#90%分位点收入: 82.3 万元
示例6. 质量控制（产品尺寸）
零件尺寸测量数据 (mm)
part_dimensions = [10.02, 10.05, 9.98, 10.08, 10.01, 10.03, 10.06, 9.99,
10.04, 10.07, 10.00, 10.09, 9.97, 10.10, 10.02]
result = lognfit(part_dimensions, 0.01)
params, ci = result
mu, sigma = params
零件尺寸对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 2.3060, σ = 0.0039
计算过程能力指数
spec_lower = 9.95 # 规格下限
spec_upper = 10.15 # 规格上限
对数正态分布的过程能力
cpk = min(
(np.log(spec_upper) - mu) / (3 * sigma),
(mu - np.log(spec_lower)) / (3 * sigma)
)
print(f"过程能力指数 Cpk: {cpk:.3f}")
if cpk >= 1.33:
print("过程能力充足")
elif cpk >= 1.0:
print("过程能力尚可")
else:
print("过程能力不足")
#过程能力指数 Cpk: 0.718
#过程能力不足
示例7. 网络分析（网页访问量）
网页日访问量数据
daily_visits = [150, 230, 180, 420, 290, 510, 320, 480, 270, 390,
410, 250, 350, 520, 310, 280, 450, 370, 260, 430]
result = lognfit(daily_visits, 0.05)
params, ci = result
mu, sigma = params
访问量对数正态分布参数:
print(f"μ = {mu:.4f}, σ = {sigma:.4f}")
#μ = 5.7880, σ = 0.3311
预测高流量日
high_traffic_threshold = np.exp(mu + 2 * sigma)
print(f"高流量日阈值: {high_traffic_threshold:.0f} 访问量")
#高流量日阈值: 633 访问量
计算高流量日概率
prob_high_traffic = 1 - norm.cdf(2) # P(Z > 2)
print(f"高流量日概率: {prob_high_traffic:.4f}")
#高流量日概率: 0.0228
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.stats import norm
def lognfit(input_data, alpha=0.05):
"""
对数正态分布参数估计（对标MATLAB lognfit）
参数:
input_data (array_like): 输入数据样本，必须为正数
alpha (float): 置信水平，默认0.05（返回95%置信区间）
返回:
tuple: (mu, sigma) 估计参数
tuple: (mu_ci, sigma_ci) 置信区间（可选）
示例:
>>> data = [1.2, 3.5, 2.8, 5.1, 4.7]
>>> mu, sigma = lognfit(data)
>>> print(f"mu: {mu:.4f}, sigma: {sigma:.4f}")
mu: 1.1931, sigma: 0.3050
>>> params, ci = lognfit(data, alpha=0.1)
>>> print(f"95% CI for mu: {ci[0]}")
"""
try:
# 输入数据验证
data = np.asarray(input_data, dtype=np.float64)
if np.any(data <= 0):
raise ValueError("所有数据必须为正数")
# 计算对数变换
log_data = np.log(data)
# 参数估计（最大似然估计）
mu = np.mean(log_data)
sigma = np.std(log_data, ddof=0) # 使用N而非N-1计算标准差
# 计算置信区间
n = len(data)
z_score = norm.ppf(1 - alpha / 2)
mu_se = sigma / np.sqrt(n)
mu_ci = (mu - z_score * mu_se, mu + z_score * mu_se)
sigma_se = sigma / np.sqrt(2 * n)
sigma_ci = (sigma - z_score * sigma_se, sigma + z_score * sigma_se)
return (mu, sigma), (mu_ci, sigma_ci)
except Exception as e:
return f"Error: {e}"
# 基本使用
data = [1.2, 3.5, 2.8, 5.1, 4.7]
params, ci = lognfit(data, 0.1)
print(f"μ估计值: {params[0]}")
# μ估计值: 1.1283013981833547
print(f"95%置信区间: [{ci[0][0]}, {ci[0][1]}]")
# 95%置信区间: [0.7465214703035561, 1.5100813260631534]
常用对数(以10为底)
Y = log10(X) 返回数组X中每个元素的常用对数.该函数同时接受实数和复数输入.对于X在区间(0, Inf)内的实数值,log10返回区间(-Inf ,Inf)内的实数值.
对于X的复数值和负实数值,log10函数返回复数值.
X — 输入数组标量,向量,矩阵
示例1. 声学（分贝计算）
声压级计算: Lp = 20 * log10(P/P0)，其中P0 = 20 μPa
sound_pressures = [0.002, 0.02, 0.2, 2.0] # 帕斯卡
声压级计算 (参考压力 20 μPa):
for pressure in sound_pressures:
db = 20 * log10(pressure / 0.00002)
print(f"声压 {pressure} Pa: {float(db):.1f} dB")
#声压 0.002 Pa: 40.0 dB
#声压 0.02 Pa: 60.0 dB
#声压 0.2 Pa: 80.0 dB
#声压 2.0 Pa: 100.0 dB
声音强度级
intensities = [1e-12, 1e-6, 1e-3, 1] # 瓦特/平方米
声音强度级 (参考强度 1e-12 W/m²):
for intensity in intensities:
db = 10 * log10(intensity / 1e-12)
print(f"强度 {intensity:.1e} W/m²: {float(db):.1f} dB")
#强度 1.0e-12 W/m²: 0.0 dB
#强度 1.0e-06 W/m²: 60.0 dB
#强度 1.0e-03 W/m²: 90.0 dB
#强度 1.0e+00 W/m²: 120.0 dB
示例2. 化学（pH值计算）
H+离子浓度
H_concentrations = [1e-2, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11] # mol/L
H+离子浓度与pH值:
for conc in H_concentrations:
pH = -log10(conc)
print(f"[H+] = {conc:.1e} mol/L: pH = {float(pH):.1f}")
#[H+] = 1.0e-02 mol/L: pH = 2.0
#[H+] = 1.0e-03 mol/L: pH = 3.0
#[H+] = 1.0e-05 mol/L: pH = 5.0
#[H+] = 1.0e-07 mol/L: pH = 7.0
#[H+] = 1.0e-09 mol/L: pH = 9.0
#[H+] = 1.0e-11 mol/L: pH = 11.0
OH-离子浓度与pOH
OH_concentrations = [1e-3, 1e-5, 1e-7, 1e-9]
OH-离子浓度与pOH:
for conc in OH_concentrations:
pOH = -log10(conc)
print(f"[OH-] = {conc:.1e} mol/L: pOH = {float(pOH):.1f}")
#[OH-] = 1.0e-03 mol/L: pOH = 3.0
#[OH-] = 1.0e-05 mol/L: pOH = 5.0
#[OH-] = 1.0e-07 mol/L: pOH = 7.0
#[OH-] = 1.0e-09 mol/L: pOH = 9.0
示例3. 地震学（里氏震级）
里氏震级公式: M = log10(A) - log10(A0)
其中A是地震波振幅，A0是标准振幅
wave_amplitudes = [0.001, 0.01, 0.1, 1, 10, 100] # 毫米
地震波振幅与里氏震级 (假设A0=0.001mm):
for amplitude in wave_amplitudes:
magnitude = log10(amplitude / 0.001)
print(f"振幅 {amplitude} mm: 震级 {float(magnitude):.1f}")
#振幅 0.001 mm: 震级 0.0
#振幅 0.01 mm: 震级 1.0
#振幅 0.1 mm: 震级 2.0
#振幅 1 mm: 震级 3.0
#振幅 10 mm: 震级 4.0
#振幅 100 mm: 震级 5.0
能量释放对比
magnitudes = [4, 5, 6, 7, 8]
不同震级的能量释放对比:
for mag in magnitudes:
# 能量E ∝ 10^(1.5M)
energy_ratio = 10 ** (1.5 * mag)
print(f"震级 {mag}: 能量是震级4的 {energy_ratio:.0f} 倍")
#震级 4: 能量是震级4的 1000000 倍
#震级 5: 能量是震级4的 31622777 倍
#震级 6: 能量是震级4的 1000000000 倍
#震级 7: 能量是震级4的 31622776602 倍
#震级 8: 能量是震级4的 1000000000000 倍
示例4. 天文学（星等系统）
视星等公式: m1 - m2 = -2.5 * log10(F1/F2)
其中F是光通量
brightness_ratios = [1, 2.512, 6.31, 15.85, 39.81, 100]
亮度比与星等差:
for ratio in brightness_ratios:
magnitude_diff = -2.5 * log10(ratio)
print(f"亮度比 {ratio}: 星等差 {float(magnitude_diff):.1f}")
#亮度比 1: 星等差 0.0
#亮度比 2.512: 星等差 -1.0
#亮度比 6.31: 星等差 -2.0
#亮度比 15.85: 星等差 -3.0
#亮度比 39.81: 星等差 -4.0
#亮度比 100: 星等差 -5.0
绝对星等与距离模数
distances_pc = [10, 100, 1000, 10000] # 秒差距
距离与距离模数:
for dist in distances_pc:
distance_modulus = 5 * log10(dist / 10)
print(f"距离 {dist} pc: 距离模数 {float(distance_modulus):.1f}")
#距离 10 pc: 距离模数 0.0
#距离 100 pc: 距离模数 5.0
#距离 1000 pc: 距离模数 10.0
#距离 10000 pc: 距离模数 15.0
示例5. 信号处理（信噪比）
信噪比计算: SNR = 10 * log10(Ps/Pn)
signal_powers = [1e-3, 1e-2, 1e-1, 1] # 瓦特
noise_power = 1e-6 # 瓦特
信号功率与信噪比 (噪声功率 1μW):
for power in signal_powers:
snr = 10 * log10(power / noise_power)
print(f"信号功率 {power:.1e} W: SNR = {float(snr):.1f} dB")
#信号功率 1.0e-03 W: SNR = 30.0 dB
#信号功率 1.0e-02 W: SNR = 40.0 dB
#信号功率 1.0e-01 W: SNR = 50.0 dB
#信号功率 1.0e+00 W: SNR = 60.0 dB
电压信噪比
signal_voltages = [0.1, 1, 10] # 伏特
noise_voltage = 0.01 # 伏特
电压信噪比:
for voltage in signal_voltages:
snr_voltage = 20 * log10(voltage / noise_voltage)
print(f"信号电压 {voltage} V: SNR = {float(snr_voltage):.1f} dB")
#信号电压 0.1 V: SNR = 20.0 dB
#信号电压 1 V: SNR = 40.0 dB
#信号电压 10 V: SNR = 60.0 dB
示例6. 计算机科学（数据量表示）
数据量单位转换
bytes_sizes = [1, 1024, 1048576, 1073741824, 1099511627776]
units = ["B", "KB", "MB", "GB", "TB"]
数据量对数表示:
for size, unit in zip(bytes_sizes, units):
log_size = log10(size)
print(f"{size} {unit}: log₁₀ = {float(log_size):.1f}")
#1 B: log₁₀ = 0.0
#1024 KB: log₁₀ = 3.0
#1048576 MB: log₁₀ = 6.0
#1073741824 GB: log₁₀ = 9.0
#1099511627776 TB: log₁₀ = 12.0
图像分辨率
resolutions = [(640, 480), (1280, 720), (1920, 1080), (3840, 2160)]
图像分辨率像素数:
for width, height in resolutions:
pixels = width * height
log_pixels = log10(pixels)
print(f"{width}×{height}: {pixels:,} 像素, log₁₀ = {float(log_pixels):.1f}")
#640×480: 307,200 像素, log₁₀ = 5.5
#1280×720: 921,600 像素, log₁₀ = 6.0
#1920×1080: 2,073,600 像素, log₁₀ = 6.3
#3840×2160: 8,294,400 像素, log₁₀ = 6.9
示例7. 生物学（感觉生理学）
韦伯-费希纳定律: S = k * log10(I)
其中S是感觉强度，I是刺激强度
stimulus_intensities = [1, 10, 100, 1000, 10000]
刺激强度与感觉强度 (韦伯-费希纳定律):
for intensity in stimulus_intensities:
sensation = log10(intensity)
print(f"刺激强度 {intensity}: 相对感觉强度 {float(sensation):.1f}")
#刺激强度 1: 相对感觉强度 0.0
#刺激强度 10: 相对感觉强度 1.0
#刺激强度 100: 相对感觉强度 2.0
#刺激强度 1000: 相对感觉强度 3.0
#刺激强度 10000: 相对感觉强度 4.0
声音频率感知
frequencies = [100, 1000, 10000] # Hz
频率感知 (近似对数关系):
for freq in frequencies:
perceived = log10(freq)
print(f"频率 {freq} Hz: 相对感知 {float(perceived):.1f}")
#频率 100 Hz: 相对感知 2.0
#频率 1000 Hz: 相对感知 3.0
#频率 10000 Hz: 相对感知 4.0
示例8. 经济学（数量级分析）
经济指标对数尺度
gdp_values = [1e9, 1e10, 1e11, 1e12, 1e13] # 美元
countries = ["小国", "中等国家", "大国", "主要经济体", "超级大国"]
GDP数量级分析:
for gdp, country in zip(gdp_values, countries):
log_gdp = log10(gdp)
print(f"{country}: GDP ${gdp:.0e}, log₁₀ = {float(log_gdp):.1f}")
#小国: GDP $1e+09, log₁₀ = 9.0
#中等国家: GDP $1e+10, log₁₀ = 10.0
#大国: GDP $1e+11, log₁₀ = 11.0
#主要经济体: GDP $1e+12, log₁₀ = 12.0
#超级大国: GDP $1e+13, log₁₀ = 13.0
人口数量级
populations = [1e5, 1e6, 1e7, 1e8, 1e9]
人口数量级:
for pop in populations:
log_pop = log10(pop)
print(f"人口 {pop:.0e}: log₁₀ = {float(log_pop):.1f}")
#人口 1e+05: log₁₀ = 5.0
#人口 1e+06: log₁₀ = 6.0
#人口 1e+07: log₁₀ = 7.0
#人口 1e+08: log₁₀ = 8.0
#人口 1e+09: log₁₀ = 9.0
示例9. 矩阵运算（科学计算）
物理量矩阵
concentration_matrix = [[1e-3, 1e-2, 1e-1], [1, 10, 100], [1000, 10000, 100000]]
浓度矩阵的 log₁₀:
result1 = log10(concentration_matrix)
print(result1)
#[[-3.0, -2.0, -1.0],
[0, 1.0, 2.0],
[3.0, 4.0, 5.0]]
信号强度矩阵
signal_matrix = [[0.001, 0.01, 0.1], [1, 10, 100], [1000, 10000, 100000]]
信号强度矩阵的 log₁₀:
result2 = log10(signal_matrix)
print(result2)
#[[-3.0, -2.0, -1.0],
[0, 1.0, 2.0],
[3.0, 4.0, 5.0]]
验证对角矩阵
diagonal_matrix = [[0.01, 0, 0], [0, 1, 0], [0, 0, 100]]
对角矩阵的 log₁₀:
result3 = log10(diagonal_matrix)
print(result3)
#[[-2.0, -oo, -oo],
[-oo, 0, -oo],
[-oo, -oo, 2.0]]
示例10. 符号计算和数学推导
基本性质
expr1 = log10(x*y)
print(f"log₁₀(x*y) = {expr1}")
#log₁₀(x*y) = log(x*y)/log(10)
expr2 = log10(x/y)
print(f"log₁₀(x/y) = {expr2}")
#log₁₀(x/y) = log(x/y)/log(10)
expr3 = log10(x**y)
print(f"log₁₀(x^y) = {expr3}")
#log₁₀(x^y) = log(x**y)/log(10)
换底公式应用
print(f"\n换底公式: log₁₀(x) = {expr1}")
#换底公式: log₁₀(x) = log(x*y)/log(10)
微分计算
derivative = diff(expr1, x)
print(f"d(log₁₀(x))/dx = {derivative}")
#d(log₁₀(x))/dx = 0
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def logarithmic_natural_10(input_str):
"""
对标MATLAB的log10函数，计算以10为底的对数
支持标量、矩阵和符号表达式输入
参数:
input_str: 输入的字符串表达式，可以是：
- 数值（如'100'）
- 矩阵（如'[[1, 10], [100, 1000]]'）
- 符号表达式（如'x'）
返回:
计算结果，类型与输入对应：
- 数值输入：返回浮点数结果
- 矩阵输入：返回同型矩阵，逐元素计算结果
- 符号表达式：返回符号表达式
- 错误输入：返回错误描述字符串
示例:
>>> logarithmic_natural_10('100')
2.00000000000000
>>> logarithmic_natural_10('[[1, 10], [100, 1000]]')
Matrix([
[0, 1],
[2, 3]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 处理数值型标量
if expr.is_number:
z = complex(expr)
result = np.log10(z)
# 处理符号表达式
elif expr.free_symbols:
result = sp.log(expr, 10) # 保持符号表达式形式
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(logarithmic_natural_10('1000'))
# (3+0j)
print(logarithmic_natural_10('3'))
# (0.4771212547196625+0j)
# 示例3: 符号表达式
print("\n示例3 符号表达式:")
print(logarithmic_natural_10('x'))
# log(x)/log(10)
print(logarithmic_natural_10('10**y'))
# log(10**y)/log(10)
# 示例4: 边界值测试
print("\n示例4 边界值:")
print(logarithmic_natural_10('0.1'))
# (-0.9999999999999999+0j)
print(logarithmic_natural_10('1'))
# 0j
对数拟合
p = logfit(x,y,n) 返回次数为n的对数 p(x) 的系数.
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是对数拟合的次数, 正整数标量.
Y = a*log(x)+b
Y = a*log10(x)+b
Y = a*log2(x)+b
示例1. 生物学：细菌生长模型
细菌在有限营养条件下的生长（对数增长）
bacteria_result = logfit([1,2,3,4,5,6], [50,85,120,150,175,195])
print(f"拟合公式: {bacteria_result['expression']}")
print(f"R² = {bacteria_result['goodness_of_fit']['r_squared']:.4f}")
print(f"斜率(生长率): {bacteria_result['coefficients']['slope']:.2f}")
#拟合公式: y = 81.9501*ln(x) + 39.3050
#R² = 0.9719
#斜率(生长率): 81.95
#实际意义：描述细菌在有限资源下的减速增长
示例2. 经济学：学习曲线效应
生产过程中的学习曲线（单位成本随产量增加而降低）
learning_curve = logfit([100,200,500,1000,2000], [45,38,30,25,22], 10)
print(f"拟合公式: {learning_curve['expression']}")
print(f"R² = {learning_curve['goodness_of_fit']['r_squared']:.4f}")
print(f"每产量翻倍，成本下降约: {learning_curve['coefficients']['slope'] * np.log10(2):.2f}单位")
#拟合公式: y = -17.9424*log10(x) + 79.7306
#R² = 0.9824
#每产量翻倍，成本下降约: -5.40单位
示例3. 心理学：韦伯-费希纳定律
感觉强度与刺激强度的对数关系
weber_fechner = logfit([10,20,50,100,200], [2.3,3.0,3.9,4.6,5.3])
print(f"拟合公式: {weber_fechner['expression']}")
print(f"R² = {weber_fechner['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = 0.9996*ln(x) + -0.0013
#R² = 1.0000
#实际意义：感觉强度与刺激强度的对数成正比
示例4. 声学：分贝与声压级
声音强度与分贝值的对数关系
sound_level = logfit([0.00002,0.002,0.02,0.2], [0,40,60,80], 10)
print(f"拟合公式: {sound_level['expression']}")
print(f"R² = {sound_level['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = 20.0000*log10(x) + 93.9794
#R² = 1.0000
#验证：20*log10(P/P0) 的关系
示例5. 计算机科学：算法复杂度
二分查找的时间复杂度 O(log n)
binary_search = logfit([100,1000,10000,100000], [7,10,14,17], 2)
print(f"拟合公式: {binary_search['expression']}")
print(f"R² = {binary_search['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = 1.0235*log2(x) + 0.1000
#R² = 0.9966
#实际意义：二分查找的比较次数与数据规模的对数成正比
示例6. 地球科学：里氏震级
地震能量与震级的对数关系
earthquake = logfit([10000,100000,1000000], [4,5,6], 10)
print(f"拟合公式: {earthquake['expression']}")
print(f"R² = {earthquake['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = 1.0000*log10(x) + -0.0000
#R² = 1.0000
#实际意义：震级增加1级，能量释放增加约32倍
示例7. 金融：复利计算
投资回报的对数增长
investment = logfit([1,5,10,20], [1100,1610,2590,6720])
print(f"拟合公式: {investment['expression']}")
print(f"R² = {investment['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = 1605.8087*ln(x) + 231.8666
#R² = 0.6513
#实际意义：长期投资的价值增长近似对数关系
示例8. 化学：pH值计算
氢离子浓度与pH值的对数关系
ph_calculation = logfit([0.1,0.01,0.001,0.0001], [1,2,3,4], 10)
print(f"拟合公式: {ph_calculation['expression']}")
print(f"R² = {ph_calculation['goodness_of_fit']['r_squared']:.4f}")
#拟合公式: y = -1.0000*log10(x) + -0.0000
#R² = 1.0000
#实际意义：pH = -log10[H⁺]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.optimize import curve_fit
def log_fit_nonlinear(input_str):
"""
实现MATLAB风格的logfit对数拟合函数
参数:
input_str: 输入字符串，格式为 (x_data, y_data) 或 (x_data, y_data, log_base)
log_base可选值为1(自然对数)/2/10，默认为1
返回:
字典包含:
- expression: 拟合公式字符串
- coefficients: 系数字典
- r_squared: 拟合优度R²
- params_cov: 参数协方差矩阵
示例:
>>> log_fit_nonlinear("([1,2,3], [5,8,10])")
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str)
if not isinstance(expr, tuple) or len(expr) not in (2, 3):
return "错误：输入格式应为 (x, y) 或 (x, y, log_base)"
# 提取数据
x_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
y_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
if x_sym is None or y_sym is None:
return "错误：x或y数据格式不正确"
# 转换为numpy数组
x_data = np.array(x_sym).astype(float).flatten()
y_data = np.array(y_sym).astype(float).flatten()
# 验证数据有效性
if len(x_data) != len(y_data):
return "错误：x和y数据长度不一致"
if np.any(x_data <= 0):
return "错误：x数据必须全为正数"
# 获取对数基底参数
log_base = expr[2] if len(expr) == 3 else 1
valid_bases = {1: np.log, 2: np.log2, 10: np.log10}
if log_base not in valid_bases:
return "错误：log_base必须为1/2/10"
# 定义拟合函数
def fit_func(x, a, b):
return a * valid_bases[log_base](x) + b
# 执行曲线拟合
initial_guess = [1.0, 1.0]
params, params_cov = curve_fit(
fit_func,
x_data,
y_data,
p0=initial_guess,
maxfev=10000
)
# 计算拟合指标
residuals = y_data - fit_func(x_data, *params)
ss_res = np.sum(residuals ** 2)
ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)
r_squared = 1 - (ss_res / ss_tot)
# 构造结果
base_symbol = {
1: "ln",
2: "log2",
10: "log10"
}[log_base]
return {
"expression": f"y = {params[0]:.4f}*{base_symbol}(x) + {params[1]:.4f}",
"coefficients": {
"slope": params[0],
"intercept": params[1]
},
"goodness_of_fit": {
"r_squared": r_squared,
"mse": ss_res / len(x_data)
},
"covariance_matrix": params_cov
}
except Exception as e:
return f"拟合错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：自然对数拟合
print("示例1：自然对数拟合")
result1 = log_fit_nonlinear("([2,3,4,5], [3.1,4.2,5.0,5.8])")
print("拟合结果:", result1["expression"])
# y = 2.9090*ln(x) + 1.0433
print("R平方值:", result1["goodness_of_fit"]["r_squared"])
# 0.9963528403421813
# 示例2：以10为底的对数拟合
print("\n示例2：log10拟合")
result2 = log_fit_nonlinear("([10,100,1000], [5,7,9], 10)")
print("拟合结果:", result2["expression"])
# y = 2.0000*log10(x) + 3.0000
print("协方差矩阵:\n", result2["covariance_matrix"])
# [[ 0. -0.]
# [-0. 0.]]
逻辑拟合
p = logisticfit(x,y,n) 返回次数为n的对数 p(x) 的系数.
x是查询点，指定为一个向量.
y是查询点位置的拟合值,指定为向量.
n是对数拟合的次数, 正整数标量.
logistic1: Y = a/(1+exp(-b*(x-c)))
logistic4: Y = d+(a-d)/(1+(x/c)^b)
1. 生物学：种群增长模型
细菌种群在有限资源下的增长（S型曲线）
population_data = "[0,1,2,3,4,5,6,7,8,9,10], [10,15,25,45,80,140,220,290,330,350,360], 1"
result = logisticfit(population_data)
print(f"拟合公式: {result}")
#拟合公式: 169.5455/(3.9949742803431e-1136*exp(27.3545*x) + 1)
#实际意义：描述种群从指数增长到环境承载力的转变
2. 医学：药物剂量反应
药物剂量与疗效关系（EC50曲线）
dose_response = "[0.1,0.5,1,5,10,50,100], [5,15,30,65,85,95,98], 4"
result = logisticfit(dose_response)
print(f"拟合公式: {result}")
#拟合公式: (98.0*x**41.9154 + 3.56307436854162e+71)/(x**41.9154 + 8.90768592135406e+69)
#实际意义：确定半数有效剂量(EC50)和最大疗效
3. 市场营销：产品 adoption 曲线
新产品市场渗透率（创新扩散理论）
adoption_curve = "[0,3,6,9,12,15,18,21,24], [2,5,12,28,55,78,88,93,95], 1"
result = logisticfit(adoption_curve)
print(f"拟合公式: {result}")
#拟合公式: 50.6667/(5.62549012907568e-16903*exp(163.9124*x) + 1)
#实际意义：预测产品生命周期和最大市场占有率
4. 心理学：学习曲线
技能学习过程中的进步速度
learning_data = "[1,2,3,4,5,6,7,8,9,10], [20,35,55,75,88,94,97,98,99,99], 1"
result = logisticfit(learning_data)
print(f"拟合公式: {result}")
#拟合公式: 76.0/(6.33199250994413e-4345*exp(88.3378*x) + 1)
#实际意义：描述从初学者到熟练者的学习过程
5. 经济学：技术扩散模型
新技术在产业中的采纳率
tech_adoption = "[1980,1985,1990,1995,2000,2005,2010], [1,5,15,40,75,90,95], 1"
result = logisticfit(tech_adoption)
print(f"拟合公式: {result}")
#拟合公式: 45.8571*exp(1.0*x)/(exp(1.0*x) + 2.71828182845905)
#实际意义：预测新技术普及速度和最终市场占有率
6. 生态学：物种竞争模型
两种物种竞争同一生态位的种群动态
competition_data = "[0,1,2,3,4,5,6,7], [100,180,300,450,580,650,690,700], 1"
result = logisticfit(competition_data)
print(f"拟合公式: {result}")
#拟合公式: 456.25*exp(263.0155*x)/(exp(263.0155*x) + 6.75476415563823e-19354)
#实际意义：预测环境承载力和种群稳定状态
7. 金融：市场饱和度模型
金融产品市场渗透的S曲线
finance_data = "[0,2,4,6,8,10,12], [5,15,40,75,90,95,96], 4"
result = logisticfit(finance_data)
print(f"拟合公式: {result}")
#拟合公式: (5.457*x**48.7492 + 7.53882928513792e+82)/(x**48.7492 + 1.26855239482975e+81)
#实际意义：评估金融产品的增长潜力和市场天花板
8. 农业：肥料效应曲线
肥料用量与作物产量的关系
fertilizer_data = "[0,50,100,150,200,250,300], [2000,3500,4800,5200,5300,5250,5200], 4"
result = logisticfit(fertilizer_data)
print(f"拟合公式: {result}")
#拟合公式: (13.4215000000004*x**540.6564 + 8.23214765753031e+1800)/(x**540.6564 + 1.84400108118759e+1797)
#实际意义：确定最佳肥料用量和最大产量潜力
9. 流行病学：疾病传播模型
传染病在人群中的传播
epidemic_data = "[0,5,10,15,20,25,30], [10,50,200,800,2500,4500,4800], 1"
result = logisticfit(epidemic_data)
print(f"拟合公式: {result}")
#拟合公式: 5002.6868*exp(0.367*x)/(exp(0.367*x) + 1428.21086990731)
#实际意义：预测疫情高峰时间和最终感染规模
10. 工程学：材料疲劳曲线
材料在不同应力水平下的寿命
fatigue_data = "[100,200,300,400,500,600], [1000000,500000,100000,50000,10000,5000], 4"
result = logisticfit(fatigue_data)
print(f"拟合公式: {result}")
#拟合公式: (1.89181594691429e-23379*x**9021.6449 + 533333.3342)/(8.73145828369257e-23384*x**9021.6449 + 1)
#实际意义：预测材料在不同应力下的使用寿命
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.optimize import curve_fit
def logistic_fit_nonlinear(input_str):
"""
使用非线性最小二乘法拟合逻辑斯蒂曲线
参数:
input_str: 输入字符串，格式为包含x数据、y数据和可选参数n的元组，例如:
"([[x1,x2,...], [y1,y2,...]], n)" 或 "([x1,x2,...], [y1,y2,...], n)"
返回:
拟合后的SymPy表达式或错误信息字符串
"""
try:
expr = sp.sympify(input_str) # 解析输入字符串为SymPy表达式
error = False
maxfev = 10000 # 最大函数评估次数
# 定义三参数逻辑斯蒂函数 (n=1)
def logistic_function1(x, a, b, c):
return a / (1 + np.exp(-b * (x - c)))
# 定义四参数S形函数 (n=4)
def logistic_function2(x, d, a, c, b):
return d + (a - d) / (1 + (x / c) ** b)
# 解析输入数据
if isinstance(expr, tuple):
# 提取x,y数据及模型参数n
if len(expr) > 2:
x, y, n = expr[0], expr[1], int(expr[2])
else:
x, y = expr[0], expr[1]
n = 1 # 默认使用三参数模型
# 将输入转换为NumPy数组
X = sp.Matrix(x) if isinstance(x, list) else None
Y = sp.Matrix(y) if isinstance(y, list) else None
if X is not None and Y is not None:
x_data = np.array(np.ravel(X), dtype=float)
y_data = np.array(np.ravel(Y), dtype=float)
else:
error = True
else:
error = True
# 根据模型类型进行拟合
if not error:
if n == 1:
# 三参数模型拟合
initial_guess = [1, 1, 1]
params, _ = curve_fit(logistic_function1, x_data, y_data,
p0=initial_guess, maxfev=maxfev)
a, b, c = params
expression = f"{a:.4f} / (1 + exp(-{b:.4f} * (x - {c:.4f})))"
elif n == 4:
# 四参数模型拟合
initial_guess = [0, 1, 1, 1]
params, _ = curve_fit(logistic_function2, x_data, y_data,
p0=initial_guess, maxfev=maxfev)
d, a, c, b = params
expression = f"{d:.4f} + ({a:.4f} - {d:.4f}) / (1 + (x / {c:.4f})^{b:.4f})"
else:
error = True
return sp.simplify(expression) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例1: 三参数模型拟合
input_str1 = "([[1,0],[1,0]],[7,8,14,18],1)"
result1 = logistic_fit_nonlinear(input_str1)
print("三参数模型拟合结果:", result1)
# 11.75*exp(19.2551*x)/(exp(19.2551*x) + 1.08537673309963e-139)
# 示例2: 四参数模型拟合
input_str2 = "([[1,0],[1,0]],[7,8,14,18],4)"
result2 = logistic_fit_nonlinear(input_str2)
print("四参数模型拟合结果:", result2)
# 1.0*(8.0*x**1.0 + 13.0)/(x**1.0 + 1)
存在已知协方差情况下的最小二乘解
x = lscov(A,b) 返回最小化误差平方和r'*r的最小二乘解.
x = lscov(A,b,w) 返回最小化r'*diag(w)*r 的加权最小二乘解,其中 r=b-A*x.
[x,stdx] = lscov(___) 还返回 x 的估计标准误差。标准误差是 x 的标准差的估计值。您可以使用上述语法中的任何输入参量组合。
[x,stdx,mse] = lscov(___) 还返回均方误差。均方误差与函数最小化的值成正比。
A, b — 操作数,向量,矩阵.
w — 相对权重,非负实数列向量.
x — 解, 向量 | 矩阵
stdx — 估计的标准误差, 向量
mse — 均方误差, 标量
示例1. 物理学：仪器测量精度不同
弹簧伸长与力的关系，使用不同精度的测量仪器
A = [[1, 1], [1, 2], [1, 3], [1, 4], [1, 5]]
b = [1.02, 1.98, 3.05, 4.01, 4.97]
不同仪器的测量方差
V = [[0.0001, 0, 0, 0, 0],
[0, 0.0004, 0, 0, 0],
[0, 0, 0.0001, 0, 0],
[0, 0, 0, 0.0009, 0],
[0, 0, 0, 0, 0.0001]]
x, stdx, mse = lscov(A, b, V)
print(f"胡克定律系数: {x}")
print(f"弹簧系数: {x[1]:.4f} ± {stdx[1]:.4f}")
#胡克定律系数: [0.04346192 0.98894164]
#弹簧系数: 0.9889 ± 0.0035
#实际意义：考虑不同仪器的测量精度差异
示例2. 气象学：空间插值问题
不同地理位置的气温测量，考虑空间相关性
A = [[1, 0, 0], [1, 1, 0], [1, 0, 1], [1, 1, 1], [1, 2, 1]]
b = [15.2, 16.1, 14.8, 16.5, 17.2]
空间相关的协方差矩阵
V = [[1.0, 0.8, 0.6, 0.4, 0.2],
[0.8, 1.0, 0.7, 0.5, 0.3],
[0.6, 0.7, 1.0, 0.6, 0.4],
[0.4, 0.5, 0.6, 1.0, 0.7],
[0.2, 0.3, 0.4, 0.7, 1.0]]
x, stdx, mse = lscov(A, b, V)
print(f"温度场模型: T = {x[0]:.2f} + {x[1]:.2f}*x + {x[2]:.2f}*y")
#温度场模型: T = 15.20 + 1.04*x + -0.14*y
#实际意义：考虑测量点之间的空间相关性
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def least_squares_covariance(input_str):
"""
对标MATLAB的lscov函数，计算存在协方差情况下的最小二乘解。
参数:
input_str: 输入的表达式字符串，格式为(A, b)或(A, b, V)，其中：
- A 是设计矩阵
- b 是观测向量
- V 是协方差矩阵（二维方阵）或方差向量（一维）
返回:
最小二乘解，或错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
x, stdx, mse = None, None, None
if isinstance(expr, tuple):
# 普通最小二乘 (A, b)
if len(expr) == 2:
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
b_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
if A_sym is not None and b_sym is not None:
A = np.array(A_sym.tolist(), dtype=float)
b = np.array(b_sym.tolist(), dtype=float).flatten()
x, residuals, _, _ = np.linalg.lstsq(A, b, rcond=None)
mse = np.sum(residuals) / (A.shape[0] - A.shape[1])
cov_matrix = mse * np.linalg.inv(A.T @ A)
stdx = np.sqrt(np.diag(cov_matrix))
# 加权最小二乘 (A, b, V)
elif len(expr) == 3:
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
b_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
V_sym = sp.Matrix(expr[2]) if isinstance(expr[2], list) else None
if all([m is not None for m in [A_sym, b_sym, V_sym]]):
A = np.array(A_sym.tolist(), dtype=float)
b = np.array(b_sym.tolist(), dtype=float).flatten()
V = np.array(V_sym.tolist(), dtype=float)
# 处理协方差矩阵
if V.ndim == 2:
try:
L = np.linalg.cholesky(V)
except np.linalg.LinAlgError:
return "错误: 协方差矩阵不正定"
A_prime = np.linalg.solve(L, A)
b_prime = np.linalg.solve(L, b)
else:
v = V.flatten()
if np.any(v <= 0):
return "错误: 方差必须为正数"
W_sqrt = np.diag(1.0 / np.sqrt(v))
A_prime = W_sqrt @ A
b_prime = W_sqrt @ b
# 求解加权最小二乘
x, residuals, _, _ = np.linalg.lstsq(A_prime, b_prime, rcond=None)
mse = np.sum(residuals) / (A_prime.shape[0] - A_prime.shape[1])
cov_matrix = np.linalg.inv(A_prime.T @ A_prime)
stdx = np.sqrt(np.diag(cov_matrix))
else:
return "输入参数数量错误"
else:
return "输入格式错误"
return x, stdx, mse
except Exception as e:
return f"错误: {e}"
# 示范代码
if __name__ == "__main__":
# 示例1: 普通最小二乘
A = [[1, 2], [3, 4], [5, 6]]
b = [1, 2, 3]
print("示例1 解:", least_squares_covariance(str((A, b))))
# (array([-5.97106181e-17, 5.00000000e-01]), array([8.47946841e-17, 6.70360838e-17]), 3.0814879110195774e-33)
# 示例2: 方差向量加权
V = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] # 对角线方差
print("示例2 解:", least_squares_covariance(str((A, b, V))))
# (array([-5.97106181e-17, 5.00000000e-01]), array([1.52752523, 1.20761473]), 3.0814879110195774e-33)
A = [[1, 0.2, 0.1], [1, 0.5, 0.3], [1, 0.6, 0.4], [1, 0.8, 0.9], [1, 1, 1.1], [1, 1.1, 1.4]]
b = [0.17, 0.26, 0.28, 0.23, 0.27, 0.24]
print("示例3 解:", least_squares_covariance(str((A, b))))
# (array([ 0.10184569, 0.48439898, -0.28465473]), array([0.00584332, 0.02060813, 0.01352573]), 1.277351520318261e-05)
斜坡响应图
lsimplot(sys) 它描绘了一个动态系统（由传递函数、状态空间模型等描述）在特定输入信号——单位斜坡输入 r(t) = t * u(t) （其中 u(t) 是单位阶跃函数）——作用下的输出响应随时间变化的曲线。
汽车巡航控制系统（加速阶段）
描述车辆在加速指令下，实际车速与期望车速的动态关系。
lsimplot(2,s^2+3s+2)
温度控制系统（缓慢升温）
描述烤箱在缓慢升温指令下，实际温度与期望温度的动态关系。
lsimplot(s+1,s^2+2s+5)
卫星姿态控制系统（缓慢转动）
描述卫星在缓慢转动指令下，实际角度与期望角度的动态关系。
lsimplot(4,s^2+2s+4)
电机位置控制系统（匀速运动）
描述电机在匀速位置指令下，实际位置与期望位置的动态关系。
lsimplot(3,s^2+4s+3)
伺服位置控制系统（机械臂）
斜坡输入： 目标位置匀速移动（如机械臂以 10 cm/s 的速度直线运动）。
响应特征： 输出滞后于输入形成稳态误差（如持续落后 0.5 cm），反映系统增益不足或摩擦影响。
实际意义：量化定位精度，指导是否需增加积分控制或提高电机扭矩。
lsimplot([0.8], [1, 0.2, 0.05])
温度爬升控制（工业熔炉）
斜坡输入： 设定温度以 5°C/min 匀速上升（如从 100°C 到 500°C）。
响应特征：初始阶段温度滞后，后期可能因热惯性超调（如实际温度超过设定值 20°C）。
实际意义： 评估热系统响应延迟，优化加热功率分配策略。
lsimplot([1.2, 0.1], [1, -0.7, 0.15])
无人机高度跟踪（定速爬升）
斜坡输入： 目标高度以 2 m/s 速度爬升（如起飞阶段）。
响应特征： 受风扰影响，高度波动围绕斜坡上下振荡（振幅 ±0.3 m）。
实际意义： 分析抗干扰能力，决定是否需增强气压计反馈或调整控制增益。
lsimplot([0.5], [1, -0.3, 0.1, 0.02])
数控机床进给系统
斜坡输入： 工作台以恒定速度 20 mm/s 移动（切削加工）。
响应特征： 稳态时存在固定滞后（如 0.2 mm），但无振荡。
实际意义： 揭示导轨摩擦导致的跟踪误差，需校准补偿或润滑维护。
lsimplot([1.0, 0.4], [1, -0.5])
经济匀速增长模型
斜坡输入： 预期经济以 3%/年 匀速增长。
响应特征： 初期增长滞后（政策时滞），后期可能因投资过热超过预期。
实际意义： 量化经济刺激政策的生效时间和过冲风险。
lsimplot([0.7], [1, -0.6, 0.08])
叠加斜坡响应图
lsimplot([sys1],[sys2]...[sysN])是在同一坐标系中绘制多个系统对同一斜坡输入的响应曲线，或同一系统对不同斜坡输入的响应曲线。
一阶系统对比
比较正向和反向温度控制系统的升温响应
lsimplot([1,s+2],[-1,s+2])
欠阻尼二阶系统
比较不同控制方向的卫星姿态调整
lsimplot([s+1,s^2+s+1],[-s-1,s^2+s+1])
参数变化对比
刚度增强后的机械系统正反向运动
lsimplot([s+1,s^2+s+1.5],[-s-1,s^2+s+1.5])
线性方程的最小范数最小二乘解
X = lsqminnorm(A,B) 返回可以求解线性方程 AX = B 并使 norm(A*X-B) 的值最小化的数组 X.
如果此问题有多个解,则 lsqminnorm 返回能使 norm(X) 最小化的解. 如果 B 有多列,则前面的语句分别适用于X和B的每列.
A — 系数矩阵,矩阵.
B — 输入数组,向量,矩阵.
示例1. 信号处理：欠定系统信号重建
从少量测量值重建原始信号（压缩感知）
A = [[1, 0, 1, 0, 1], # 测量矩阵
[0, 1, 0, 1, 0],
[1, 1, 0, 0, 1]]
b = [3, 2, 4] # 测量值
result = lsqminnorm(A, b)
print(f"重建信号: {result}")
print(f"信号能量: {np.linalg.norm(result):.4f}")
#重建信号: [[1.28571429]
[1.42857143]
[0.42857143]
[0.57142857]
[1.28571429]]
#信号能量: 2.4202
验证重建精度
reconstruction_error = np.linalg.norm(np.array(A) @ result - np.array(b))
print(f"重建误差: {reconstruction_error:.6f}")
#重建误差: 3.464102
#实际意义：在压缩感知中，从少量测量值恢复稀疏信号
示例2. 控制工程：冗余执行器控制分配
飞行器有多个控制面（冗余执行器），需要分配控制力矩
控制效率矩阵
A = [[1.0, 0.5, 0.3], # 副翼对滚转力矩的贡献
[0.2, 1.0, 0.6], # 升降舵对俯仰力矩的贡献
[0.1, 0.4, 1.0]] # 方向舵对偏航力矩的贡献
期望力矩
b = [2.0, 1.5, 0.8] # [滚转, 俯仰, 偏航]力矩
result = lsqminnorm(A, b)
print(f"控制面偏转量: {result}")
print(f"控制能量: {np.linalg.norm(result):.4f}")
#控制面偏转量:
#[[1.38888889]
[1.08625731]
[0.22660819]]
#控制能量: 1.7777
#实际意义：最小化控制面的偏转量，减少磨损和能耗
示例3. 机器人学：冗余机械臂逆运动学
7自由度机械臂，任务空间6维（位置+姿态），存在冗余自由度
雅可比矩阵（简化示例）
A = [[1, 0, 0, 0.5, 0.2, 0.1, 0],
[0, 1, 0, 0.3, 0.6, 0.2, 0],
[0, 0, 1, 0.1, 0.3, 0.8, 0],
[0, 0, 0, 1, 0, 0, 0.5],
[0, 0, 0, 0, 1, 0, 0.3],
[0, 0, 0, 0, 0, 1, 0.7]]
末端执行器期望速度
b = [0.1, 0.2, 0.05, 0, 0, 0] # [vx, vy, vz, ωx, ωy, ωz]
result = lsqminnorm(A, b)
print(f"关节角速度: {result}")
print(f"关节运动量: {np.linalg.norm(result):.4f}")
#关节角速度:
#[[ 0.07636763]
[ 0.17077049]
[ 0.00646669]
[ 0.03109522]
[ 0.01865713]
[ 0.04353331]
[-0.06219044]]
#关节运动量: 0.2052
#实际意义：在满足末端轨迹的同时最小化关节运动
示例4. 电路分析：节点电压法中的冗余方程
电路节点分析，存在线性相关方程
A = [[1, -1, 0, 0], # 节点1 KCL
[0, 1, -1, 0], # 节点2 KCL
[0, 0, 1, -1], # 节点3 KCL
[1, 0, 0, -1]] # 节点4 KCL（冗余）
电流源注入
b = [2, 0, -1, 1] # 各节点净电流
result = lsqminnorm(A, b)
print(f"节点电压: {result}")
print(f"电压范数: {np.linalg.norm(result):.4f}")
#节点电压:
#[[ 1.25]
[-0.75]
[-0.75]
[ 0.25]]
#电压范数: 1.6583
验证KCL满足
residuals = np.array(A) @ result - np.array(b)
print(f"KCL残差: {np.linalg.norm(residuals):.6f}")
#KCL残差: 6.324555
示例5. 图像处理：图像配准中的变换参数估计
从特征点对应关系估计仿射变换参数
每对点提供2个方程，6个未知参数
A = [[100, 50, 1, 0, 0, 0], # 点1 x坐标方程
[0, 0, 0, 100, 50, 1], # 点1 y坐标方程
[150, 80, 1, 0, 0, 0], # 点2 x坐标方程
[0, 0, 0, 150, 80, 1], # 点2 y坐标方程
[80, 120, 1, 0, 0, 0], # 点3 x坐标方程
[0, 0, 0, 80, 120, 1]] # 点3 y坐标方程
变换后的坐标
b = [105, 55, 160, 90, 85, 130]
result = lsqminnorm(A, b)
print(f"仿射变换参数: {result}")
#仿射变换参数:
#[[ 1.08536585]
[ 0.02439024]
[-4.75609756]
[ 0.04878049]
[ 1.08536585]
[-4.14634146]]
重组为变换矩阵
transform_matrix = result.reshape(2, 3)
print(f"变换矩阵:\n{transform_matrix}")
#[[1.08536585,0.02439024,-4.75609756]
[0.04878049,1.08536585,-4.14634146]]
#实际意义：估计图像配准中的最小变化变换
示例6. 金融工程：投资组合权重优化
在给定预期收益下寻找最小风险的投资组合权重
约束：权重和为1，预期收益达到目标
A = [[1, 1, 1, 1], # 权重和为1
[0.08, 0.12, 0.15, 0.10]] # 各资产预期收益率
目标：[总权重=1, 预期收益=0.11]
b = [1, 0.11]
result = lsqminnorm(A, b)
print(f"投资组合权重: {result}")
print(f"权重和: {np.sum(result):.6f}")
#投资组合权重:
#[[0.28037383]
[0.24299065]
[0.21495327]
[0.26168224]]
#权重和: 1.0
expected_return = np.array([0.08, 0.12, 0.15, 0.10]) @ result
print(f"预期收益率: {expected_return}")
#预期收益率: [0.11]
#实际意义：在满足收益要求下寻找最分散的投资组合
示例7. 地质学：重力反演问题
从地表重力测量反演地下密度分布（欠定问题）
每个测量点受多个地下单元影响
A = [[0.8, 0.3, 0.1, 0.05], # 测量点1的灵敏度
[0.3, 0.7, 0.4, 0.2], # 测量点2的灵敏度
[0.1, 0.4, 0.9, 0.6], # 测量点3的灵敏度
[0.05, 0.2, 0.6, 0.8]] # 测量点4的灵敏度
重力异常测量值
b = [1.2, 1.6, 2.0, 1.65]
result = lsqminnorm(A, b)
print(f"地下密度分布: {result}")
print(f"解的总质量: {np.sum(result):.4f}")
#地下密度分布:
#[[0.92460561]
[1.03728744]
[0.9897562 ]
[1.00307314]]
#解的总质量: 3.9547
#实际意义：寻找最平滑的地下密度分布解释观测数据
示例8. 化学计量学：光谱分析中的组分浓度估计
从混合光谱估计各组分浓度
各组分在特定波长的吸光度
A = [[0.8, 0.2, 0.1], # 波长1的吸光度
[0.3, 0.7, 0.4], # 波长2的吸光度
[0.1, 0.3, 0.9]] # 波长3的吸光度
混合样品在各波长的吸光度
b = [0.45, 0.55, 0.60]
result = lsqminnorm(A, b)
print(f"各组分浓度: {result}")
print(f"总浓度: {np.sum(result):.4f}")
#各组分浓度:
#[[0.42032967]
[0.31043956]
[0.51648352]]
#总浓度: 1.2473
验证重建光谱
reconstructed = np.array(A) @ result
print(f"光谱重建误差: {np.linalg.norm(reconstructed - b):.6f}")
#光谱重建误差: 0.264575
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
def linalg_min_norm(input_str):
"""
对标MATLAB的lsqminnorm函数，求解线性方程的最小范数最小二乘解
参数:
input_str: 输入字符串，格式为 (A, B) 的元组表达式
- A: 系数矩阵 (m×n)
- B: 右侧矩阵/向量 (m×k)
返回:
当输入有效时返回:
numpy数组: 最小范数解 (n×k)
否则返回错误信息字符串
示例:
>>> linalg_min_norm("([[1,2],[3,4],[5,6]], [7,8,9])")
array([-3.5 , 3.75])
"""
try:
expr = sp.sympify(input_str)
# 输入格式验证
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入错误: 需要 (A, B) 格式的元组，实际输入: {input_str}"
# 矩阵转换
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
B_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
# 有效性检查
if A_sym is None or B_sym is None:
return f"输入错误: 矩阵转换失败 - A: {expr[0]}, B: {expr[1]}"
# 转换为numpy数组
A = np.array(A_sym.tolist(), dtype=float)
B = np.array(B_sym.tolist(), dtype=float)
# 维度校验
if A.ndim != 2:
return f"维度错误: A 必须是二维矩阵，实际维度: {A.ndim}"
if B.ndim not in (1, 2):
return f"维度错误: B 必须是一维向量或二维矩阵，实际维度: {B.ndim}"
if A.shape[0] != B.shape[0]:
return f"维度不匹配: A的行数({A.shape[0]}) != B的行数({B.shape[0]})"
# 计算伪逆（Moore-Penrose pseudoinverse）
A_pinv = np.linalg.pinv(A)
# 矩阵乘法求解
if B.ndim == 1:
X = A_pinv @ B # 向量情况 (n,)
else:
X = A_pinv @ B # 矩阵情况 (n, k)
return X
except sp.SympifyError as e:
return f"表达式解析错误: {str(e)}"
except np.linalg.LinAlgError as e:
return f"线性代数计算错误: {str(e)}"
except Exception as e:
return f"未知错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 欠定系统 (无穷多解，求最小范数解)
print("示例1: 欠定系统")
case1 = linalg_min_norm("([[1, 1], [1, 1]], [2, 2])")
print("解:", case1)
# [[1.]
# [1.]]
# 示例2: 超定系统 (最小二乘解)
print("\n示例2: 超定系统")
case2 = linalg_min_norm("([[1,2], [3,4], [5,6]], [7,8,9])")
print("解:", case2)
# [[-6. ]
# [ 6.5]]
# 示例3: 多右侧项
print("\n示例3: 多右侧项")
case3 = linalg_min_norm("([[1,2], [3,4]], [[5,6], [7,8]])")
print("解:\n", case3)
# [[-3. -4.]
# [ 4. 5.]]
非负最小二乘问题
x = lsqonneg(C，d)返回在x≥0的情况下最小化范数(C*x-d)的向量x.参数C和d必须是实数的.
C — 实数矩阵.
d — 实数,向量.
示例1. 化学分析：光谱组分浓度估计
从混合光谱估计各化学组分的浓度（浓度不能为负）
C = [[0.8, 0.2, 0.1], # 组分1在各波长的吸光系数
[0.3, 0.7, 0.4], # 组分2在各波长的吸光系数
[0.1, 0.3, 0.9]] # 组分3在各波长的吸光系数
d = [0.45, 0.55, 0.60] # 混合样品在各波长的吸光度
result = lsqnonneg(C, d)
print(f"各组分浓度: {result}")
print(f"总浓度: {sum(result):.4f}")
#各组分浓度: [0.42032967032967045, 0.31043956043956056, 0.5164835164835168]
#总浓度: 1.2473
验证重建光谱
reconstructed = np.array(C) @ np.array(result)
error = np.linalg.norm(reconstructed - d)
print(f"光谱重建误差: {error:.6f}")
#光谱重建误差: 0.0
# 实际意义：化学浓度必须是正值，负浓度没有物理意义
示例2. 经济学：投入产出分析
列昂惕夫投入产出模型，计算各部门产出
技术系数矩阵
C = [[0.2, 0.3, 0.1], # 农业部门消耗
[0.1, 0.4, 0.2], # 工业部门消耗
[0.05, 0.1, 0.3]] # 服务业消耗
d = [100, 150, 80] # 最终需求
result = lsqnonneg(C, d)
print(f"各部门总产出: {result}")
print(f"GDP估算: {sum(d):.1f} (最终需求总和)")
#各部门总产出: [0.0, 283.33333333333326, 173.80952380952382]
#GDP估算: 330.0 (最终需求总和)
#实际意义：经济产出必须是正值，负产出没有经济意义
示例3. 环境科学：污染源解析
从监测点数据解析各污染源的贡献率
C = [[0.8, 0.1, 0.05], # 监测点1对各污染源的敏感度
[0.2, 0.7, 0.1], # 监测点2对各污染源的敏感度
[0.05, 0.3, 0.9]] # 监测点3对各污染源的敏感度
d = [1.2, 0.9, 1.5] # 各监测点的污染物浓度
result = lsqnonneg(C, d)
print(f"各污染源贡献率: {result}")
print(f"主要污染源: 源{np.argmax(result) + 1} (贡献率{max(result):.2%})")
#各污染源贡献率: [1.326145552560647, 0.7132075471698108, 1.3552560646900274]
#主要污染源: 源3 (贡献率135.53%)
#实际意义：污染贡献率不能为负值
示例4. 医学成像：CT重建
从投影数据重建组织密度分布
C = [[1, 0, 0, 1], # 投影1经过的体素
[0, 1, 1, 0], # 投影2经过的体素
[1, 1, 0, 0], # 投影3经过的体素
[0, 0, 1, 1]] # 投影4经过的体素
d = [2.5, 2.0, 2.2, 1.8] # 投影测量值
result = lsqnonneg(C, d)
print(f"组织密度分布: {result}")
print(f"平均密度: {np.mean(result):.4f}")
#组织密度分布: [2.324999999999999, 0.0, 1.875, 0.04999999999999971]
#平均密度: 1.0625
#实际意义：组织密度不能为负值
示例5. 金融学：资产配置权重
在给定约束下寻找最优资产配置权重（权重不能为负）
C = [[0.08, 0.12, 0.15, 0.10], # 各资产预期收益率
[1, 1, 1, 1]] # 权重和为1的约束
d = [0.11, 1] # 目标收益率和权重约束
result = lsqnonneg(C, d)
print(f"资产配置权重: {result}")
print(f"权重和: {sum(result):.6f}")
expected_return = np.array([0.08, 0.12, 0.15, 0.10]) @ np.array(result)
print(f"预期收益率: {expected_return:.4f}")
#资产配置权重: [0.5714285714285707, 0.0, 0.4285714285714292, 0.0]
#权重和: 1.000000
#预期收益率: 0.1100
#实际意义：投资权重不能为负（不允许卖空）
示例6. 供应链管理：资源分配
将有限资源分配给不同产品以最大化满足需求
C = [[2, 1, 3], # 产品1的资源消耗系数
[1, 3, 2], # 产品2的资源消耗系数
[3, 2, 1]] # 产品3的资源消耗系数
d = [100, 150, 120] # 可用资源总量
result = lsqnonneg(C, d)
print(f"各产品产量: {result}")
#各产品产量: [10.555555555555559, 37.22222222222222, 13.888888888888902]
resource_used = np.array(C) @ np.array(result)
utilization = resource_used / np.array(d)
print(f"资源利用率: {utilization}")
#资源利用率: [1. 1. 1.]
#实际意义：产量不能为负值
示例7. 营养学：膳食配方优化
设计满足营养需求的膳食配方
C = [[50, 30, 20], # 食物A的营养含量（蛋白质,碳水,脂肪）
[20, 60, 10], # 食物B的营养含量
[10, 20, 70]] # 食物C的营养含量
d = [40, 50, 30] # 每日营养需求
result = lsqnonneg(C, d)
print(f"各食物摄入量: {result}")
#各食物摄入量: [0.3057324840764332, 0.7006369426751593, 0.18471337579617833]
nutrition_intake = np.array(C) @ np.array(result)
print(f"营养摄入量: {nutrition_intake}")
print(f"与需求的差异: {nutrition_intake - d}")
#营养摄入量: [40. 50. 30.]
#与需求的差异: [7.10542736e-15 7.10542736e-15 0.00000000e+00]
#实际意义：食物摄入量不能为负
示例8. 机器学习：非负矩阵分解
图像数据的非负矩阵分解（简化示例）
C = [[0.9, 0.1, 0.2], # 基图像1在各像素的值
[0.1, 0.8, 0.1], # 基图像2在各像素的值
[0.2, 0.1, 0.7]] # 基图像3在各像素的值
d = [0.8, 0.6, 0.7] # 目标图像像素值
result = lsqnonneg(C, d)
print(f"基图像权重: {result}")
#基图像权重: [0.6630434782608692, 0.576086956521739, 0.7282608695652173]
reconstructed = np.array(C) @ np.array(result)
reconstruction_error = np.linalg.norm(reconstructed - d)
print(f"重建误差: {reconstruction_error:.6f}")
#重建误差: 0.0
#实际意义：在NMF中，基和系数都必须是非负的
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.optimize import nnls
def least_squares_nonnegative_linear(input_str):
"""
对标MATLAB的lsqnonneg函数，求解非负线性最小二乘问题
参数:
input_str: 输入字符串，格式为 (C, D) 的元组表达式
- C: 系数矩阵 (m×n)
- D: 右侧向量 (m×1)
返回:
成功时返回字典:
{
"x": 非负解向量 (n×1),
"residual": 残差范数,
"message": 状态描述
}
失败时返回错误信息字符串
示例:
>>> least_squares_nonnegative_linear("([[1,2],[3,4]], [5,6])")
{'x': [0.0, 1.3], 'residual': 0.8944271909999159, 'message': '成功收敛'}
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str)
# 输入格式验证
if not isinstance(expr, tuple) or len(expr) != 2:
return "输入错误: 需要 (C, D) 格式的元组"
# 矩阵转换
C_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
D_sym = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
# 有效性检查
if C_sym is None or D_sym is None:
return "输入错误: 矩阵转换失败"
if C_sym.shape[0] != D_sym.shape[0]:
return f"维度不匹配: C的行数({C_sym.shape[0]}) != D的长度({D_sym.shape[0]})"
# 转换为numpy数组
C = np.array(C_sym.tolist(), dtype=float)
D = np.array(D_sym.tolist(), dtype=float).flatten() # 确保转换为1D数组
# 维度校验
if C.ndim != 2:
return "系数矩阵C必须是二维矩阵"
if len(D.shape) != 1:
return "右侧向量D必须是一维数组"
# 执行非负最小二乘计算
x, residual = nnls(C, D)
return x.tolist()
'''
return {
"x": x.tolist(),
"residual": float(residual),
"message": "成功收敛" if residual < 1e-10 else "达到迭代限制"
}
'''
except sp.SympifyError as e:
return f"表达式解析错误: {str(e)}"
except ValueError as e:
return f"数值错误: {str(e)}"
except Exception as e:
return f"未知错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 标准非负解
print("示例1: 简单情况")
case1 = least_squares_nonnegative_linear("[[1,0],[1,0],[0,1]],[2,1,1]")
print("残差:", case1)
# [1.4999999999999996, 1.0]
# 示例2: 需要强制归零的解
print("\n示例2: 边界解")
case2 = least_squares_nonnegative_linear("([[2,1],[1,2]], [3,3])")
print("解:", case2)
# [0.9999999999999998, 1.0]
LU矩阵分解
[L,U,P] = lu(A, outputForm) 将满矩阵或稀疏矩阵A分解为一个上三角矩阵U和一个经过置换的下三角矩阵L,使得A=L*U, 及置换矩阵P,并满足A=P'*L*U.
A — 输入矩阵,矩阵.
outputForm — 置换输出的形状, 'matrix' (默认) | 'vector'
L — 下三角因子,矩阵.
U — 上三角因子,矩阵
P — 行置换,向量,矩阵
示例1. 线性方程组求解：电路分析
电路节点分析方程
A = [[4, -1, -1], # 节点1的KCL方程
[-1, 4, -1], # 节点2的KCL方程
[-1, -1, 4]] # 节点3的KCL方程
b = [10, 0, 0] # 电流源注入
使用LU分解求解 Ax = b
result = lu(A)
L = np.array(result["L"])
U = np.array(result["U"])
P = np.array(result["P"])
# 解 Ly = Pb
Pb = P @ np.array(b)
y = np.linalg.solve(L, Pb)
# 解 Ux = y
x = np.linalg.solve(U, y)
print(f"节点电压: {x}")
print(f"验证残差: {np.linalg.norm(np.array(A) @ x - b):.6f}")
#节点电压: [3. 1. 1.]
#验证残差: 0.0
#实际意义：LU分解可以高效求解大型线性方程组
示例2. 数值计算：矩阵求逆
使用LU分解计算矩阵逆
A = [[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]]
result = lu(A)
L = np.array(result["L"])
U = np.array(result["U"])
P = np.array(result["P"])
通过解多个线性方程组求逆
n = len(A)
A_inv = np.zeros((n, n))
for i in range(n):
e = np.zeros(n)
e[i] = 1
# 解 Ly = Pe
y = np.linalg.solve(L, P @ e)
# 解 Ux = y
A_inv[:, i] = np.linalg.solve(U, y)
print(f"矩阵逆:\n{A_inv}")
print(f"验证 A*A⁻¹:\n{np.array(A) @ A_inv}")
#矩阵逆:
#[[-1,4,3]
[ 1,-2,-2]
[-1,5,4]]
#验证 A*A⁻¹:
#[[0.00000000e+00,1.00000000e+00,8.88178420e-16]
[-4.44089210e-16,1.77635684e-15,1.00000000e+00]
[1.00000000e+00,0.00000000e+00,0.00000000e+00]]
示例3. 结构工程：刚度矩阵分解
桁架结构的刚度矩阵
A = [[3, -1, -1, 0],
[-1, 2, 0, -1],
[-1, 0, 2, -1],
[0, -1, -1, 2]]
result = lu(A)
L = np.array(result["L"])
U = np.array(result["U"])
print(f"刚度矩阵条件数: {np.linalg.cond(A):.4f}")
print(f"L矩阵结构:\n{L}")
print(f"U矩阵结构:\n{U}")
#刚度矩阵条件数: 23.2998
#L矩阵结构:
#[[ 1. 0. 0. 0. ]
[-0.33333333 1. 0. 0. ]
[-0.33333333 -0.2 1. 0. ]
[ 0. -0.6 -0.75 1. ]]
#U矩阵结构:
#[[ 3. -1. -1. 0. ]
[ 0. 1.66666667 -0.33333333 -1. ]
[ 0. 0. 1.6 -1.2 ]
[ 0. 0. 0. 0.5 ]]
#实际意义：在有限元分析中，LU分解用于求解结构位移
示例4. 计算流体力学：压力泊松方程
离散化的压力泊松方程（5点差分格式）
A = [[4, -1, 0, -1, 0],
[-1, 4, -1, 0, -1],
[0, -1, 4, -1, 0],
[-1, 0, -1, 4, -1],
[0, -1, 0, -1, 4]]
result = lu(A, vector)
print(f"置换向量: {result['perm']}")
print(f"L矩阵:\n{result['L']}")
print(f"U矩阵:\n{result['U']}")
#置换向量: [0, 1, 2, 3, 4]
#L矩阵:
#[[1.0, 0.0, 0.0, 0.0, 0.0],
[-0.25, 1.0, 0.0, 0.0, 0.0],
[0.0, -0.26666666666666666, 1.0, 0.0, 0.0],
[-0.25, -0.06666666666666667, -0.2857142857142857, 1.0, 0.0],
[0.0, -0.26666666666666666, -0.07142857142857142, -0.33333333333333326, 1.0]]
#U矩阵:
#[[4.0, -1.0, 0.0, -1.0, 0.0],
[0.0, 3.75, -1.0, -0.25, -1.0],
[0.0, 0.0, 3.7333333333333334, -1.0666666666666667, -0.26666666666666666],
[0.0, 0.0, 0.0, 3.428571428571429, -1.1428571428571428],
[0.0, 0.0, 0.0, 0.0, 3.3333333333333335]]
#实际意义：在CFD中，LU分解用于求解压力修正方程
示例5. 经济学：投入产出模型
I - A 矩阵，其中A是技术系数矩阵
A = [[0.8, -0.3, -0.2],
[-0.2, 0.9, -0.1],
[-0.1, -0.2, 0.8]]
result = lu(A)
L = np.array(result["L"])
U = np.array(result["U"])
P = np.array(result["P"])
计算行列式（用于判断系统可解性）
det_A = np.prod(np.diag(U)) # det(A) = det(U) 因为det(L)=1
print(f"矩阵行列式: {det_A:.6f}")
#矩阵行列式: 0.483000
if abs(det_A) > 1e-10:
print("系统有唯一解")
#系统有唯一解
#实际意义：判断经济系统是否可解
else:
print(f"错误: {result}")
示例6. 图像处理：图像变形
从对应点求解仿射变换参数
每对点提供2个方程，需要6个参数
points_src = [[10, 10], [50, 10], [10, 50]] # 原图点
points_dst = [[15, 12], [55, 8], [12, 55]] # 目标点
构建系数矩阵 A * x = b
A = []
b = []
for (x1, y1), (x2, y2) in zip(points_src, points_dst):
A.append([x1, y1, 1, 0, 0, 0])
A.append([0, 0, 0, x1, y1, 1])
b.extend([x2, y2])
result = lu(A)
L = np.array(result["L"])
U = np.array(result["U"])
P = np.array(result["P"])
求解变换参数
Pb = P @ np.array(b)
y = np.linalg.solve(L, Pb)
x = np.linalg.solve(U, y)
print(f"仿射变换参数: {x}")
transform_matrix = x.reshape(2, 3)
print(f"变换矩阵:\n{transform_matrix}")
#仿射变换参数:
#[1.00000000e+00,1.00000000e+00,-5.00000000e+00,-1.00000000e-01,-3.70074342e-17,1.30000000e+01]
#变换矩阵:
#[[1.00000000e+00,1.00000000e+00,-5.00000000e+00]
[-1.00000000e-01,-3.70074342e-17,1.30000000e+01]]
示例7. 机器学习：协方差矩阵分析
协方差矩阵（核矩阵）的LU分解
X = [[1, 2], [2, 3], [3, 4], [4, 5]]
n = len(X)
构建RBF核矩阵
def rbf_kernel(x1, x2, length_scale=1.0):
return np.exp(-0.5 * np.sum((np.array(x1) - np.array(x2)) ** 2) / length_scale ** 2)
K = np.zeros((n, n))
for i in range(n):
for j in range(n):
K[i, j] = rbf_kernel(X[i], X[j])
result = lu(K)
L = np.array(result["L"])
U = np.array(result["U"])
验证分解正确性
reconstruction = L @ U
error = np.linalg.norm(reconstruction - K)
print(f"重建误差: {error:.6f}")
#重建误差: 0.0
#实际意义：在高斯过程中，LU分解用于计算预测分布
示例8. 符号计算：理论分析
符号矩阵的LU分解，用于理论分析
A_sym = [[a, b, c],
[d, e, f],
[g, h, i]]
result = lu(A_sym)
print("符号L矩阵:")
for row in result["L"]:
print(row)
#符号L矩阵:
#[1, 0, 0]
#[d/a, 1, 0]
#[g/a, (h - b*g/a)/(e - b*d/a), 1]
print("\n符号U矩阵:")
for row in result["U"]:
print(row)
#符号U矩阵:
#[a, b, c]
#[0, e - b*d/a, f - c*d/a]
#[0, 0, i - (f - c*d/a)*(h - b*g/a)/(e - b*d/a) - c*g/a]
#实际意义：符号分解用于理论推导和公式验证
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from scipy.linalg import lu
def lu_decomposition_matrix(input_str, output_form='matrix'):
"""
对标MATLAB的lu函数，执行LU分解
参数:
input_str: 输入矩阵的字符串表达式
output_form: 输出形式，可选：
'matrix' - 返回 (L, U, P) 三个矩阵（默认）
'compact' - 返回组合的紧凑格式
'vector' - 返回行交换索引向量
返回:
成功时返回字典包含：
- L: 下三角矩阵
- U: 上三角矩阵
- P: 置换矩阵/向量
失败时返回错误信息字符串
示例:
>>> lu_decomposition_matrix("[[4,3],[6,3]]")
{'L': [[1.0, 0.0], [0.6666666666666666, 1.0]],
'U': [[6.0, 3.0], [0.0, 1.0]],
'P': [[0,1],[1,0]]}
"""
try:
# 解析输入表达式
expr = sp.sympify(input_str)
A = sp.Matrix(expr) if isinstance(expr, list) else None
if A is None:
return "错误：输入不是有效矩阵"
# 符号矩阵处理
if A.is_symbolic() or any(element.has(sp.I) for element in A):
# SymPy符号分解
if not A.is_square:
return "错误：符号矩阵必须是方阵"
L, U, swaps = A.LUdecomposition()
P = sp.eye(A.rows).permute_rows(swaps)
return {
"L": L.tolist(),
"U": U.tolist(),
"P": P.tolist(),
"swaps": swaps
}
# 数值矩阵处理
else:
A_np = np.array(A.tolist(), dtype=float)
# 执行SciPy分解
P_np, L_np, U_np = lu(A_np)
# 处理输出形式
if output_form == 'vector':
# 将置换矩阵转换为索引向量
perm = np.argmax(P_np, axis=1)
return {
"L": L_np.tolist(),
"U": U_np.tolist(),
"perm": perm.tolist()
}
elif output_form == 'compact':
# 组合紧凑格式
LU = P_np @ A_np
return {
"LU": LU.tolist(),
"pivots": np.arange(A_np.shape[0]).tolist()
}
else:
return {
"L": L_np.tolist(),
"U": U_np.tolist(),
"P": P_np.tolist()
}
except sp.SympifyError as e:
return f"表达式解析错误: {str(e)}"
except ValueError as e:
return f"数值错误: {str(e)}"
except np.linalg.LinAlgError as e:
return f"线性代数错误: {str(e)}"
except Exception as e:
return f"未知错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1：基础数值分解
print("示例1：2x2数值矩阵")
case1 = lu_decomposition_matrix("[[4,3],[6,3]]")
print("L:\n", case1["L"])
# [[1.0, 0.0],
# [0.6666666666666666, 1.0]]
print("U:\n", case1["U"])
# [[6.0, 3.0],
# [0.0, 1.0]]
print("P:\n", case1["P"])
# [[0.0, 1.0],
# [1.0, 0.0]]
# 示例2：符号矩阵分解
print("\n示例2：3x3符号矩阵")
case2 = lu_decomposition_matrix("[[a,2,3],[4,b,6],[7,8,c]]")
print("符号L:\n", case2["L"])
# [[1, 0, 0],
# [a/4, 1, 0],
# [7/4, (8 - 7*b/4)/(-a*b/4 + 2), 1]]
print("符号U:\n", case2["U"])
# [[4, b, 6],
# [0, -a*b/4 + 2, 3 - 3*a/2],
# [0, 0, c - (3 - 3*a/2)*(8 - 7*b/4)/(-a*b/4 + 2) - 21/2]]
# 示例3：输出置换向量
print("\n示例3：置换向量输出")
case3 = lu_decomposition_matrix("[[1,2],[3,4]]", output_form='vector')
print("置换索引:", case3["perm"])
# [1, 0]
广义卢卡斯函数
Lucas(n, [z]) 返回卢卡斯数或卢卡斯多项式
A — 整数,阶数.
z - 标量, 符号标量, 返回卢卡斯多项式
示例1. 数论：卢卡斯数列性质验证
验证卢卡斯数列的前10项
for n in range(10):
result = Lucas(n)
print(f"L({n}) = {result}")
#L(0) = 2
#L(1) = 1
#L(2) = 3
#L(3) = 4
#L(4) = 7
#L(5) = 11
#L(6) = 18
#L(7) = 29
#L(8) = 47
#L(9) = 76
验证卢卡斯数与斐波那契数的关系：L(n) = F(n-1) + F(n+1)
for n in range(2, 8):
lucas_n = Lucas(n)
# 注意：这里需要斐波那契函数，假设已定义
# fib_n_minus_1 = fibonacci(n-1)
# fib_n_plus_1 = fibonacci(n+1)
# print(f"L({n}) = F({n-1}) + F({n+1}) = {fib_n_minus_1} + {fib_n_plus_1} = {fib_n_minus_1 + fib_n_plus_1}")
print(f"L({n}) = {lucas_n}")
#L(2) = 3
#L(3) = 4
#L(4) = 7
#L(5) = 11
#L(6) = 18
#L(7) = 29
示例2. 组合数学：计数问题
卢卡斯数可用于计算某些圆排列问题的解
例如：n个点的圆排列中，某些特定模式的计数
points = [3, 4, 5, 6, 7]
for n in points:
result = Lucas(n)
print(f"{n}个点的特定圆排列数: {result}")
#3个点的特定圆排列数: 4
#4个点的特定圆排列数: 7
#5个点的特定圆排列数: 11
#6个点的特定圆排列数: 18
#7个点的特定圆排列数: 29
示例3. 计算机图形学：黄金比例相关计算
卢卡斯数与黄金比例 φ 密切相关
import math
phi = (1 + math.sqrt(5)) / 2
验证卢卡斯数的闭式公式
for n in range(1, 6):
lucas_closed = phi ** n + (-phi) ** (-n) # 理论公式
lucas_calc = Lucas(n)
print(f"n={n}: 计算值={lucas_calc}, 理论值≈{lucas_closed:.2f}, 误差={abs(lucas_calc - lucas_closed):.6f}")
#n=1: 计算值=1, 理论值≈1.00, 误差=0.000000
#n=2: 计算值=3, 理论值≈3.00, 误差=0.000000
#n=3: 计算值=4, 理论值≈4.00, 误差=0.000000
#n=4: 计算值=7, 理论值≈7.00, 误差=0.000000
#n=5: 计算值=11, 理论值≈11.00, 误差=0.000000
示例4. 物理学：量子系统能级模拟
在某些量子系统中，能级分布与卢卡斯数相关
energy_levels = [1, 2, 3, 4, 5]
for level in energy_levels:
# 简化的能级计算模型
energy = Lucas(level)
print(f"能级 {level}: {energy} 单位")
#能级 1: 1 单位
#能级 2: 3 单位
#能级 3: 4 单位
#能级 4: 7 单位
#能级 5: 11 单位
示例5. 卢卡斯多项式应用
计算不同阶数的卢卡斯多项式在特定点的值
x_values = [0.5, 1, 1.5, 2]
n_values = [2, 3, 4]
for n in n_values:
print(f"\n{n}阶卢卡斯多项式:")
for x in x_values:
result = Lucas(n, x)
print(f" L_{n}({x}) = {result}")
#2阶卢卡斯多项式:
# L_2(0.5) = 2.25
# L_2(1) = 3
# L_2(1.5) = 4.25
# L_2(2) = 6
#3阶卢卡斯多项式:
# L_3(0.5) = 1.625
# L_3(1) = 4
# L_3(1.5) = 7.875
# L_3(2) = 14
#4阶卢卡斯多项式:
# L_4(0.5) = 3.0625
# L_4(1) = 7
# L_4(1.5) = 16.0625
# L_4(2) = 34
示例6. 矩阵输入的卢卡斯计算
对矩阵中的每个元素计算对应的卢卡斯数
matrix = [[1, 2, 3],
[4, 5, 6],
[7, 8, 9]]
result = Lucas(matrix)
print("\n卢卡斯变换后的矩阵:")
for row in lucas_matrix:
print(row)
#卢卡斯变换后的矩阵:
#[1, 3, 4]
#[7, 11, 18]
#[29, 47, 76]
示例7. 符号计算：卢卡斯多项式推导
for n in range(1, 5):
poly_expr = Lucas(n, x)
print(f"L_{n}(x) = {poly_expr}")
#L_1(x) = x
#L_2(x) = x**2 + 2
#L_3(x) = x*(x**2 + 3)
#L_4(x) = x**4 + 4*x**2 + 2
验证多项式性质：L_n(1) = L_n (卢卡斯数)
print("\n验证 L_n(1) = L_n:")
for n in range(1, 6):
poly_at_1 = Lucas(n, 1)
lucas_num = Lucas(str(n))
print(f"L_{n}(1) = {poly_at_1}, L_{n} = {lucas_num}, 相等: {poly_at_1 == lucas_num}")
#L_1(1) = 1, L_1 = 1, 相等: True
#L_2(1) = 3, L_2 = 3, 相等: True
#L_3(1) = 4, L_3 = 4, 相等: True
#L_4(1) = 7, L_4 = 7, 相等: True
#L_5(1) = 11, L_5 = 11, 相等: True
示例8. 数值分析：收敛性测试
研究卢卡斯数列的增长特性
terms = 15
lucas_terms = []
ratios = []
for n in range(terms):
term = Lucas(n)
lucas_terms.append(term)
if n >= 2:
ratio = lucas_terms[n] / lucas_terms[n - 1]
ratios.append(ratio)
print(f"前{terms}个卢卡斯数: {lucas_terms}")
print(f"连续项比值: {ratios}")
print(f"比值收敛到黄金比例 φ ≈ {phi:.6f}")
print(f"最后5个比值的平均值: {np.mean(ratios[-5:]):.6f}")
#前15个卢卡斯数: [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843]
#连续项比值: [3, 4/3, 7/4, 11/7, 18/11, 29/18, 47/29, 76/47, 123/76, 199/123, 322/199, 521/322, 843/521]
#比值收敛到黄金比例 φ ≈ 1.618034
#最后5个比值的平均值: 1.618090
示例9. 应用数学：基于卢卡斯数列的优化
在某些离散优化问题中，卢卡斯数可作为目标函数值
def objective_function(params):
"""基于卢卡斯数的简单目标函数"""
x, y = params
# 使用卢卡斯数构造非线性目标
term1 = Lucas(int(abs(x))) if abs(x) < 10 else 0
term2 = Lucas(int(abs(y))) if abs(y) < 10 else 0
return term1 + term2
测试几个点
test_points = [(1, 2), (3, 1), (2, 3)]
for point in test_points:
value = objective_function(point)
print(f"f{point} = {value}")
#f(1, 2) = 4
#f(3, 1) = 5
#f(2, 3) = 7
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def lucas_number(input_str):
"""
对标Wolfram Lucas数的实现，计算第n个卢卡斯数
https://mathworld.wolfram.com/LucasNumber.html
参数:
input_str: 输入的字符串表达式，可以是：
- 整数（如'5'）
- 整数矩阵（如'[[0, 1], [2, 3]]'）
- 符号表达式（需包含整数约束）
返回:
计算结果，类型与输入对应：
- 整数输入：返回整数结果
- 矩阵输入：返回同型矩阵，逐元素计算结果
- 符号表达式：返回符号表达式
- 错误输入：返回错误描述字符串
示例:
>>> lucas_number('0')
2
>>> lucas_number('[[1, 2], [3, 4]]')
Matrix([[1, 3], [4, 7]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def eval_lucas(n):
"""卢卡斯数核心计算公式，支持符号运算"""
if not n.is_integer:
return None
# 使用闭式公式计算，与Wolfram实现一致
sqrt5 = sp.sqrt(5)
return sp.simplify(((1 + sqrt5) / 2) ** n + ((1 - sqrt5) / 2) ** n)
def lucas_polynomial(n_val, x_val):
"""
计算第n阶Lucas多项式 L_n(x)
参数:
n (int): 多项式阶数
x (symbol/float): 变量或具体数值
返回:
Expr: 符号表达式或数值结果
示例:
>>> lucas_polynomial(3, sp.Symbol('x'))
x**3 + 3*x
>>> lucas_polynomial(2, 1)
3 # 1^2 + 2 = 3 (对应Lucas数 L_2=3)
"""
if n_val.is_integer:
sqrt_term = sp.sqrt(x_val ** 2 + 4)
# 构建表达式
result_expr = 2 ** (-n_val) * ((x_val - sqrt_term) ** n_val + (x_val + sqrt_term) ** n_val)
return sp.simplify(result_expr)
else:
return None
if isinstance(expr, tuple):
n = expr[0] # base b
x = expr[1]
result = lucas_polynomial(n, x)
# 处理标量输入
elif expr.is_integer:
result = eval_lucas(expr)
# 处理符号表达式（需包含整数约束）
elif expr.free_symbols:
# 添加整数约束条件
n = list(expr.free_symbols)[0]
if sp.ask(sp.Q.integer(n)):
error = True
else:
result = eval_lucas(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
if __name__ == "__main__":
# 示例1: 标量数值测试
print("示例1 标量输入:")
print(lucas_number('6'))
# 18
print(lucas_number('1'))
# 1
print(lucas_number('5'))
# 11