16901讲义笔记:加权残差法 这部分内容的主题是有限元方法(FEM),虽然我们是从数学角度来探讨这 个问题,但有很多有趣的物理穿插其中,特别是流体的问题,(例如:以扩散效 应为主的物理现象 一个简单的例子 我们希望解决下面这个问题: 轴向载 荷均匀 *桁架竖直放置,且其轴向载荷均匀 *左右端固定 模型方程 1,且v(0)=0,v(1)=0 进行两次积分,再结合边界条件,我们就可以得到方程的精确解 →v(x)=x2+ax+b v(0)=0=b→b= v(1)=7+a=0→a v(x)=x(x-1 加权残数法是一种逼近控制方程的解的方法,下面我们就来看在这个问题中,它 是怎样进行的。 设方程的解是你选择的函数的线性组合。」
16.901 讲义笔记：加权残差法 这部分内容的主题是有限元方法（FEM），虽然我们是从数学角度来探讨这 个问题，但有很多有趣的物理穿插其中，特别是流体的问题，（例如：以扩散效 应为主的物理现象） 一个简单的例子 我们希望解决下面这个问题： 轴向载 荷均匀 ﹡桁架竖直放置，且其轴向载荷均匀 ﹡左右端固定 模型方程： 2 2 1 d v dx = ，且v v (0) 0, (1) 0 = = 进行两次积分，再结合边界条件，我们就可以得到方程的精确解 1 2 ( ) 2 ⇒ = ++ v x x ax b v b (0) 0 = = ⇒ b = 0 1 (1) 0 2 v a =+= ⇒ 1 2 a = − 1 ( ) ( 1) 2 vx xx = − 加权残数法是一种逼近控制方程的解的方法，下面我们就来看在这个问题中，它 是怎样进行的。 第一步：设方程的解是你选择的函数的线性组合
例如,在这个问题中,我们可以选择正弦函数的形式,也许是 U(r)=a, sin Tx+a, sin 2x+a, sin 3Tx 这里的a1a2a3均为常数,可以通过某种方法确定。sinx,sin2xx和sin3x就 是所选的函数。 *我选择这些函数是因为他们自然的满足v(0)=0,v(1)=0的边界条件 现在,我们定义一个残差,它只是对控制方程进行了很小的整理: R(v) 方程的精确解就是任何时候(例如:0≤x≤1)都满足R(v)=0和边界条件。 让我们用U替换方程中的ν,即: R(U)=2-1a, sin zx+ a2 sin 2mx+a, sin 3Tx R(=-z-asinTx-27 Fa, sin2Tx-(3r a, sin3rx-1 因为我们只有三个未知数(a,a2a3),所以不可能任意地令上式等于零。但有 些办法来解决这个问题: 1,配方法:选择N=3个位置,在这些位置上令R(U=0。如果进展顺利,我们 能够得到合理的a值。 2,最小二乘法:选择a,使得「R取最小 3,加权残差法:选择N=3个加权函数,找出满足[a2R(Ukx=0 「aO)kx=0∫aRO=0的a,这里的a就是权重函数 对于这个问题,我们用上述三种方法来分别讨论。 配方法 配点法的关键就是如何选择x,使得在该点强制满足R(U)=0,(无疑地)这儿会 有一些微妙的问题,例如,如果我们选择 2
例如，在这个问题中，我们可以选择正弦函数的形式，也许是 12 3 Ux a x a x a x ( ) sin sin 2 sin 3 =+ + π π π 这里的 均为常数，可以通过某种方法确定。sin 123 aa a , , π x, sin 2 sin 3 π π x x 和 就 是所选的函数。 ﹡我选择这些函数是因为他们自然的满足v v (0) 0, (1) 0 = = 的边界条件 现在，我们定义一个残差，它只是对控制方程进行了很小的整理： 2 2 () 1 d v R v dx = − 方程的精确解就是任何时候（例如：0 ≤ x ≤1）都满足 R v() 0 = 和边界条件。 让我们用 U 替换方程中的 v，即： { } 2 2 12 3 ( ) sin sin 2 sin 3 1 d R U a xa xa x dx = ++ − π π π () ()2 12 3 2 2 R( ) sin 2 sin2 3 sin3 U ax a x a x =− − π ππ ππ π − −1 因为我们只有三个未知数（ ），所以不可能任意地令上式等于零。但有 些办法来解决这个问题： 123 aaa , , 1， 配方法：选择N=3 个位置，在这些位置上令R(U)=0。如果进展顺利，我们 能够得到合理的ai值。 2， 最小二乘法：选择ai，使得 取最小。 1 2 0 R dx ∫ 3， 加权残差法：选择N=3 个加权函数，找出满足 ， ， 1 1 0 ω R U dx () 0 = ∫ 1 2 0 ω R U dx () 0 = ∫ 1 3 0 ω R U dx () 0 = ∫ 的ai，这里的ωi 就是权重函数 对于这个问题，我们用上述三种方法来分别讨论。 配方法 配点法的关键就是如何选择 x ，使得在该点强制满足 R(U)=0，（无疑地）这儿会 有一些微妙的问题，例如， 如果我们选择： 1 2 3 0 1 2 1 x x x = = =
在x=0处会出问题。让我们来看一下为什么会是这样 在x=0处(注意:x=0时, sin mzx=0) →R(U(0)=-1这表明不能令其为0,因为a的值对方程不起作用。 如果把前面的选择改为 x2 441 x3 我们可以得到下面的方程: 丌 -z'asin-(2r)a, sin*-(3T)a, sin -1=0 在x2处 4n21(27)如mx-(3x)4m2-1=0 在x3处 9丌 丌a.Sin Ta,SI 3兀(3)a,sin 1=0 或者写成矩阵方程,即 4-(2r sin?-(37)" 3丌 T SIn (2z)"sinT -(37)sin T sin 3兀-(2z)sin 2-(3x)m1 求逆(我用 MATLAB)可得 a1=-0.1223,a2=0,a3=-0.0023 注意:a2=0并不令人意外,因为这个问题是关于x=1/2完全对称的,但 sin2丌x不是关于x=1/2对称的。如图
在 x=0 处会出问题。让我们来看一下为什么会是这样。 在 x=0 处（注意：x=0 时，sin mpx=0） ⇒ R U( (0) 1 ) = − ，这表明不能令其为 0，因为ai的值对方程不起作用。 如果把前面的选择改为： 1 2 3 1 4 1 2 3 4 x x x = = = 我们可以得到下面的方程： 在x1处： ( ) ( ) 123 2 2 2 3 sin 2 3 sin sin 424 aaa ππ π −− − π π π −1 0 = 在x2处： 1 23 () () 2 2 2 3 sin 2 3 sin sin 2 2 aaa π π −− − π π ππ −1 0 = 在x3处： ( ) ( ) 12 3 2 2 2 3 9 sin sin sin 3 2 3 424 aa a π π π 1 0 π −− − π π − = 或者写成矩阵方程，即： ( ) ( ) () () ( ) ( ) 2 2 2 1 2 2 2 2 3 2 2 2 3 sin 2 sin 3 sin 424 1 3 sin 2 sin 3 sin 1 2 2 1 339 sin 2 sin 3 sin 424 a a a π π π πππ π π π ππ π πππ ππ π ⎡ ⎤ −− − ⎢ ⎥⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ −− − ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ −− − ⎣ ⎦ 求逆（我用 MATLAB）可得： 1 23 a aa =− = =− 0.1223, 0, 0.0023 注意：a2=0 并不令人意外，因为这个问题是关于x=1/2 完全对称的，但 sin 2π x不是关于x=1/2 对称的。如图
最小二乘法 现在,我们就要试图去寻找a使得「R2最小,从前面的内容我们有: R(U)=-xa, sin zx-(27)a, sin 2zx-(37)a,sin3zX-I 对上式求积分,会有很多下面这种形式的项: 1/2,如果n=m sin mxsin nzxdx= 0.如果n≠m 把等于零的项都去掉,很容易看出: ∫h=amnx+(2)9m2x丁+(ax)5m3x+ 2(x)a, sin zx+2(2m)a, sin 2 rx+2(3x)a, sin 3rx+1 dr (看不见)给出: 0 1a2al0(2x)0 4x.0.12 t a2 00(3n) a' Da+g a 为使之最小化,我们对a求偏导,使得 af 0
最小二乘法 现在，我们就要试图去寻找 使得 i a 1 2 0 R dx ∫ 最小，从前面的内容我们有： 12 3 ( ) ( ) 2 2 2 RU a x a x a x ( ) sin =− − ππ π π 2 3 π π sin 2 − sin 3 −1 对上式求积分，会有很多下面这种形式的项： 1 0 1/ 2, sin sin 0, n m m x n xdx n m π π ⎧ = = ⎨ ⎩ ≠ ∫ 如果 如果 把等于零的项都去掉，很容易看出： { ( ) ( ) 1 1 2 2 2 2 2 2 2 12 3 0 0 R dx a x a x a x =+ + ⎡ ⎤ ππ π π π π sin 2 sin 2 3 sin3 ⎡ ⎤ ⎡ ∫ ∫ ⎣ ⎦ ⎣ ⎦ ⎣ ⎤ + ⎦ ( ) ( ) ( ) } 22 2 12 3 22 2 π ππ ππ π ax a x a x sin 2 sin 2 3 sin3 1 + + + dx （看不见）给出： [ ] ( ) ( ) [ ] 4 1 1 4 123 2 2 4 3 3 0 0 1 , , 0 2 0 4 , 0, 12 1 2 003 a a aa a a a a a π π π π ⎡ ⎤ π ⎡ ⎤ ⎡ ⎢ ⎥ ⎤ ⎢ ⎥ ⎢ = + ⎥ + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎢ ⎥ ⎣ ⎦ ⎦ 1 1 2 T T = + a Da g a G G GG + 为使之最小化，我们对a 求偏导，使得 G 1 0 a ∂f = ∂ 2 0 a ∂f = ∂ 3 0 a ∂f = ∂
其中∫=aDa,注意,常数项1不影响解 这是一个很著名的问题,其解的形式为: 解得:(对于我们这个问题,D是对角阵,所以很容易求解) a1=-0.1290a2=0a3=-00048 加权残差法 我们最后考虑的是加权残差法,在这种方法里,我们令: .oR(Ubx=01=123 O)可以有很多种合理选择,其中最著名的是伽辽金加权残差法,他用同一个函 数去逼近U(x): @,=sinx @,=sin 2Tx @,=sin 3Tx O, R(U)dx=J(sin nx)I-za, sinTx-(2r)a2sin2Tx (37)a, sin 3Tx-l]dx 如果利用前面 Sinmixsinn丌xd 的结果,则有 JooR(XDr=I 「a2R(x==,(2 O,R(UYx=-(37) 这时令它们为零,有 4 4 或 a1=-0.1290a2=0a3=-0.0048
其中 1 2 T f = a Da K K ，注意，常数项 1 不影响解。 这是一个很著名的问题，其解的形式为： Da = − G Gg ( ) g 1 a D − ⇒ =− G G 解得：（对于我们这个问题，D 是对角阵，所以很容易求解） 1 a = −0.1290 2 a = 0 3 a = −0.0048 加权残差法 我们最后考虑的是加权残差法，在这种方法里，我们令： 1 0 () 0 ωi R U dx = ∫ i =1,2,3 ωi 可以有很多种合理选择，其中最著名的是伽辽金加权残差法，他用同一个函 数去逼近 ： U x( ) 1 ω = sinπ x 2 ω = sin 2π x 3 ω = sin 3π x ( ) () 1 1 2 2 1 2 0 0 ( ) sin [ sin 2 sin 2 ωi R U dx x a x a x = ππ π π π − − ∫ ∫ ( )2 3 − 3 sin 3 1] π π a x − dx 如果利用前面 1 的结果，则有： 0 sin sin m x n xdx π π ∫ 1 2 1 1 0 1 2 ( ) 2 R U dx a π π ω =− − ∫ 1 2 2 2 0 1 ( ) (2 ) 2 ω R U dx a = − π ∫ 1 2 3 3 0 3 1 2 ( ) (3 ) 2 R U dx a π π ω =− − ∫ 这时令它们为零，有： 1 3 4 a π = − 2 a = 0 3 2 4 (3 ) a π = − 或 1 a = −0.1290 2 a = 0 3 a = −0.0048
这和最小二乘的结果是完全一样的 Comparison of collocation, least squares, and weighted residuals 0.04 0.08 0.12 least square/mwr(dash-dot)
这和最小二乘的结果是完全一样的