基于Newton迭代动态调整的Runge-Kutta法的构造及其在Blasius方程中的应用

Construction of the Runge-Kutta Method Based on Newton’s Iterative Dynamic Adjustment and Its Application to the Blasius Equation

摘要: 本文提出了一种新型数值求解方法，该方法将四阶Runge-Kutta法与Newton迭代法相结合，旨在高效求解流体力学中的Blasius方程边值问题。首先，我们将Blasius方程转化为一组一阶微分方程组，并采用四阶Runge-Kutta法(RK4)进行数值求解。随后，引入Newton迭代法动态调整初始条件，以确保满足边界条件的要求。实验结果表明，与传统的打靶法和结合打靶法的四阶Runge-Kutta法(SRK)进行对比实验，新方法在迭代次数和计算时间上均展现显著优势，同时求解精度也得到提升。

Abstract: In this paper, we propose a novel numerical solution method that combines the fourth-order Runge-Kutta method with the Newton iterative method, aiming to efficiently solve the margin problem of Blasius equation in fluid mechanics. Firstly, we transform the Blasius equation into a set of first-order differential equations and solve it numerically using the fourth-order Runge-Kutta method (RK4). Subsequently, the Newton iteration method is dynamically introduced to adjust the initial conditions to ensure that the boundary conditions are satisfied. The experimental results show that the new method exhibits significant advantages in terms of the number of iterations and computation time, as well as improved solution accuracy, in comparison experiments with the traditional shooting method and the Runge-Kutta method (SRK) combined with the improved shooting method.

1. 引言

在流体力学领域，边界层现象对于深入理解流体在固体表面附近的流动特性具有至关重要的意义。作为描述无粘流体在平板上稳态边界层流动的核心方程，Blasius方程[1]在该领域占据了不可替代的重要地位。尽管Blasius方程的解已广为人知，但在实际工程应用中，复杂的边界条件和多变的流体特性常常使得数值求解方法成为更为实用且必要的选择。

传统的数值求解方法，如打靶法、有限差分法、有限元法和谱方法等，各自具有独特的优势和局限性。打靶法结构简单，但强烈依赖于初始猜测解。有限差分法虽然实现简单且易于理解，但在处理非线性方程和高阶边界条件时，常常面临收敛速率慢和稳定性不足的挑战。有限元法则以其强大的解的逼近能力而著称，但其实现过程相对复杂，且对参数调整的要求较高。谱方法虽然在某些特定情况下能提供高效的精度，但其高度依赖于高频率的基函数，这在一定程度上限制了其适用范围。

为了克服这些传统方法的局限性，研究者们一直在不断探索新的数值求解策略。在打靶法的发展历程中，该方法最初被广泛应用于解决各种边值问题，通过不断调整初始条件来逼近满足边界条件的解[2]。随着计算技术的不断进步，打靶法也逐渐得到了改进和优化。多重打靶法(Multiple Shooting Method)作为打靶法的一种重要变体，在1984年由Bock和Plitt [3]首次提出并应用于优化控制问题的求解。该方法通过将整个求解区间划分为多个子区间，并在每个子区间上分别应用打靶法，从而显著提高了计算效率和收敛性。此后，多重打靶法在数值求解领域得到了广泛应用，并不断发展和完善[4]-[6]。

与此同时，Runge-Kutta法[7]作为一种常用的求解常微分方程的数值方法，也在不断改进和优化中。该方法基于离散化的思想，通过逐步逼近的方式计算方程的数值解，具备收敛速度快、精度高等优点，尤其适用于求解非线性方程和高阶方程[8]。近年来，Sadasivuni等人[9]提出了一种创新的数值求解方法，将Runge-Kutta法与打靶法相结合，利用打靶法确定边值问题的初始条件，并借助Runge-Kutta法求解相应的常微分方程，实现了更快的收敛性与更高的求解精度。

尽管这些方法取得了显著成果，但在调整初始条件时仍需进行多次迭代，并且Runge-Kutta法在处理某些复杂问题时仍存在一定的局限性。为此，本研究提出了一种结合四阶Runge-Kutta法与Newton迭代法的新型数值求解方法，旨在高效求解Blasius方程的边值问题。这种方法不仅提高了收敛效率，还确保了求解精度，对工程应用和科学计算具有重要的实用价值。

