边界条件为控制函数的双曲方程的数值方法
A Numerical Method for Hyperbolic Equations with Boundary Conditions of Control Functions
DOI: 10.12677/aam.2025.143109, PDF, HTML, XML,    科研立项经费支持
作者: 侯 仟, 潘桂林, 杜绍洪*:重庆交通大学数学与统计学院,重庆
关键词: 双曲方程控制函数解析解数值解收敛速度Hyperbolic Equation Control Function Analytical Solution Numerical Solution Convergence Rate
摘要: 在这篇文章里,我们考虑源自于控制消振问题的双曲方程,由于它的边值条件是需求解的控制函数,导致它不同于通常的偏微分方程的定解问题。在没有关于初始位移和初始速度的任何假设下,给出了控制函数的解析形式,截取级数的前2(N + 1)项作为数值逼近解,获得了数值逼近的收敛速度,数值实验验证了理论结果。
Abstract: In this paper, we are concerned with the development of numerical method for a class of hyperbolic equations, which are originated from the problem of vibration control and elimination in practice. Since its boundary condition is the control function which is needed to be solved, it is different from the ordinary initial value and boundary value problems of partial differential equations. Without any assumptions on initial displacement and initial velocity, the analytical form of the control function is given, and a numerical approximation solution is obtained by intercepting the first 2(N + 1) terms of the series, and a convergence rate of numerical approximation is analyzed. Numerical experiments confirm the theoretical results.
文章引用:侯仟, 潘桂林, 杜绍洪. 边界条件为控制函数的双曲方程的数值方法[J]. 应用数学进展, 2025, 14(3): 229-236. https://doi.org/10.12677/aam.2025.143109

1. 引言

导弹、宇宙飞船在高速飞行的过程中容易产生剧烈的振动,这些振动将影响飞行精度、破坏飞行器的稳定性等。最终导致飞行物偏离预定的运行轨道。例如,导弹在发射过程中,在发射速度与强烈气流的共同作用下会引起导弹的剧烈振动,这些振动不仅会损害导弹结构的稳定性,同时会使导弹无法按照预定的轨道运行,最终导致无法精准击中目标。因此,控制消振问题就变得极为重要。为了解决飞行器的振动问题,考虑将飞行器视为杆件,在两端施加相同的控制来消除纵振动的物理模型可归结为以下的双曲方程的定解问题[1]

