29、拉普拉斯变换与逆变换在常微分方程求解中的应用

拉普拉斯变换与逆变换在常微分方程求解中的应用

一、拉普拉斯变换与逆变换概述

拉普拉斯变换及其逆变换是解决常微分方程（ODE）的有效工具。在求解ODE时，一般步骤如下：
1. 对给定的ODE进行拉普拉斯变换，将时域问题转换为复频域（s域）问题。
2. 求解s域中的代数方程，得到关于Y(s)的表达式。
3. 对得到的Y(s)表达式进行逆拉普拉斯变换，以获得原问题在时域的最终解。在这一步，通常需要将Y(s)表达式分解为部分分式，以便使用逆拉普拉斯变换对应表。

下面是求解ODE的拉普拉斯变换和逆变换的流程图：

graph TD;
    A[给定ODE] --> B[拉普拉斯变换];
    B --> C[s域代数方程];
    C --> D[求解Y(s)];
    D --> E[部分分式分解];
    E --> F[逆拉普拉斯变换];
    F --> G[时域解y(t)];

二、拉普拉斯变换的具体应用示例

2.1 二阶非齐次微分方程的拉普拉斯变换

考虑一个二阶非齐次微分方程：
[
\frac{d^{2}y}{dt^{2}}+A\frac{dy}{dt}+C = e^{nt}, \quad \frac{dy(0)}{dt}=k, \quad y(0)=m
]
对该方程两边进行拉普拉斯变换：
- 已知(\mathcal{L}{e^{nt}}=\frac{1}{s - n})
- (\mathcal{L}\left{\frac{d^{2}y}{dt^{2}}\right}=s^{2}Y(s)-s y(0)-\frac{dy(0)}{dt}=s^{2}Y(s)-sm - k)
- (\mathcal{L}\left{A\frac{dy}{dt}\right}=A\left(sY(s)-y(0)\right)=A\left(sY(s)-m\right))
- (\mathcal{L}{C}=C)

将上述结果代入原方程的拉普拉斯变换式，得到：
[
s^{2}Y(s)-sm - k+A\left(sY(s)-m\right)+C=\frac{1}{s - n}
]
求解(Y(s))：
[
Y(s)=\frac{sm + k+Am - C+\frac{1}{s - n}}{s^{2}+As + C}=\frac{(sm + k+Am - C)(s - n)+1}{(s - n)(s^{2}+As + C)}
]
将(Y(s))分解为部分分式并进行逆拉普拉斯变换，可得到原微分方程的(y(t))：
[
y(t)=\frac{e^{nt}}{An - n^{2}}\left(Cn - kn - mn\right)+\frac{An - Cn+kn + mn - A(k + mn)+C}{A - n}e^{At}
]

2.2 MATLAB中的laplace()函数

MATLAB的符号数学工具箱提供了内置的 laplace() 函数，用于计算任何代数表达式或微分方程的拉普拉斯变换。其一般语法如下：

F = laplace(f);
F = laplace(f, t);
F = laplace(f, var1, var2);

以下是一些具体示例：
- 计算(x(t)=\sin(2t))的拉普拉斯变换：

syms t
xt = sin(2*t);
Xs = laplace(xt);
disp(Xs); % 输出: 2/(s^2 + 4)

• 计算(y(t)=\sin(Kt))的拉普拉斯变换：

syms t K
yt = sin(K*t);
Ys = laplace(yt);
disp(Ys); % 输出: K/(K^2 + s^2)

• 计算(y(x)=ax^{3}+b)的拉普拉斯变换：

syms x a b
yx = a*x^3 + b;
Ys = laplace(yx);
disp(Ys); % 输出: (6*a)/s^4 + b/s

也可以将变换结果表示为关于(t)的形式：

Yt = laplace(yx, x, t);
disp(Yt); % 输出: (6*a)/t^4 + b/t

2.3 ilaplace()函数及其应用

ilaplace() 函数用于计算拉普拉斯变换的逆变换，其语法与 laplace() 函数类似：

f = ilaplace(F);
f = ilaplace(F, s);
f = ilaplace(F, var1, var2);

下面通过几个ODE练习，展示如何使用 laplace/ilaplace 函数，并将其求解结果与 dsolve 函数的结果进行比较。

2.3.1 示例：比较laplace/ilaplace与dsolve求解一阶ODE

求解(\dot{y}+2y = 0, \quad y(0)=0.5)，以下是使用 laplace/ilaplace 和 dsolve 求解的脚本（ ODE_Laplace.m ）：

