传输线方程是一阶双曲型偏微分方程，求解此类型方程没有现成的差分格式，因此需要把传输线方程转换成一阶拟线性偏微分方程组，然后采用Lax差分格式，偏心格式或者Lax-Wendroff格式对方程组进行求解，此种方法计算量比较大，且转换为一阶拟线性偏微分方程组的推导过程比较复杂^[1-2]。本文采用常微分方程数值解法求解传输线的电磁瞬态过程，求解分为两个步骤，一是将传输线方程转换为常微分方程组，二是在初始条件下求解常微分方程组。求解常微分方程组的数值解法有很多，常用的有Euler方法、显示单步法、龙格−库塔（Runge-Kutta，RK）方法、线性多步法等^[1-6]。为了提高计算精度与计算效率，本文将RK方法与空间高阶展开结合应用至传输线方程的求解当中，并扩展了此种算法的应用范围。

1 高精度龙格−库塔方法理论

1.1 改进RK方法的迭代方程

无源的传输线方程如下

1$ \left\{ \begin{array}{l} - \dfrac{{\partial u}}{{\partial {\textit{z}}}} = {R_0}i + {L_0}\dfrac{{\partial i}}{{\partial t}} \\ - \dfrac{{\partial i}}{{\partial {\textit{z}}}} = {G_0}u + {C_0}\dfrac{{\partial u}}{{\partial t}} \end{array} \right. $

两端的边界条件如下

2$ {u_0} = - {i_0}{R_s} $

3$ {u_l} = {i_l}{R_l} $

式中： $u$为线上电压； ${i}$为线上电流； ${R_0}$为单位长度电阻； ${L_0}$为单位长度电感； ${G_0}$为单位长度电导； ${C_0}$为单位长度电容。

将传输线分为N段，每一个节点上都有包含电压值和电流值，整个方程组共有2N个方程，含有 $2N + 2$个变量，添加2个边界条件后，只需求解2N个独立变量。

为了得到空间四阶中心差分近似，采用泰勒级数展开进行处理。

4$ \begin{split} & f(x + \dfrac{1}{2}\Delta x) = f\left( x \right) + \left( { \pm \dfrac{1}{2}\Delta x} \right)\dfrac{\partial }{{\partial x}}f\left( x \right) + \dfrac{1}{{2!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^2}\dfrac{{{\partial ^2}}}{{\partial {x^2}}}f\left( x \right) + \\ & \dfrac{1}{{3!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^3}\dfrac{{{\partial ^3}}}{{\partial {x^3}}}f\left( x \right) + \dfrac{1}{{4!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^4}\dfrac{{{\partial ^4}}}{{\partial {x^4}}}f\left( x \right) + \dfrac{1}{{5!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^5}\dfrac{{{\partial ^5}}}{{\partial {x^5}}}f\left( x \right) + \cdots \end{split} $

5$ \begin{split} & f(x + \dfrac{3}{2}\Delta x) = f\left( x \right) + \left( { \pm \dfrac{3}{2}\Delta x} \right)\dfrac{\partial }{{\partial x}}f\left( x \right) + \dfrac{1}{{2!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^2}\dfrac{{{\partial ^2}}}{{\partial {x^2}}}f\left( x \right) + \\ & \dfrac{1}{{3!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^3}\dfrac{{{\partial ^3}}}{{\partial {x^3}}}f\left( x \right) + \dfrac{1}{{4!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^4}\dfrac{{{\partial ^4}}}{{\partial {x^4}}}f\left( x \right) + \dfrac{1}{{5!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^5}\dfrac{{{\partial ^5}}}{{\partial {x^5}}}f\left( x \right) + \cdots \end{split} $

根据式（4）和式（5）可以得到多项式如下

6$ 27\left[ {f\Bigg(x + \dfrac{1}{2}\Delta x\Bigg) - f\Bigg(x - \dfrac{1}{2}\Delta x\Bigg)} \right] - \left[ {f\Bigg(x + \dfrac{3}{2}\Delta x\Bigg) - f\Bigg(x - \dfrac{3}{2}\Delta x\Bigg)} \right] $

采用式（6）代替式（1）中的空间微分，得到空间四阶中心差分后的方程

7$ \left\{ \begin{array}{l} \dfrac{{\partial {u_k}}}{{\partial t}} \approx \dfrac{{\left( {27{i_{k - 1}} + {i_{k + 1}} - 27{i_k} - {i_{k - 2}}} \right)}}{{24\Delta {\textit{z}}{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_k} \\ \dfrac{{\partial {i_k}}}{{\partial t}} \approx \dfrac{{\left( {27{u_k} + {u_{k + 2}} - 27{u_{k + 1}} - {u_{k - 1}}} \right)}}{{24\Delta {\textit{z}}{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_k} \end{array} \right. $

式中： ${u_k} = {u_1},{u_2} \ldots {u_N}$， ${i_k} = {i_0},{i_1} \ldots {i_{N - 1}}$； $N = l/\Delta {\textit{z}}$。

由于在空间高阶展开时，求一个点的电压值和电流值需要邻近多个网格点的值，因此在边界以及边界附近采用二阶的展开，对式（1）在 $\left( {k + 1} \right)\Delta {\textit{z}}$处采取后向欧拉公式进行离散，对式（1）在 $k\Delta {\textit{z}}$处采取前向欧拉公式进行离散得到

8$ \left\{ \begin{array}{l} \dfrac{{\partial {u_{k + 1}}}}{{\partial t}} \approx \dfrac{{{i_k} - {i_{k + 1}}}}{{\Delta {\textit{z}}{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_{k + 1}} \\ \dfrac{{\partial {i_k}}}{{\partial t}} \approx \dfrac{{{u_k} - {u_{k + 1}}}}{{\Delta {\textit{z}}{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_k} \end{array} \right. $

将边界条件代入式（8）可以得到边界上的迭代公式

9$ \left\{ \begin{array}{l} \dfrac{{\partial {u_N}}}{{\partial t}} \approx \dfrac{{{i_{N{\rm{ - 1}}}} - {R_l}^{ - 1}{u_N}}}{{\Delta {\textit{z}}{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_N} \\ \dfrac{{\partial {i_0}}}{{\partial t}} \approx \dfrac{{ - {R_{\rm{s}}}{i_0} - {u_N}}}{{\Delta {\textit{z}}{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_0} \end{array} \right. $

式中： ${R_{\rm{s}}}$， ${R_l}$分别为源端和终端电阻。

将电报方程的空间坐标离散后总共得到2N个相互独立的状态变量，并得到2N阶的线性常态状态方程组：

10$ \frac{{{\rm{d}}{{X}}}}{{{\rm{d}}t}} = {{FX}} + {{s}}(t) $

其中 ${{s}}\left( t \right)$是离散后的源项， ${{X}} = \left[ \begin{array}{l} u \\ i \\ \end{array} \right]$，F为常数矩阵。

11$ {{F}} = \left[\!\!\!\!\!\! {\begin{array}{*{20}{c}} {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\! \ddots \!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24C\Delta {\textit{z}}}}} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{C\Delta {\textit{z}}}}} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{{R_l}C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}} \\ {\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - {R_s}}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {\dfrac{1}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {\dfrac{{ - 1}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\! \ddots \!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \end{array}} \!\!\!\!\!\!\right] $

4级4阶的RK公式如下

12$ \left\{ \begin{array}{l} {y_{n + 1}} = {y_n} + \dfrac{1}{6}h\left( {{k_1} + 2{k_2} + 2{k_3} + {k_4}} \right) \\ {k_1} = f\left( {{x_n},{y_n}} \right) \\ {k_2} = f\left( {{x_n} + \dfrac{1}{2}h,{y_n} + \dfrac{1}{2}h{k_1}} \right) \\ {k_3} = f\left( {{x_n} + \dfrac{1}{2}h,{y_n} + \dfrac{1}{2}h{k_2}} \right) \\ {k_4} = f\left( {{x_n} + h,{y_n} + h{k_3}} \right) \end{array} \right. $

其中， ${k_1},{k_2},{k_3},{k_4}$为RK方法的系数，h为步长， $x \in \left[ {{x_0},b} \right]$，b可以是有限的，也可以是 $ + \infty $，y为 $\left[ {{x_0},b} \right]$到R的连续可微函数。

考虑外界平面波入射，由于RK方法中电流点电压点在同一空间位置，可以得到源项如下

13$ \left\{ \begin{array}{l} - \dfrac{\partial }{{\partial {\textit{z}}}}{E_{\rm{T}}}\left( {{\textit{z}},t} \right) + {E_{\rm{L}}}\left( {{\textit{z}},t} \right) = (\dfrac{{{A_{\rm{T}}}}}{{{v_{\textit{z}}}}} - \dfrac{{{A_{\rm{L}}}}}{{{v_x}}})\left[ {{\xi _0}\left( {{t^{n + 3/2}} - \dfrac{{k\Delta {\textit{z}}}}{{{v_{\textit{z}}}}}} \right) - {\xi _0}\left( {{t^{n + 1/2}} - \dfrac{{k\Delta {\textit{z}}}}{{{v_{\textit{z}}}}}} \right)} \right]/\Delta t \\ - C\dfrac{\partial }{{\partial {\textit{z}}}}{E_{\rm{T}}}\left( {{\textit{z}},t} \right) = - C{A_{\rm{T}}}\left[ {{\xi _0}\left( {{t^{n + 3/2}} - \dfrac{{k\Delta {\textit{z}}}}{{{v_{\textit{z}}}}}} \right) - {\xi _0}\left( {{t^{n + 1/2}} - \dfrac{{k\Delta {\textit{z}}}}{{{v_{\textit{z}}}}}} \right)} \right]/\Delta t \end{array} \right. $

其中， ${E_{\rm{T}}}\left( {{\textit{z}},t} \right)$和 ${E_{\rm{L}}}\left( {{\textit{z}},t} \right)$分别为入射场垂直和平行于传输线导体的分量， ${A_{\rm{T}}}$， ${A_{\rm{L}}}$计算如下

14$ \left\{ \begin{array}{l} {A_{\rm{T}}} = \left( {{e_x}{x_k} + {e_y}{y_k}} \right) \\ {A_{\rm{L}}} = \left( {\dfrac{{{x_k}}}{{{v_x}}} + \dfrac{{{y_k}}}{{{v_y}}}} \right){e_{\textit{z}}} \end{array} \right. $

${e_x}$， ${e_y}$， ${e_{\textit{z}}}$为方向矢量， $v_x,{v_y},{v_{\textit{z}}}$为波沿坐标轴的传播速度

15$ \left\{ \begin{array}{l} {e_x} = \sin {\theta _{\rm{E}}}\sin {\theta _{\rm{p}}} \\ {e_y} = - \sin {\theta _{\rm{E}}}\cos {\theta _{\rm{p}}}\cos {\phi _{\rm{p}}} - \cos {\theta _{\rm{E}}}\sin {\phi _{\rm{p}}} \\ {e_z} = - \sin {\theta _{\rm{E}}}\cos {\theta _{\rm{p}}}\sin {\phi _{\rm{p}}} + \cos {\theta _{\rm{E}}}\cos {\phi _{\rm{p}}} \end{array} \right. $

16$ \left\{ \begin{array}{l} v{}_x = - \dfrac{v}{{\cos {\theta _{\rm{p}}}}} \\ {v_y} = - \dfrac{v}{{\sin {\theta _{\rm{p}}}\cos {\phi _{\rm{p}}}}} \\ {v_{\textit{z}}} = - \dfrac{v}{{\sin {\theta _{\rm{p}}}\sin {\phi _{\rm{p}}}}} \end{array} \right. $

其中， ${\theta _{\rm{p}}},{\varphi _{\rm{p}}},{\theta _{\rm{E}}}$分别为入射波的入射角、方位角、极化角。

1.2 收敛性分析

根据文献[7]给出的定理，假设初值问题的单步方法的增量函数 $\varphi \left( {x,y;h} \right)$关于x，y和h是连续的，并对于y满足Lipschitz条件，即 $\left| {\varphi \left( {x,y;h} \right) - \varphi \left( {x,{\textit{z}};h} \right)} \right| \leqslant L\left| {y - {\textit{z}}} \right|$，那么单步法收敛的充分必要条件是单步法是相容的。为了证明显示RK方法的收敛性，只要证明其相容性以及增量函数满足Lipschitz条件就可以了。

以经典的2级2阶RK方法（中点公式）为例

17$ \varphi \left( {x,y;h} \right) = f\left( {x + \dfrac{1}{2}h,y + \dfrac{1}{2}hf(x,y)} \right) $

18$ \mathop {\lim }\limits_{h \to 0} \varphi \left( {x,y;h} \right) = f\left( {x,y} \right) $

上式表明中点公式满足相容性。

19$ \begin{split} & {\rm{ }}\left| {\varphi \left( {x,{y_1};h} \right) - \varphi \left( {x,{y_2};h} \right)} \right| = \left| {f\Bigg(x + \dfrac{h}{2},{y_1} + \dfrac{1}{2}hf(x,{y_1})\Bigg) - f\Bigg(x + \dfrac{h}{2},{y_2} + \dfrac{1}{2}hf(x,{y_2})\Bigg)} \right| \leqslant \\ & {L_1}\left| {{y_1} - {y_2} + \dfrac{1}{2}h\left[ {f(x,{y_1}) - f(x,{y_2})} \right]} \right| \leqslant {L_1}\left| {{y_1} - {y_2}} \right| + \dfrac{1}{2}h{L_1} \cdot {L_1}\left| {{y_1} - {y_2}} \right| = {L_1}\left( {1 + \dfrac{1}{2}h{L_1}} \right)\left| {{y_1} - {y_2}} \right| \end{split} $

当 $h \leqslant {h_0}$，取 $L = {L_1}\left( {1 + \frac{1}{2}{h_0}{L_1}} \right)$。此时 $\varphi \left( {x,y;h} \right)$对y满足Lipschitz条件，从而可知中点公式收敛，其他RK方法皆可根据此方法证明收敛。

1.3 稳定性分析

为了简单起见，一般把数值方法用于试验方程来进行讨论^[8]。

20$ y'\left( x \right) = \lambda y\left( x \right) $

其中， ${\rm{Re}} \left( \lambda \right) < 0$，若 $\lambda $为实数，则 $\lambda < 0$。

将经典的4级4阶龙格—库塔方法应用至试验方程得到

21$ {y_{n + 1}} = \left[ {1 + \lambda h + \frac{{{{\left( {\lambda h} \right)}^2}}}{2} + \frac{{{{\left( {\lambda h} \right)}^3}}}{6} + \frac{{{{\left( {\lambda h} \right)}^4}}}{{24}}} \right]{y_n} $

22$ E\left( {\lambda h} \right) = {1 + \lambda h + \frac{{{{\left( {\lambda h} \right)}^2}}}{2} + \frac{{{{\left( {\lambda h} \right)}^3}}}{6} + \frac{{{{\left( {\lambda h} \right)}^4}}}{{24}}} $

如果 $\left| {E\left( {\lambda h} \right)} \right| < 1$，那么 $\lambda h$的区域为经典的4级4阶龙格库塔方法的绝对稳定性区域，见图1。其中横轴为实数轴，竖轴为虚数轴。

图 1. Stability region of the Runge-Kutta method

Fig. 1. Stability region of the Runge-Kutta method

下载图片 查看所有图片

1.4 计算内存

RK数值解法需要存储一个系数矩阵F，以及4个中间变量 ${k_1},{k_2},{k_3},{k_4}$。为了尽量减小计算内存，采用稀疏矩阵的存储方式对F进行存储，每得到一个中间变量就把它放入累加器中，可以得到

23$ {\rm{R}}{{\rm{K}}_4}\text{-}{\rm{H}}{{\rm{O}}_4}\text{数值解法所需内存} = \left[ {\left( {2N} \right) \times 12 + 2N + 2N} \right] \times 4 $

2 龙格−库塔数值方法有效性验证及算例

2.1 双导体传输线

采用经典算例进行计算^[8-10]，一根半径 ${r_{\rm{w}}} = 0.254\;{\rm{ mm}}$、长 $L = 1\;{\rm{ m}}$的导线位于无限大接地平面以上，距地高度 $h = 2\;{\rm{ cm}}$。终端的阻性负载分别为 ${R_{\rm{S}}} = 500\;{\rm{ \Omega }}$和 ${R_{\rm{L}}} = 1\;000\;{\rm{ \Omega }}$。入射的均匀平面波从顶部入射， ${\theta _{\rm{E}}} = 0^\circ $， ${\theta _{\rm{p}}} = 0^\circ $， ${\phi _{\rm{p}}} = 0^\circ $。入射电场 $\xi \left( t \right)$的时域波形为具有不同上升时间的梯形周期脉冲序列，重复频率 $1\;{\rm{ MHz}}$，幅值 $1\;{\rm{V/m}}$，占空比50%。脉冲上升时间设置为 $10\;{\rm{ ns}}$，如图2所示。

图 2. Transmission line model

Fig. 2. Transmission line model

下载图片 查看所有图片

计算结果如图3所示，可以看出当空间步长为 $0.02\;{\rm{m}}$时，RK数值方法的计算结果与FDTD方法的计算结果在此算例上相同，验证了本文提出的RK方法的正确性。当空间步长为 $0.02\;{\rm{m}}$时，此时RK数值方法依旧稳定，但是由于步长过大导致了计算结果误差较大。

图 3. Voltage at left end of twin conductor transmission line

Fig. 3. Voltage at left end of twin conductor transmission line

下载图片 查看所有图片

2.2 多导体传输线

针对多导体平行均匀传输线，设置以下算例，如图4（a）和图4（b）所示。导线长2 m，导线半径为0.190 5 mm，绝缘层厚度 ${r_{\rm{w}}} = 0.254\;{\rm{mm}}$，导线中心之间的距离为1.3 mm。均匀平面波沿着z方向入射，电场的极化方向沿x方向。单位长度的参数矩阵 $L = \left[\! \begin{array}{l} 0.748\;5\quad {\rm{0}}{\rm{.240\;8}} \\ {\rm{0}}{\rm{.240\;8}}\quad 0{\rm{.748\;5}} \\ \end{array} \!\right]{\rm{ \mu H/m}}$， $C = \left[\! \begin{array}{l} 24.982\quad {\rm{ - 6}}{\rm{.266}} \\ {\rm{ - 6}}{\rm{.266}}\quad 24{\rm{.982}} \\ \end{array} \!\right]{\rm{ pF/m}}$。传输线终端接 $500\;{\rm{\Omega }}$电阻， ${R_{\rm{S}}} = {R_{\rm{L}}} = \left[\! \begin{array}{l} {\rm{500}} \quad 0 \\ {\rm{0 }} \quad 500 \\ \end{array} \!\right]{\rm{ }}\Omega $，入射电场 $\xi \left( t \right)$的时域波形如图4（c）所示。入射波入射角度分别为 ${\theta _{\rm{E}}} = {90^ \circ },{\theta _{\rm{p}}} = {90^ \circ },$${\phi _{\rm{p}}} = - {90^ \circ }$。

图 4. Model of multiconductor transmission lines

Fig. 4. Model of multiconductor transmission lines

下载图片 查看所有图片

计算结果如图5所示，结果验证了改进RK方法针对场激励多导体传输线计算的正确性。

图 5. Voltage at left end of multiconductor transmission line

Fig. 5. Voltage at left end of multiconductor transmission line

下载图片 查看所有图片

2.3 RK方法与FDTD方法的计算精度

针对集中激励源激励的传输线进行计算，采用图6（a）所示的传输线模型，采用图6（b）所示的电压源，线长 $0.8\;{\rm{ m}}$， ${L_0} = 309\;{\rm{ nH/m}}$， ${C_0} = 144\;{\rm{ pF/m}}$, ${R_{\rm{s}}} = 50\;{\rm{ \Omega }},{R_{\rm{L}}} = 50\;{\rm{ \Omega }}$。

图 6. Lossless transmission line model

Fig. 6. Lossless transmission line model

下载图片 查看所有图片

根据传输线基本理论

24$ {Z_0} = \sqrt {{L_0}/{C_0}} $

25$ {u_{{R_{\rm{s}}}}}\left( t \right) = \frac{{{Z_0}}}{{{R_{\rm{s}}} + {Z_0}}}{V_{\rm{s}}}\left( t \right) $

26$ \beta = \frac{{{u_{\rm{L}}}}}{{{u_{{R_{\rm{s}}}}}}} = \frac{{2{R_{\rm{L}}}}}{{{R_{\rm{L}}} + {Z_0}}} $

其中， ${Z_0}$为传输线的特征阻抗， ${u_{{R_{\rm{s}}}}}\left( t \right)$为激励源加载至传输线上的电压， $\;\beta $为透射系数， ${u_{\rm{L}}}$为负载电压，对式（18）进行求解可以得到负载电压（计算过程中保留有效小数点后8位）

27$ {u_{\rm{L}}}\left( t \right) = \frac{{{Z_0}}}{{{R_{\rm{s}}} + {Z_0}}}{V_{\rm{s}}}\left( t \right)\frac{{2{R_{\rm{L}}}}}{{{R_{\rm{L}}} + {Z_0}}} $

根据式（20）可以画出负载电压波形图，如图7所示。

图 7. Waveform of load

Fig. 7. Waveform of load

下载图片 查看所有图片

FDTD方法选取 $\Delta {\textit{z}} = 0.000\;8\;{\rm{ m}}$， $\Delta t = 5 \times {10^{ - 12}}\;{\rm{ s}}$，RK数值方法选取 $\Delta {\textit{z}} = 0.005\;{\rm{ m}}$， $\Delta t = 1 \times {10^{ - 11}}\;{\rm{ s}}$，计算结果如图8所示。

图 8. Terminal voltage response

Fig. 8. Terminal voltage response

下载图片 查看所有图片

RK方法的计算时间为1.1 s，FDTD的计算时间为4.8 s。从图8可以看出，RK方法与FDTD方法在波形转换的前沿处都存在振荡现象，其中RK方法的振荡现象更为严重，这是由于RK方法的空间步长比FDTD方法的空间步长更大，且RK方法对于激励源的光滑性有更高的要求。从波形平稳处截取10 ns点的波形与解析解进行对比，RK方法与解析解的误差为0.000 449 79，FDTD方法与解析解的误差为0.016 700 00。

上述算例的计算结果表明，当激励源的函数足够光滑，在同等的计算时间下，RK方法的计算精度更高。

3 结　论

为了改善传输线方程的求解精度，本文将RK方法引入至外界场激励下传输线方程的求解。提高了空间展开阶数，扩展了RK方法应用范围，并对其收敛特性与绝对稳定性做了分析。理论分析表明，激励源的函数足够光滑时，在相同的计算条件下，此方法的计算精度高于传统FDTD方法，通过数值算例验证，同等精度下，改进RK方法消耗的计算时间少于传统FDTD方法。