{ u tt a 2 u xx =00<x<l, t>0 u| x=0 = u| x=l =v( t ) u| t=0 =φ( x ) u t | t=0 =ψ( x ) (1)

其中, a 是正常数, v( t ) 是控制函数, u( t,x ) 是飞行器的位移。由于控制是在飞行器起飞后进行的,因此,可假设 v( 0 )= v ( 0 )=0 。控制目标是通过在两端施加相同的控制力 v( t ) ,在物理实现中需要将其转化为控制输入,如推进系统的推力和方向。在规定的时间内,将飞行器的纵向振动衰减至零,即选择适当的控制函数 v( t ) ,使得振动 u 在指定的时间 T 之后消除,即 u( T,x )= u t ( T,x )=0

对于绝大多数偏微分方程(组)的定解问题,由于无法求出解析解,求近似解的各种数值方法便应运而生。如今,数值方法已成为大规模科学计算和求解各类工程实践问题的强有力工具,它涉及有限差分方法[2]、有限元方法[3] (包括自适应有限元方法[4])、谱方法[5],以及基于深度神经网络的数值逼近方法[6],这些数值方法包括有限差分方法、有限元方法以及谱方法的研究和应用已持续了数十年的时间,而基于深度神经网络的逼近方法的研究才刚刚起步。

然而,这些成熟的数值方法和基于物理信息的神经网络都不能用来求解问题(1),由于双曲方程(1)的边值条件是需求解的控制函数,导致方程(1)不能用数值方法来离散。因此,广泛使用的这些数值方法并不适宜于求解问题(1),而经典的分离变量方法[7]将成为一个好的选择。

采用分离变量方法来求(1)的形式解,通常需假设 φ( x ) (初始位移)和 ψ( x ) (初始速度)关于 x=l/2 对称。在这篇文章里,我们去掉了这个假设,给出了它的解析解,然后截取级数的前 2( N+1 ) 项作为它的数值解,分析了逼近解的收敛速度。

本文组织如下:在第二节,我们证明了级数截断后的收敛速度;在第三节,我们给出了解析解;在第四节,我们给出了控制函数的一个数值逼近收敛阶;在第五节,我们进行了数值实验证明了理论分析的正确性。

2. 预备引理

在这篇文章里,用 ( , ) 表示 [ 0,l ] 上的 L 2 内积, H k ( [ 0,l ] ) 是标准的索伯列夫空间, H k | | H k 分别是 H k ( [ 0,l ] ) 上的范数和半范数。当 k=0 时,索伯列夫空间 H 0 ( [ 0,l ] )= L 2 ( [ 0,l ] ) ,范数 H 0 = L 2 ;当 k= 时,范数 H =

引理1 设 f( x ) [ 0,l ] 上的周期函数,并有傅里叶展开

f( x )= a 0 2 + n=1 a n cos 2nπx l + b n sin 2nπx l (2)

其中, a n b n f( x ) 的傅里叶系数。 f N ( x ) f( x ) 的截取前 2N+1 项后的逼近,即

f N ( x )= a 0 2 + n=1 N a n cos 2nπx l + b n sin 2nπx l (3)

如果 f( x ) H k ( [ 0,l ] ) ( k>1 ) ,那么成立如下误差估计

f( x ) f N ( x ) l k c( k ) 2 k1 π k N k1 | f | H k (4)

f( x ) f N ( x ) l k1/2 c( k ) 2 k3/2 π k N k1 | f | H k (5)

其中, c( k )= 1 1 k + 1 2 k + 1 3 k +

证明:使用分部积分,并且注意到表达式有 f (n) ( 0 )= f (n) ( l ), 0nk ,得到

a n = 2 l 0 l f ( x )cos 2nπx l dx ={ ( 1 ) (k+1)/2 l k1 2 k1 ( nπ ) k 0 l sin 2nπx l f (k) ( x )dx k ( 1 ) k/2 l k1 2 k1 ( nπ ) k 0 l cos 2nπx l f (k) ( x )dx  k (6)

据(6),有

| a n | l k1/2 2 k1/2 1 ( nπ ) k | f | H k (7)

同样得到

| b n |=| 2 l 0 l f( x )sin 2nπx l dx | l k1/2 2 k1/2 1 ( nπ ) k | f | H k (8)

从(7)和(8),得到

f( x ) f N ( x ) = n=N+1 a n cos 2nπx l + b n sin 2nπx l l k | f | H k 2 k1 π k n=N+1 1 n k l k c( k ) 2 k1 π k 1 N k1 | f | H k (9)

类似地得到(5)。

3. 解析解的推导

定理1 设 v( t ) 是(1)的解,对于 0<tl/a ,它具有如下解析形式

v( t )= 1 2 n=0 ( φ 2n+1 sin ( 2n+1 )πat l + l ( 2n+1 )πa ψ 2n+1 cos ( 2n+1 )πat l ) (10)

其中, φ 2n+1 ψ 2n+1 分别是函数 φ( x ) ψ( x ) 按照 { sin ( 2n+1 )πx l } 展开的傅里叶系数。

证明:令 u( t,x )=w( t,x )+v( t ) ,则 w( t,x ) 满足以下定解问题

{ w tt a 2 w xx = v ( t )0<x<l, t>0 w| x=0 =0   w| x=l =0 w| t=0 =φ( x ) w t | t=0 =ψ( x ) (11)

于是

w( t,x )= k=1 sin kπ l x [ φ k cos kπ l at + l kπa 0 t 2( cos( kπ )1 ) kπ v ( τ )sin kπa( tτ ) l dτ + l kπa ψ k sin kπ l at ] (12)

这里

φ k  = 2 l 0 l φ( x )sin kπx l dx (13)

ψ k  = 2 l 0 l ψ( x )sin kπx l dx (14)

利用分部积分公式,我们可以得到

w( t,x )= k=1 sin kπ l x [ φ k cos kπ l at 2a( cos( kπ )1 ) l 0 t v( τ )sin kπa( tτ ) l dτ + l kπa ψ k sin kπ l at ]v( t ) (15)

于是

u( t,x )= k=1 sin kπ l x [ φ k cos kπ l at 2a( cos( kπ )1 ) l 0 t v( τ )sin kπa( tτ ) l dτ + l kπa ψ k sin kπ l at ] (16)

因为 u( T,x )= u t ( T,x )=0 ,所以

φ k cos kπ l aT+ l kπa ψ k sin kπ l aT= 2a( cos( kπ )1 ) l 0 T v( τ )sin kπa( Tτ ) l dτ (17)

kπa l φ k sin kπ l aT+ ψ k cos kπ l aT= 2a( cos( kπ )1 ) l kπa l 0 T v( τ )cos kπa( Tτ ) l dτ (18)

φ k cos kπ l aT+ l kπa ψ k sin kπ l aT= 2a( cos( kπ )1 ) l 0 T v( τ )sin kπa( Tτ ) l dτ (19)

φ k sin kπ l aT+ l kπa ψ k cos kπ l aT= 2a( cos( kπ )1 ) l 0 T v( τ )cos kπa( Tτ ) l dτ (20)

等式(19)与(20)表明右端可视为 v( t ) 关于三角函数系

{ sin kπa( Tt ) l cos kπa( Tt ) l (21)

的傅里叶展开,而左边是已知的,从而便可以写出未知函数 v( t ) 的表达式。

k=2n 时,

φ k cos kπ l aT+ l kπa ψ k sin kπ l aT=0 (22)

φ k sin kπ l aT+ l kπa ψ k cos kπ l aT=0 (23)

那么 φ 2n = ψ 2n =0 ,即偶数项不出现。

k=2n+1 时,

φ 2n+1 cos ( 2n+1 )π l aT+ l ( 2n+1 )πa ψ 2n+1 sin ( 2n+1 )π l aT= 4a l 0 T v( τ )sin ( 2n+1 )πa( Tτ ) l dτ (24)

φ 2n+1 sin ( 2n+1 )π l aT+ l ( 2n+1 )πa ψ 2n+1 cos ( 2n+1 )π l aT= 4a l 0 T v( τ )cos ( 2n+1 )πa( Tτ ) l dτ (25)

我们取 T= l a ,令 t= l a τ ,有

φ 2n+1 = 4a l 0 l a v( l a t )sin ( 2n+1 )π l a tdt (26)

l ( 2n+1 )πa ψ 2n+1 = 4a l 0 l a v( l a t )cos ( 2n+1 )π l a tdt (27)

然后,令 v ˜ ( t )=v( l a t ) ,有

1 2 φ 2n+1 = 2 l a 0 l a v ˜ ( t )sin ( 2n+1 )π l a tdt (28)

1 2 l ( 2n+1 )πa ψ 2n+1 = 2 l a 0 l a v ˜ ( t )cos ( 2n+1 )π l a tdt (29)

v ˜ ( t ) 关于三角函数系

{ sin ( 2n+1 )π l a t cos ( 2n+1 )π l a t (30)

的傅里叶展开为

v ˜ ( t )= n=0 1 2 φ 2n+1 sin ( 2n+1 )π l a t 1 2 l ( 2n+1 )πa ψ 2n+1 cos ( 2n+1 )π l a t (31)

于是,当 0<tT=l/a 时,可得

v( t )= n=0 ( 1 2 φ 2n+1 sin ( 2n+1 )π l a ( l a t ) 1 2 l ( 2n+1 )πa ψ 2n+1 cos ( 2n+1 )π l a ( l a t ) ) (32)

化简(32),便得到(10)。

4. 数值逼近

上节推导的控制函数 v( t ) 是一个无穷级数,在实践中,自然截取前 2( N+1 ) 项作为它的数值逼近,即

v N ( t )= 1 2 n=0 N ( φ 2n+1 sin ( 2n+1 )πat l + l ( 2n+1 )πa ψ 2n+1 cos ( 2n+1 )πat l ) (33)

其中,N的选取是任意的,但随着N的增大,计算量也会随之增加,因此,实践中只需选取适当的N即可。如当控制较为复杂时,需要尽可能多的级数项作为精确解,此时N的选取尽可能大;当控制较为简单时,此时较低的N就足够了。由于 φ 2n+1 ψ 2n+1 通常不能精确计算,实践中常使用高精度的复合高斯求积公式,而不影响截断误差的阶。

由于振动在指定的时间之后消除,这意味着在T时刻之后无需再加控制,从而可理解为函数值 v( T )= v ( T )=0 。据延拓定理,可以将原本在有限的时间内有定义且在之后为零的控制函数延拓为一个周期函数。因此,可假设 v( t ) [ 0,T ] 上的周期函数,便可得截断误差估计。

定理2 设 v( t ) v N ( t ) 分别满足(10)和(33)。假设边界控制函数 v( t ) [ 0,T ] 上的周期函数,并且 v( t ) H k ( [ 0,T ] ) ,那么成立如下估计

v( t ) v N ( t ) =O( N (k1) ) (34)

v( t ) v N ( t ) =O( N (k1) ) (35)

证明:使用引理1,便可得这两个估计。

5. 数值实验

本节我们用两个数值算例来测试数值算法的收敛性。根据(10),控制函数 v( t ) 是一个无穷级数,其精确解是未知的。因此,根据(33),我们截取前1000项作为精确解来测试数值解按不同范数的收敛速度。为此,把区间 [ 0,l/a ] 分成 l/ ( aτ ) 等分,这里 τ 为网格步长。用 t j =jτ 表示第 j 个节点坐标,定义离散极大范误差 v( t ) v N ( t ) ,τ 和离散 L 2 范误差 v( t ) v N ( t ) 0,τ 如下:

v( t ) v N ( t ) ,τ := max 1jl/(aτ) | v( t j ) v N ( t j ) | (36)

v( t ) v N ( t ) 0,τ := ( j=1 l/ ( aτ ) ( v( t j ) v N ( t j ) ) 2 τ ) 1 2 (37)

5.1. 算例1

在(1)中考虑 a=0.25 l=π/2 ,初始位移 φ( x )=sinx+cosx1 ,初始速度 ψ( x )=0 。这个选择给了 x[ 0,π/2 ] t( 0,2π ] 。在计算中,网格步长 τ=0.001π ,傅里叶系数是精确计算的。

表1是数值误差分别按离散极大范和离散L2范的计算结果,从中可以看到:数值误差随截断项数N的增大而减少,并且按离散极大范具有2阶收敛速度,而按离散L2范具有2.4阶的收敛速度,这个收敛结果遵循了理论估计。

Table 1. Numerical error and convergence rate of example 1

1. 算例1的数值误差和收敛速度

N

v( t ) v N ( t ) ,τ

收敛阶

v( t ) v N ( t ) 0,τ

收敛阶

10

6.304× 10 5

1.245× 10 2

15

2.981× 10 5

1.847

4.891× 10 3

2.304

20

1.731× 10 5

1.889

2.481× 10 3

2.360

25

1.129× 10 5

1.915

1.455× 10 3

2.391

30

7.947× 10 6

1.926

9.376× 10 4

2.410

5.2. 算例2

在(1)中考虑 a=1 l=2π ,并且初始位移 φ( x )=cosx1 ,初始速度 ψ( x )=cosx1 。这个选择给了 x[ 0,2π ] t( 0,2π ] 。在计算中,网格步长仍取 τ=0.001π ,傅里叶系数仍是精确计算的。表2是数值误差分别按离散极大范和离散L2范的计算结果,从中观察到了与算例1一样的结果,再一次保持了与理论分析的一致性,也进一步表明了数值方法的有效性。

Table 2. Numerical error and convergence rate of example 2

2. 算例2的数值误差和收敛速度

N

v( t ) v N ( t ) ,τ

收敛阶

v( t ) v N ( t ) 0,τ

收敛阶

10

1.039× 10 3

2.009× 10 1

15

4.864× 10 4

1.873

7.857× 10 2

2.315

20

2.810× 10 4

1.907

3.978× 10 2

2.366

25

1.825× 10 4

1.934

2.332× 10 2

2.394

30

1.283× 10 4

1.933

1.502× 10 2

2.413

6. 结论及展望

本文研究了一种边界条件为控制函数的双曲方程数值解法,在没有关于初始速度和初始位移的任何假设下,利用分离变量法给出了控制函数的解析形式,截取前 2( N+1 ) 项作为近似解,得到数值逼近的收敛速度,证明数值方法是有效的。然而,仍有部分问题值得进一步研究,本文中的研究目标是纵向振动在一定的时间后衰减为零,因此可用延拓定理假设控制函数为周期函数,而在实际的工程应用中只需纵向振动衰减到一定范围内即可,此时控制函数不一定为周期函数,数值方法还可以进一步改进。

基金项目

重庆市教委科学技术研究计划重大项目(KJZD-M202300705);重庆交通大学研究生科研创新项目(2023S0107)。

NOTES

*通讯作者。

参考文献

[1] 刘盾. 实用数学物理方程[M]. 重庆: 重庆大学出版社, 1999.
[2] 李荣华, 冯果忱. 微分方程数值解法[M]. 北京: 高等教育出版社, 2001.
[3] Susanne, C., Brenner, L. and Ridgway, S. (1994) The Mathematical Theory of Finite Element Methods. Springer-Verlag.
[4] Du, S.H., Lin, R.C. and Zhang, Z.M. (2021) Robust Recovery-Type A Posteriori Error Estimators for Streamline Upwind/Petrov Galerkin Discretizations for Singularly Perturbed Problems. Applied Numerical Mathematics, 168, 23-40.
https://doi.org/10.1016/j.apnum.2021.05.020
[5] Shen, J. and Tang, T. (2006) Spectral and High-Order Methods with Applications. Science Press.
[6] He, J.C., et al. (2020) ReLU Deep Neural Networks and Linear Finite Elements. Journal of Computational Mathematics, 38, 502-527.
https://doi.org/10.4208/jcm.1901-m2018-0160
[7] 王元明. 数学物理方程与特殊函数[M]. 北京: 高等教育出版社, 2009.