求解第一类奇异初值问题的新方法

《Journal of Computational and Applied Mathematics》：New Approach to Solving Singular Initial Value Problems of the First Kind

　　

摘要
提出了一种解决一阶非线性常微分方程组初值问题的新方法，该方程组在初始点 x=0 处具有第一类奇点。利用自变量替换 x=e^t，我们将区间 [0, a] 上的原始初值问题转化为区间 (?∞,ln a] 上的等价问题。为了在有限网格 {tn∈(?∞,ln a], n=0,1,…,N, tN=ln a} 的网格节点 t0 处找到近似解，构造了新的一阶方法。对于其他网格节点处的解，可以使用经典的数值方法。对于所提出的方法，我们开发了一种新的有效算法，用于在给定容差下数值求解奇点 x=0 附近的奇异初值问题。

引言
具有第一类奇点的常微分方程（ODE）初值问题（IVP）在数学建模领域及其方法开发中引起了研究人员的关注（见例如 [1], [2], [3], [4], [5], [6], [7], [8], [9]）。这类问题在物理学和天体物理学中具有重要相关性。它们被用于建模雪崩积聚 [5], [10]、恒星结构、球形气体云的热行为、等温气体球和热离子电流 [2], [3], [4], [11]。使用奇异边值问题（BVPs）可以描述更大一类现象 [6], [12], [13], [14]。在这方面，对奇异 IVPs 的特别兴趣也与以下事实有关：应用打靶法或高阶三点差分格式（见 [[13], [14]）来求解奇异 BVPs 需要求解相关的奇异 IVPs。这也产生了构建高效奇异 IVP 方法和算法的需求。

奇点的存在不允许正确应用标准方法。在 [15], [16] 中已经确定，将 Runge-Kutta 和线性多步法应用于奇异 IVPs 会导致这些方法的阶数降低。这一事实不允许正确使用任何已知的用于实际误差估计和步长自动选择的算法。

在本文中，我们考虑以下奇异 IVP
dudx=M(x)u(x)x+f(x,u(x)),x∈(0,a],f(x,u):[0,a]×Rs→Rs,u:[0,a]→Rs,M(x)∈Rs×s,u(0)=A,
其中 M(0)u(0)=0, A∈ker(M(0))。设矩阵 M(x) ∈ C^(1)[0, a]。这意味着该矩阵可以写成以下形式
M(x)=M(0)+xB(x)
其中 B(x) ∈ C[0, a]。还假设 f(x,u)∈C([0,a]×Rs) 且 f(x, u) 在 [0,a]×Rs 上关于 u 满足 Lipschitz 条件。在 [17] 中表明，在这种限制下，初值问题 (1), (2) 在 [0,a]×Rs 上具有唯一连续可微解。

应该注意，对于求解奇异 IVPs，可以使用配点法（见例如 [12], [18], [19]）。通常，配点法应用于奇异边值问题（BVPs），后者比 IVPs 更难以数值求解。众所周知（见 [20, p.212], 定理 7.7），对于 IVPs，配点法等价于隐式 Runge-Kutta 方法。然而，上述方法是隐式的，需要巨大的计算成本通过牛顿迭代法求解非线性代数方程组，特别是在大维度 ODE 系统的情况下。

在 [21] 中对于奇异 IVP
1xλddx[xλK(x)dudx]=?f(x,u),x∈(0,a],f(x,u):(0,a]×Rm→Rm,u:[0,a]→Rm,K(x)∈Rm×m,λ>0,u(0)=A,du(0)dx=0,
我们提出了泰勒级数方法和显式 Runge-Kutta 类型方法，它们不会丢失阶数。注意，该问题可以转化为双维度（s=2m）的问题 (1), (2)，具有矩阵 M(x)∈R2m×2m
M(x)=(0xK?1(x)0?λI),
其中 I∈Rm×m 是单位矩阵，0∈Rm×m 是零矩阵。

在本文中，我们将 [21], [22] 中提出的方法扩展到更一般的问题 (1), (2)，该问题将问题 (1*), (2*) 作为特例。通过替换 x=e^t，我们将问题 (1), (2) 简化为无限区间 (?∞,ln a] 上的 IVP。为了找到该 IVP 在奇点附近，即某些有限不规则网格 ω^h={tn∈(?∞,ln a], n=0,1,…,N, tN=ln a} 的节点 t0 处的数值解，我们构造了泰勒级数方法和四阶 Runge-Kutta 类型方法。为了找到网格其他节点处的数值解，我们应用显式四阶 Runge-Kutta 方法。这种方法保证了整个区间上的四阶精度（见误差估计与收敛性，定理 1）。此外，在计算奇点附近的解之后，我们不仅可以使用的显式 Runge-Kutta 方法，还可以使用线性多步法（见 [21]）。还开发了一种用于在给定容差下数值求解奇点附近 IVPs 的误差估计和步长选择的有效算法，该算法基于嵌入式 Runge-Kutta 类型方法。与隐式 Runge-Kutta 方法不同，所提出的方法在数值积分的每一步都不需要额外的迭代过程。我们将通过我们的方法获得的计算结果与基于隐式 Runge-Kutta 配点类型方法的结果（即与代码 RADAU5（见 [23]），它是这些方法的最佳实现之一）进行了比较。数值例子展示了所提出方法的显著优势。

部分片段
奇点附近的单步法
问题 (1)–(3) 可以写为
dudx=M(0)u(x)x+f^(x,u(x)),x∈(0,a],u(0)=A,
其中 f^(x,u(x))=f(x,u(x))+B(x)u(x)。利用替换 x=e^t，我们将问题 (4), (5) 简化为以下形式
dU(t)dt=M(0)U(t)+etf?(t,U(t)),limt→?∞U(t)=A,
其中 U(t)=u(et),f?(t,U(t))=f^(x,u(x))。
类似于 [22]，在半无限区间 (?∞,b], b≤ln a 上，我们选择有限不规则网格 ω^h={tn∈(?∞,b], n=0,1,…,N, tN=b}，步长 hn=tn?tn?1>0, n=1,2,…,N。
我们在...处假设以下双重不等式

误差估计与收敛性
初值问题 (6), (7) 可以重写为
dU(t)dt=F(t,U),?∞其中 f(t,u)=m(0)u(t)+etf?(t,u(t))。
根据所提出的方法，问题 (26) 的数值解可以通过逐步过程获得
y0=a+(i?m(0))?1h0g1+(i?12m(0))?1h0∑i=16b1igi+(i?13m(0))?1h0∑i=16b2igi+(i?14m(0))?1h0∑i=16b3igi,
g1=f^(0,a),
g2=f^(c2h0,a+h0a21g1),
g3=f^(c3h0,a+h0(a31g1+a32g2)),
g4=f^(c4h0,a+h0(a41g1+a42g2+a43g3)),
g5=f^(c5h0,a+h0(a51g1+a52g2+a53g3+a54g4)),
g6=f^(c6h0,a+

奇异 ivp 数值求解算法
直接实现具有表 1 给定系数的 runge-kutta 类型方法 (23) 需要计算值
g1=f^(0,a),
g2=f^(13h0,a+13h0g1),
g3=f^(25h0,a+h0(425g1+625g2)),
g4=f^(h0,a+h0(14g1?3g2+154g3)),
g5=f^(23h0,a+h0(227g1+109g2?5081g3+881g4)),
g6=f^(45h0,a+h0(350g1+35g2+115g3+215g4?350g5)),
求解四个线性代数方程组
(i?m(0))g1=g1,
(2i?m(0))g2=2(?4g1+42524g3?223g4?3518g5+752g6),
(3i?m(0))g3=3(7912g1?302572g3+1769g4+9278g5?100g6),
(4i?m(0))g4=4(?5516g1+

数值示例
例 1
为了阐明理论，我们考虑以下线性奇异 ivp（见 [1]）
dudx=1x(01?2?8x+x2?3)(u1u2)+(012xe1?x),x∈(0,1],u1(0)=0,u2(0)=0.
该问题的精确解是
u1(x)=x2e1?x,u2(x)=(2?x)x2e1?x.
我们以给定容差求解问题 (33)，在网格 ω^h 的节点 t0（奇点附近）使用我们的算法。在网格的其他节点，我们应用了 dormand-prince 5(4) 的嵌入公式（代码 dopri5 [20]）。求解结果见表 2，其中 nstep 和 nfun。

b.y. datsko|m.v. kutniv
热舒夫理工大学，8 powstancow warszawy, rzeszow, 35-959, 波兰 a,limt→?∞u(t)=A, 其中 f(t,u)=M(0)U(t)+etf?(t,U(t))。 根据所提出的方法，问题 (26) 的数值解可以通过逐步过程获得 y0=A+(I?M(0))?1h0g1+(I?12M(0))?1h0∑i=16b1igi+(I?13M(0))?1h0∑i=16b2igi+(I?14M(0))?1h0∑i=16b3igi, g1=f^(0,A), g2=f^(c2h0,A+h0a21g1), g3=f^(c3h0,A+h0(a31g1+a32g2)), g4=f^(c4h0,A+h0(a41g1+a42g2+a43g3)), g5=f^(c5h0,A+h0(a51g1+a52g2+a53g3+a54g4)), g6=f^(c6h0,A+ 奇异 ivp 数值求解算法 直接实现具有表 1 给定系数的 runge-kutta 类型方法 (23) 需要计算值 g1=f^(0,A), g2=f^(13h0,A+13h0g1), g3=f^(25h0,A+h0(425g1+625g2)), g4=f^(h0,A+h0(14g1?3g2+154g3)), g5=f^(23h0,A+h0(227g1+109g2?5081g3+881g4)), g6=f^(45h0,A+h0(350g1+35g2+115g3+215g4?350g5)), 求解四个线性代数方程组 (i?m(0))g1=g1, (2i?m(0))g2=2(?4g1+42524g3?223g4?3518g5+752g6), (3i?m(0))g3=3(7912g1?302572g3+1769g4+9278g5?100g6), (4i?m(0))g4=4(?5516g1+ 数值示例 例 1 为了阐明理论，我们考虑以下线性奇异 ivp（见 [1]） dudx=1x(01?2?8x+x2?3)(u1u2)+(012xe1?x),x∈(0,1],u1(0)=0,u2(0)=0. 该问题的精确解是 u1(x)=x2e1?x,u2(x)=(2?x)x2e1?x. 我们以给定容差求解问题 (33)，在网格 ω^h 的节点 t0（奇点附近）使用我们的算法。在网格的其他节点，我们应用了 dormand-prince 5(4) 的嵌入公式（代码 dopri5 [20]）。求解结果见表 2，其中 nstep 和 nfun。 b.y. datsko|m.v. kutniv 热舒夫理工大学，8 powstancow warszawy, rzeszow, 35-959,>其中 f(t,u)=m(0)u(t)+etf?(t,u(t))。
根据所提出的方法，问题 (26) 的数值解可以通过逐步过程获得
y0=a+(i?m(0))?1h0g1+(i?12m(0))?1h0∑i=16b1igi+(i?13m(0))?1h0∑i=16b2igi+(i?14m(0))?1h0∑i=16b3igi,
g1=f^(0,a),
g2=f^(c2h0,a+h0a21g1),
g3=f^(c3h0,a+h0(a31g1+a32g2)),
g4=f^(c4h0,a+h0(a41g1+a42g2+a43g3)),
g5=f^(c5h0,a+h0(a51g1+a52g2+a53g3+a54g4)),
g6=f^(c6h0,a+

奇异 ivp 数值求解算法
直接实现具有表 1 给定系数的 runge-kutta 类型方法 (23) 需要计算值
g1=f^(0,a),
g2=f^(13h0,a+13h0g1),
g3=f^(25h0,a+h0(425g1+625g2)),
g4=f^(h0,a+h0(14g1?3g2+154g3)),
g5=f^(23h0,a+h0(227g1+109g2?5081g3+881g4)),
g6=f^(45h0,a+h0(350g1+35g2+115g3+215g4?350g5)),
求解四个线性代数方程组
(i?m(0))g1=g1,
(2i?m(0))g2=2(?4g1+42524g3?223g4?3518g5+752g6),
(3i?m(0))g3=3(7912g1?302572g3+1769g4+9278g5?100g6),
(4i?m(0))g4=4(?5516g1+

数值示例
例 1
为了阐明理论，我们考虑以下线性奇异 ivp（见 [1]）
dudx=1x(01?2?8x+x2?3)(u1u2)+(012xe1?x),x∈(0,1],u1(0)=0,u2(0)=0.
该问题的精确解是
u1(x)=x2e1?x,u2(x)=(2?x)x2e1?x.
我们以给定容差求解问题 (33)，在网格 ω^h 的节点 t0（奇点附近）使用我们的算法。在网格的其他节点，我们应用了 dormand-prince 5(4) 的嵌入公式（代码 dopri5 [20]）。求解结果见表 2，其中 nstep 和 nfun。

b.y. datsko|m.v. kutniv
热舒夫理工大学，8 powstancow warszawy, rzeszow, 35-959, 波兰>