% ODE_Laplace.m
clearvars; clc; close all
% Step #1. Define symbolic variables' names
syms t s Y
ODE1 = 'D(y)(t)=-2*y(t)';
% Step #2. Laplace Transforms
LT_A = laplace(ODE1, t, s);
% Step #3. Substitute ICs and initiate an unknown Y
LT_A = subs(LT_A, {'laplace(y(t), t, s)', 'y(0)'}, {Y, 0.5});
% Step #4. Solve for Y (unknown)
Y = solve(LT_A, Y);
display('Laplace Transforms of the given ODE with ICs');
disp(Y);
% Step #5. Evaluate Inverse Laplace Transform
Solution_Laplace = ilaplace(Y);
display('Solution found using Laplace Transforms: ');
pretty(Solution_Laplace);
% Step #6. Compute numerical values and plot them
t = 0:0.01:2.5;
LTsol = eval(vectorize(Solution_Laplace));
figure, semilogx(t, LTsol, 'ro-');
xlabel('t');
ylabel('solution values');
title('laplace/ilaplace vs dsolve ');
grid on;
hold on;
% Compare with dsolve solution method
Y_d = dsolve('Dy=-2*y', 'y(0)=0.5');
display('Solution with dsolve');
pretty(Y_d);
Y_sol = eval(vectorize(Y_d));
plot(t, Y_sol, 'b-', 'linewidth', 2);
grid minor;
legend('laplace+ilaplace', 'dsolve');
hold off;
axis tight;

执行该脚本后，得到的输出如下：

Laplace Transforms of the given ODE with ICs
1/(2*s + 4)
Solution found using Laplace Transforms: 
  exp(-2 t)
  ---------
      2
Solution with dsolve
  exp(-2 t)
  ---------
      2

从结果可以看出， laplace/ilaplace 和 dsolve 得到的解是一致的。

2.3.2 示例：二阶非齐次ODE的求解

求解(\ddot{y}+\dot{y}=\sin(t), \quad y(0)=1, \quad \dot{y}(0)=2)，以下是求解脚本（ Laplace_vs_Dsolve.m ）：

% Laplace_vs_Dsolve.m
% %%   LAPLACE TRANSFORMS
% Given: y"+y'=sin(t) with ICs: [1, 2] for y(0) & y'(0).
clearvars, clc, close all
syms t s Y
ODE2nd = 'D(D(y))(t)+D(y)(t)-sin(t)';
% Step 1. Laplace Transforms
LT_A = laplace(ODE2nd, t, s);
% Step 2. Substitute ICs and initiate an unknown Y
LT_A = subs(LT_A, {'laplace(y(t), t, s)', 'y(0)', 'D(y)(0)'}, {Y, 1, 2});
% Step 3.  Solve for Y unknown
Y = solve(LT_A, Y);
display('Laplace Transforms of the given ODE with ICs');
disp(Y);
Solution_Laplace = ilaplace(Y);
display('Solution found using Laplace Transforms: ');
pretty(Solution_Laplace);
t = 0:0.01:13;
LTsol = eval(vectorize(Solution_Laplace));
figure, plot(t, LTsol, 'ro-');
xlabel('t');
ylabel('solution values');
title('laplace/ilaplace vs. dsolve: ddy+dy=sin(t)');
hold on;
% dsolve solution method
Y = dsolve('D2y+Dy=sin(t)', 'y(0)=1, Dy(0)=2', 't');
display('Solution with dsolve:  ');
pretty(Y);
Y_sol = eval(vectorize(Y));
plot(t, Y_sol, 'b-', 'linewidth', 2);
grid minor;
legend('laplace+ilaplace', 'dsolve');
hold off;

计算得到的解析解如下：

Laplace Transforms of the given ODE with ICs
(s + 1/(s^2 + 1) + 3)/(s^2 + s)
Solution found using Laplace Transforms: 
      cos(t)   sin(t)   5 exp(-t)
  4 - ------ - ------ - ---------
        2        2         2
Solution with dsolve: 
      cos(t)   sin(t)   5 exp(-t)
  4 - ------ - ------ - ---------
        2        2         2

从结果和绘制的图形可以看出， laplace/ilaplace 和 dsolve 得到的解完全一致。

2.4 无法得到解析解的情况

对于一些二阶非齐次和非线性ODE，可能无法使用 laplace/ilaplace 和 dsolve 得到解析解。例如，考虑方程(2\ddot{y}+3\dot{y}^{3}-\cos(100t)|y| - 2 = 0, \quad y(0)=0, \quad \dot{y}(0)=0)，以下是求解脚本（ Lap_inv_Lap.m ）：