2. 问题描述

Blasius方程是流体力学中描述无粘性流体在平板上流动的边界层方程，其形式如下：

$f^{‴}+\frac{1}{2}f{f}^{″}=0$ ， (1)

其中，$f\left(\eta \right)$ 表示流函数，$\eta$ 是垂直于平板的坐标。该方程的边界条件为：

$f\left(0\right)=0$ , $f^{\prime }\left(0\right)=0$ , $f^{\prime }\left(\infty \right)=1$ . (2)

这些边界条件具有明确的物理意义：

$f\left(0\right)=0$ 表示在平板表面流体的流速为零，这是由于流体与固体表面之间的粘性作用导致的流体粘附现象。

$f^{\prime }\left(0\right)=0$ 表示在平板表面流体的流函数为零，这是无滑移条件的直接体现。

$f^{\prime }\left(\infty \right)=1$ 表示在边界层外，流体的流速梯度趋向于自由流速，即远离固体表面的区域，流体的流动不再受固体表面的影响。

尽管Blasius方程的近似解已为人所熟知，但在实际工程应用中，由于边界条件的复杂性和流体特性的多变性，数值求解方法成为更为实用且必要的选择。本研究的目标是通过数值方法求解Blasius方程，以期为流体力学领域的研究和工程应用提供精确的数值解。

3. 方法

本节将详细介绍用于求解Blasius方程的结合打靶法的四阶Runge-Kutta法及结合Newton迭代的新型Runge-Kutta法。

3.1. 结合打靶法的四阶Runge-Kutta法

四阶Runge-Kutta法(RK4)是一种广泛用于求解常微分方程初值问题的数值方法。该方法通过在每个步长内计算多个中间斜率来近似下一个点的值，具有较高的精度和稳定性，特别适合处理Blasius方程这类非线性问题。

3.1.1. 方法步骤

考虑微分方程的通用形式：

$y^{\prime }=f\left(x,y\right)$ . (3)

四阶Runge-Kutta法的步骤如下：

• 初始化：给定初始条件$x_0$ 和$y_0$ ，设定步长h。

• 迭代计算：对于每一步计算，使用以下公式：

$k_1=f\left(x_n,y_n\right)$ ,

$k_2=f\left(x_n+\frac{h}{\text{2}},y_n+\frac{h{k}_1}{2}\right)$ ,

$k_3=f\left(x_n+\frac{h}{\text{2}},y_n+\frac{h{k}_2}{2}\right)$ ,

$k_4=f\left(x_n+h,y_n+h{k}_3\right)$ ,

$y_{n+1}=y_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right).$ (4)

3.1.2. 应用于Blasius方程

为了将RK4方法应用于Blasius方程，我们首先将(1)转化为一组一阶微分方程。设定：

$\begin{array}{l}{y}_1=f,\\ y_2=f^{\prime },\\ y_3=f^{″}.\end{array}$ (5)

则Blasius方程可重写为：

$\begin{array}{l}\frac{\text{d}{y}_1}{\text{d}\eta }=y_2,\\ \frac{\text{d}{y}_2}{\text{d}\eta }=y_3,\\ \frac{\text{d}{y}_3}{\text{d}\eta }=-\frac{1}{2}{y}_1{y}_3,\end{array}$ (6)

在求解过程中，我们设定初始条件$y_1\left(0\right)=0$ ，$y_2\left(0\right)=0$ ，而$y_3\left(0\right)=f^{″}\left(0\right)$ 为未定参数。为了满足${y^{\prime }}_2\left(\infty \right)=1$ ，Raghu Karthik Sadasivuni等人通过采用打靶法进行动态调整。通过迭代的方式调整初始条件$f^{″}\left(0\right)$ ，使最终解逐步接近目标值。具体步骤如下：

• 设定初始条件$f^{″}\left(0\right)$ 的初始猜测(例如1.0)。

