Laplace数值逆运算的讨论
■ 问题提出
在博文 逆Laplace数值逆变换 给出了数值计算Laplace逆变换的简易程序。其中存在以下几个问题需要讨论:
- 问题1: 程序实现过程原理以及优化;
- 问题2: 运算参数: σ \sigma σ, ω max \omega_{\max} ωmax, N i n t N_{{\mathop{\rm int}} } Nint对于积分数值的影响。从上篇博文中明显看到有些计算出的结果大大偏离的实际值。比如其中的函数 sin ( t ) \sin \left( t \right) sin(t)所计算出来的幅值超过了1.
01程序实现原理以及优化
1.Laplace逆运算
所以,算法的核心是对
F
(
s
)
F\left( s \right)
F(s)进行傅里叶反变换,然后在乘以
e
σ
t
e^{\sigma t}
eσt。
由于确认变换后的函数 f ( t ) f\left( t \right) f(t)是实函数,因此,为了节省计算时间,只对傅里叶反变换的积分,进行正半轴的积分,同时积分的上限有参数 ω max \omega _{\max } ωmax决定。对积分的数值取齐实部,再乘以2便可以得到 f ( t ) f\left( t \right) f(t)。
2.Python积分程序实现优化
使用梯形积分来实现函数的积分,可以获得更精确的积分值。理论分析可知:
其中
Δ
x
=
(
b
−
a
)
/
N
\Delta x = \left( {b - a} \right)/N
Δx=(b−a)/N 以及
x
i
=
a
+
i
Δ
x
x_i = a + i\Delta x
xi=a+iΔx。那么积分误差上限为:
其中,
∣
f
′
′
(
x
)
∣
≤
K
2
\left| {f^{''} \left( x \right)} \right| \le K_2
∣∣∣f′′(x)∣∣∣≤K2,
x
∈
[
a
,
b
]
x \in \left[ {a,b} \right]
x∈[a,b]。
▲ 图2.1 梯形积分方法示意图
计算 T N ( f ) T_N \left( f \right) TN(f)可以由两种方式:
方式1: 对左、右黎曼积分加权平均:
方式2: 利用如下的公式计算:
最后实现的代码为:
def trapz(f, a, b, N=50):
x = linspace(a, b, N+1)
y = f(x)
y_right = y[1:]
y_left = y[:-1]
dx = (b-a) / N
T = dx/2 * sum(y_right + y_left)
return T
利用上述公式,对 sin ( t ) \sin \left( t \right) sin(t)进行积分测试:
printf(trapz(sin, 0, pi/2, 1000))
所得到结果为:0.9999997943832332。
可以看到使用n=1000对应的结果后面的积分精度达到了小数点后面6位9的小数点的位数。
3.实现Laplace逆运算
通过上面梯形积分方法,实现Laplace数值逆变换,具体的子程序如下面所示。
#------------------------------------------------------------
def invlt(t, fs, sigma, omiga, nint):
omigadim = linspace(0, omiga, nint+1, endpoint=True)
y = [(exp(1j*o*t) * fs(sigma+1j*o)).real for o in omigadim]
y_left = y[:-1]
y_right = y[0:]
T = sum(y_right + y_left) * omiga/nint
return exp(sigma*t) * T/ pi / 2
#------------------------------------------------------------
def fs(s):
return 1/(s*s+1)
#------------------------------------------------------------
sigma = 0.2
omiga=200
nint=omiga*50
tdim = linspace(0, 2*pi* 3, 200)
ft = [invlt(t, fs, sigma, omiga, nint) for t in tdim]
02一些基本函数的实验
下面通过对一些基本常见函数的laplace变换,来测试一下上述程序的性能。
Ⅰ.sin(t)
sigma=0.2, omiga=200, nint=omiga*50
Ⅱ.exp(-t)
sigam=-1+0.1, omiga=200, nint=omiga*50
Ⅲ.u(t)
Ⅳ.u(t-1)
Ⅴ.周期化脉冲信号
Ⅵ.三角形
▲ omiga从10变化到1000
※ 结论
通过原理分析,可以获得建议的Laplace数值逆运算的正确的PYTHON算法程序。
这个程序是直接对Laplace反变换公式利用梯形积分方法获得计算结果。通过对几种常见的信号Laplace的反变换,验证了这个算法的正确性。
通过在此过程,可以看到,对于参数sigma, omiga, nint对于计算结果还是有很大的影响。另外,对于时间t,只能在比较小的范围内有效,当t超过一定长度,前面所计算的结果都会出现比较大的误差。
▲ 对于方波进行Laplace数值反变换的结果
§ 参考文献
相关文章
- shell脚本--数值比较
- 【FAST API】路径参数和数值校验¶
- Java实现 LeetCode 363 矩形区域不超过 K 的最大数值和
- postgresql常用数据类型:数值、日期、字符串类型
- 【第3版emWin教程】第13章 emWin6.x数值显示
- 华为OD机试 - Excel单元格数值统计(Java & JS & Python)
- ML之回归预测之Lasso:利用Lasso算法对对红酒品质wine数据集解决回归(实数值评分预测)问题—优化模型【增加新(组合)属性】
- EL之Boosting之GB(DTR):利用梯度提升法解决回归(对多变量的数据集+实数值评分预测)问题
- EL之RF(随机性的Bagging+DTR):利用随机选择属性的bagging方法解决回归(对多变量的数据集+实数值评分预测)问题
- 常微分方程初值问题数值解法[完整公式](Python)
- 【C】 数值统计
- 高等数值计算方法学习笔记第4章第三部分【数值积分(数值微分)】
- VL11 4位数值比较器电路