% Lap_inv_Lap.m
clearvars, clc, close all
syms t s Y
ODE2nd = '2*D(D(y))(t)+3*(D(y)(t))^3-cos(100*t)*abs(y(t))-2';
% Step 1. Laplace Transforms
LT_A = laplace(ODE2nd, t, s);
% Step 2. Substitute ICs and initiate an unknown Y
LT_A = subs(LT_A, {'laplace(y(t), t, s)', 'y(0)', 'D(y)(0)'}, {Y, 0, 0});
% Step 3.  Solve for Y unknown
Y = solve(LT_A, Y);
display('Laplace Transforms of the given ODE with ICs');
disp(Y); 
Solution_Laplace = ilaplace(Y);
display('Solution found using Laplace Transforms: ');
pretty(Solution_Laplace);

执行该脚本后，输出结果表明无法得到明确的解析解。

2.5 拉普拉斯变换在处理不连续强迫函数的ODE中的应用

考虑一个二阶非齐次ODE：(2\ddot{y}+3\dot{y}-2y = g(t), \quad y(0)=0, \quad \dot{y}(0)=0)，其中不连续强迫函数(g(t))定义为：
[
g(t)=\begin{cases}
2, & 0\leq t<2 \
5, & 2\leq t<10 \
0, & t\geq 10
\end{cases}
]
该方程的拉普拉斯变换为：
[
2\left(s^{2}Y(s)-sy(0)-\dot{y}(0)\right)+3\left(sY(s)-y(0)\right)-2Y(s)=\frac{5\left(e^{-2s}-e^{-10s}\right)}{s}
]
求解(Y(s))：
[
Y(s)=\frac{5\left(e^{-2s}-e^{-10s}\right)}{s\left(2s^{2}+3s - 2\right)}
]
使用 ilaplace() 函数计算逆变换，以下是求解脚本（ Lap_4_non_homog.m ）：

% Lap_4_non_homog.m
syms t s
F = 5*(exp(-2*s)-exp(-10*s))/s;
Y = 2*s^2 + s + 2;
TF = F/Y;
TFt = ilaplace(TF);
pretty(TFt);
Sol = vectorize(TFt);
t = linspace(0, 20, 400);
S = eval(Sol);
plot(t, S, 'bo-');
grid minor;
title('Differential Equation with Discontinuous Forcing Fcn');
grid on;
xlabel('time');
ylabel('y(t) solution');
shg;

执行该脚本后，得到的解和绘制的图形表明，使用拉普拉斯变换（ laplace/ilaplace ）计算具有不连续外部强迫函数的非齐次ODE的解析解相对容易、快速且轻松。

综上所述，拉普拉斯变换和逆变换在求解ODE中具有重要作用，特别是在处理具有不连续强迫函数的非齐次ODE时表现出明显的优势。同时， laplace/ilaplace 函数和 dsolve 函数在大多数情况下可以得到一致的解析解，但对于一些复杂的ODE，可能无法得到解析解，此时需要使用数值方法进行求解。

三、拉普拉斯变换与逆变换的优势总结

3.1 处理不连续强迫函数的便利性

在处理具有不连续强迫函数的非齐次常微分方程时，拉普拉斯变换展现出了显著的优势。如前面提到的二阶非齐次 ODE (2\ddot{y}+3\dot{y}-2y = g(t))，其中 (g(t)) 是不连续的分段函数。通过拉普拉斯变换，将时域的不连续问题转换为 s 域的代数问题，避免了在时域中直接处理不连续函数的复杂性。具体操作步骤如下：
1. 对原方程进行拉普拉斯变换，根据拉普拉斯变换的性质，将方程中的各项转换为 s 域的表达式。
2. 代入初始条件，将方程化简为关于 (Y(s)) 的代数方程。
3. 求解 (Y(s))，得到 s 域的解。
4. 对 (Y(s)) 进行逆拉普拉斯变换，得到时域的解 (y(t))。

在 MATLAB 中，可以使用 laplace() 和 ilaplace() 函数方便地完成这些操作，如 Lap_4_non_homog.m 脚本所示。这种方法相对容易、快速且轻松，能够有效地处理不连续强迫函数的问题。

3.2 与 dsolve 函数的一致性

通过多个示例可以看出， laplace/ilaplace 函数和 dsolve 函数在大多数情况下可以得到一致的解析解。例如，在求解一阶 ODE (\dot{y}+2y = 0, \quad y(0)=0.5) 和二阶非齐次 ODE (\ddot{y}+\dot{y}=\sin(t), \quad y(0)=1, \quad \dot{y}(0)=2) 时，两种方法得到的解完全相同。这表明拉普拉斯变换和逆变换是一种可靠的求解 ODE 的方法，与传统的解析求解方法（如 dsolve 函数）具有一致性。