• 通过RK4方法计算$f^{\prime }\left(\eta \right)$ 的值。

• 计算$f^{\prime }\left(\infty \right)$ (实际上是$f^{\prime }\left({\eta }_{end}\right)$ 的近似值)。

• 根据计算结果与目标值之间的差异，利用二分法动态调整初始条件，重复此过程，直到结果满足预设精度。

为方便起见，我们将这种使用打靶法动态调整的Runge-Kutta简记为SRK。

3.2. Newton法动态调整的四阶Runge-Kutta法

在求解Blasius方程的过程中，满足边界条件$f^{\prime }\left(\infty \right)=1$ 是一个关键挑战。为了克服这一难题，我们引入了Newton法来动态调整初始条件$f^{″}\left(0\right)$ ，以确保最终解能够满足边界条件。

3.2.1. 方法步骤

Newton法[10]，也称为牛顿–拉弗森方法(Newton-Raphson method)，是一种迭代求解方程根的数值方法。该方法通过线性化非线性方程的根来寻找方程的近似解。对于给定的非线性方程$f\left(x\right)=0$ ，Newton法的迭代公式为：

$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}$ , (7)

其中，$f\left(x_n\right)$ 是目标函数在$x_n$ 处的值，$f^{\prime }\left(x_n\right)$ 是目标函数在$x_n$ 处的导数。

3.2.2. 目标函数与导数的计算

在本文中，目标函数为：

$g\left(f^{″}\left(0\right)\right)=f^{\prime }\left(\infty \right)-1$ . (8)

我们的目标是通过调整初始条件$f^{″}\left(0\right)$ 使得$g\left(f^{″}\left(0\right)\right)=0$ ，即$f^{\prime }\left(\infty \right)=1$ 。

为了实施Newton迭代，我们需要计算目标函数的导数$g^{\prime }\left(f^{″}\left(0\right)\right)$ 。在实际应用中，直接计算导数可能较为复杂，因此我们采用中心差分法来近似导数：

$g^{\prime }\left(f^{″}\left(0\right)\right)\approx \frac{g^{\prime }\left(f^{″}\left(0\right)+\Delta \right)-g^{\prime }\left(f^{″}\left(0\right)-\Delta \right)}{2\Delta }$ , (9)

其中$\Delta$ 为一个小的增量，以确保差分近似的准确性。本文中，我们取$\Delta ={10}^{-6}$ 。

3.2.3. 迭代更新与收敛条件

根据Newton法的更新公式，我们迭代更新初始条件：

$f^{″}{\left(0\right)}_{new}=f^{″}{\left(0\right)}_{old}-\frac{g\left(f^{″}{\left(0\right)}_{old}\right)}{g^{\prime }\left(f^{″}{\left(0\right)}_{old}\right)}$ . (10)

迭代过程持续进行，直到满足预设的收敛条件。收敛条件通常设置为：

$\left|g\left(f^{″}\left(0\right)\right)\right|<\epsilon$ ,

其中$\epsilon$ 为预设的容忍度，用于控制迭代的精度。

通过将四阶Runge-Kutta法与Newton法结合，我们可以在运行过程中有效调整初始条件，迅速满足边界条件，从而提高求解的准确性和效率。初始阶段利用Runge-Kutta法获取近似解，在后续步骤中通过Newton法进行动态微调，以求得更精确的结果。为方便，我们将这种使用Newton法结合的四阶Runge-Kutta法简记为NRK。

4. 数值实验

在本节中，我们将通过使用改进的Runge-Kutta法结合Newton迭代法求解Blasius方程边值问题，并与传统打靶法及SRK方法进行对比。

4.1. 数值实验设置

我们设定的Blasius方程如下：

$f^{‴}+\frac{1}{2}f{f}^{″}=0$ ,

定义边界条件为：

$f\left(0\right)=0,f^{\prime }\left(0\right)=0,f^{\prime }\left(\infty \right)=1$ .

在数值实验中，设定结束的${\eta }_{end}=10$ ，步长h选择为 0.1，初始猜测$f^{″}\left(0\right)=1$ 。预设的容忍度$\epsilon ={10}^{-9}$ ，最大迭代次数为$ite{r}_{\text{max}}=1000$ 。我们分别使用四阶Runge-Kutta法结合Newton法和传统的打靶法进行求解。

4.2. 运行结果

本节将展示采用改进的Runge-Kutta法结合Newton迭代法求解Blasius方程边值问题的实验结果。我们将对比分析新方法与传统打靶法及结合打靶法的Runge-Kutta法(SRK)的性能差异。具体来说，我们将关注迭代次数、计算时间以及在边界条件$f^{″}\left(0\right)$ 处的计算精度。这些指标将帮助我们评估新方法在收敛速度和求解精度上的优势。

表1为三种格式在同一条件的情况下的迭代次数、计算时间和$f^{″}\left(0\right)$ 的计算结果。从表中可以看出，打靶法在收敛至预期值时需要大量的迭代，而SRK方法仅需要20次迭代即可达到高精度的解，相比而言，本文的NRK法迭代次数更少，显著减少了计算时间。而从最终的$f^{″}\left(0\right)$ 的结果可以看出，两种Runge-Kutta都取得了与经典结论0.33206一致的结果，而打靶法却没有达到。

Table 1. Comparison of iteration times, computation time, and calculation results of $f^{″}\left(0\right)$

表1. 迭代次数、计算时间及$f^{″}\left(0\right)$ 的计算结果对比

方法	迭代次数	计算时间(s)	$f^{″}\left(0\right)$
打靶法	27,001	51.9142	0.3308222675879946562
SRK	20	0.0359	0.3320573462609329729
NRK	2	0.0309	0.3320573459459646437

考虑到打靶法需要两万多步才能达到要求，我们并没有将结果呈现在文章之中。表2为两种Runge-Kutta法的迭代误差，对比可以清晰地看到两种改进的Runge-Kutta方法在收敛速度和精度上的差异：

Table 2. Comparison of iteration step errors

表2. 迭代步误差对比

迭代步	SRK法	NRK法
1	1.08540902e+00	1.08540902e+00
2	2.77606469e−01	2.22044605e−16
3	8.50898571e−02
4	2.75948094e−02
5	9.11493542e−03
6	3.02912691e−03
7	1.00869116e−03
8	3.36117394e−04
9	1.12026579e−04
10	3.73407980e−05
11	1.24467776e−05
12	4.14890857e−06
13	1.38296759e−06
14	4.60988977e−07
15	1.53662966e−07
16	5.12209852e−08
17	1.70736609e−08
18	5.69122016e−09
19	1.89707317e−09
20	6.32358166e−10

NRK法的误差在第二次迭代就已经达到了足够的要求，而SRK则需要二十步，且NRK的误差要远小于SRK。说明Runge-Kutta法与Newton迭代法的结合在求解Blasius方程时提供了显著的收敛速度优势。

5. 结论

本文成功地开发并验证了一种结合四阶Runge-Kutta法和Newton迭代法的新型数值求解方法，用于求解流体力学中的关键问题——Blasius方程的边值问题。通过将Blasius方程转化为一阶微分方程组，并应用改进的Runge-Kutta法，结合Newton迭代法动态调整初始条件，本研究显著提高了求解效率和精度。

实验结果表明，与传统的打靶法和结合打靶法的Runge-Kutta法(SRK)相比，新方法在迭代次数和计算时间上均有显著减少，同时在求解精度上也更为准确。具体来说，新方法(NRK)在仅两次迭代后即可达到与经典解析解0.33206一致的结果，而打靶法需要两万多次迭代，SRK方法也需要二十次迭代。此外，NRK法的迭代误差在第二次迭代后已降至足够小，远低于SRK法。

本研究的结论强调了结合Runge-Kutta法和Newton迭代法在处理边值问题时的优越性。该方法不仅提高了收敛速度，还确保了求解的高精度，这对于工程应用和科学计算中的流体力学问题具有重要的实用价值。未来的工作可以进一步探索该方法在其他类型的边界层方程和更广泛的流体力学问题中的应用，以及进一步优化算法以适应更复杂的流动条件。