3.3 操作的简便性

MATLAB 的符号数学工具箱提供的 laplace() 和 ilaplace() 函数，使得拉普拉斯变换和逆变换的操作变得非常简便。只需要定义符号变量，使用相应的函数进行变换和逆变换，就可以得到所需的结果。例如，计算 (x(t)=\sin(2t)) 的拉普拉斯变换，只需要以下几行代码：

syms t
xt = sin(2*t);
Xs = laplace(xt);
disp(Xs);

这种简便的操作方式，降低了使用拉普拉斯变换求解 ODE 的门槛，使得工程师和科研人员可以更方便地应用这种方法解决实际问题。

四、拉普拉斯变换与逆变换的应用场景

4.1 电路分析

在电路分析中，常微分方程经常用于描述电路中的电流和电压变化。例如，对于一个包含电阻、电感和电容的 RLC 电路，其电压和电流的关系可以用二阶非齐次常微分方程来描述。拉普拉斯变换可以将这些时域的电路方程转换为 s 域的代数方程，从而方便地求解电路的响应。具体步骤如下：
1. 根据电路的基本定律（如基尔霍夫定律），列出电路的时域方程。
2. 对时域方程进行拉普拉斯变换，将其转换为 s 域的代数方程。
3. 求解 s 域的代数方程，得到电压或电流的 s 域表达式。
4. 对 s 域表达式进行逆拉普拉斯变换，得到时域的电压或电流响应。

4.2 控制系统设计

在控制系统设计中，常微分方程用于描述系统的动态行为。拉普拉斯变换可以将系统的时域模型转换为 s 域的传递函数模型，从而方便地进行系统的分析和设计。例如，对于一个机械弹簧 - 质量 - 阻尼系统，其运动方程可以用二阶非齐次常微分方程来描述。通过拉普拉斯变换，可以得到系统的传递函数，进而分析系统的稳定性、响应特性等。具体步骤如下：
1. 根据系统的物理原理，列出系统的时域方程。
2. 对时域方程进行拉普拉斯变换，将其转换为 s 域的传递函数。
3. 根据传递函数，分析系统的稳定性、响应特性等。
4. 根据分析结果，设计控制器，对系统进行优化。

4.3 信号处理

在信号处理中，拉普拉斯变换可以用于分析和处理信号的频谱特性。例如，对于一个非周期信号，可以通过拉普拉斯变换将其转换为 s 域的表达式，从而分析信号的频率成分和衰减特性。具体步骤如下：
1. 定义信号的时域表达式。
2. 对时域信号进行拉普拉斯变换，得到信号的 s 域表达式。
3. 分析 s 域表达式，得到信号的频率成分和衰减特性。
4. 根据分析结果，对信号进行滤波、调制等处理。

五、总结与展望

5.1 总结

拉普拉斯变换和逆变换是求解常微分方程的重要工具，具有以下优点：
- 能够处理具有不连续强迫函数的非齐次常微分方程，避免了在时域中直接处理不连续函数的复杂性。
- 与 dsolve 函数在大多数情况下可以得到一致的解析解，表明其可靠性。
- 操作简便，通过 MATLAB 的 laplace() 和 ilaplace() 函数可以方便地完成变换和逆变换。

5.2 展望

虽然拉普拉斯变换和逆变换在求解常微分方程中具有重要作用，但对于一些复杂的 ODE，可能无法得到解析解。在这种情况下，需要使用数值方法进行求解。未来的研究可以进一步探索拉普拉斯变换和逆变换与数值方法的结合，提高求解复杂 ODE 的能力。此外，随着计算机技术的不断发展，拉普拉斯变换和逆变换在更多领域的应用也将得到进一步拓展，如生物医学工程、航空航天等领域。

同时，为了更好地应用拉普拉斯变换和逆变换，建议读者在实际操作中多进行练习，熟悉 MATLAB 的 laplace() 和 ilaplace() 函数的使用方法。可以通过以下步骤进行练习：
1. 选择一些简单的常微分方程，使用拉普拉斯变换和逆变换进行求解，并与 dsolve 函数的结果进行比较。
2. 尝试处理具有不连续强迫函数的非齐次常微分方程，掌握拉普拉斯变换在处理这类问题中的优势。
3. 结合实际应用场景，如电路分析、控制系统设计等，将拉普拉斯变换和逆变换应用到实际问题中，提高解决实际问题的能力。

总之，拉普拉斯变换和逆变换是一种强大的数学工具，在常微分方程求解和实际应用中具有重要的价值。通过不断学习和实践，读者可以更好地掌握这一工具，为解决实际问题提供有力的支持。