两点边值问题的GalerKin方法
1.问题的提出
考虑一个定义在
[
a
,
b
]
[a,b]
[a,b]上的二阶常微分方程边值问题:
−
u
¨
(
t
)
+
q
(
t
)
u
(
t
)
=
f
(
t
)
,
a
<
t
<
b
u
(
a
)
=
0
,
u
(
b
)
=
0
-\ddot u(t)+q(t)u(t) = f(t),a<t<b\\ u(a) = 0,u(b) = 0
−u¨(t)+q(t)u(t)=f(t),a<t<bu(a)=0,u(b)=0
其中
q
,
f
∈
C
[
a
,
b
]
q,f \in C[a,b]
q,f∈C[a,b]为已知的函数,
∀
t
∈
[
a
,
b
]
,
q
(
t
)
≥
0
\forall t\in[a,b],q(t) \geq0
∀t∈[a,b],q(t)≥0。
2.求解方法(试函数法)
我们定义以下
n
n
n个线性无关的基函数:
ϕ
1
(
t
)
,
ϕ
2
(
t
)
,
.
.
.
ϕ
n
(
t
)
ϕ
i
(
a
)
=
ϕ
i
(
b
)
=
0
\phi_1(t),\phi_2(t),...\phi_n(t)\\\phi_i(a) = \phi_i(b)=0
ϕ1(t),ϕ2(t),...ϕn(t)ϕi(a)=ϕi(b)=0
我们将上述问题中的
u
(
t
)
u(t)
u(t)表述为上述基函数的线性组合。
u
(
t
)
∼
U
(
t
)
=
∑
j
=
1
n
x
j
ϕ
j
(
t
)
u(t)\sim U(t) = \sum_{j = 1}^nx_j\phi_j(t)
u(t)∼U(t)=j=1∑nxjϕj(t)
现在的关键问题就是要先设定一个这样的基函数
ϕ
i
(
t
)
\phi_i(t)
ϕi(t):
例如一般可以选择以下两个函数
ϕ i ( t ) = s i n ( i π t − a b − a ) , i = 1 , 2... n \phi_i(t) = sin(i\pi\frac{t-a}{b-a}),i = 1,2...n ϕi(t)=sin(iπb−at−a),i=1,2...n
ϕ i ( t ) = ( t − a ) ( b − t ) t i , i = 1 , 2... n \phi_i(t) = (t-a)(b-t)t^i,i = 1,2...n ϕi(t)=(t−a)(b−t)ti,i=1,2...n
我们将
U
(
t
)
U(t)
U(t)带入上述的二阶常微分方程,求出残余量:
R
(
t
;
x
)
=
−
∑
j
=
1
n
x
j
ϕ
¨
j
(
t
)
+
∑
j
=
1
n
x
j
q
(
t
)
ϕ
j
(
t
)
−
f
(
t
)
R(t;x) = -\sum_{j = 1}^nx_j\ddot \phi_j(t)+\sum_{j= 1}^nx_jq(t)\phi_j(t)-f(t)
R(t;x)=−j=1∑nxjϕ¨j(t)+j=1∑nxjq(t)ϕj(t)−f(t)
那么如果我们可以找到满足条件
x
=
(
x
1
,
x
2
,
.
.
x
n
)
T
x = (x_1,x_2,..x_n)^T
x=(x1,x2,..xn)T使得
R
(
t
;
x
)
=
‾
0
R(t;x)\overline = 0
R(t;x)=0。但是这实际上是不可能找到这样的函数的,那么我们就只能退而求其次,用最小二乘法求
x
∗
x^{*}
x∗:
F
(
x
∗
)
=
m
i
n
x
∈
R
n
∫
a
b
R
(
t
;
x
)
2
d
x
F(x^*) = min_{x\in R^n}\int_a^bR(t;x)^2dx
F(x∗)=minx∈Rn∫abR(t;x)2dx
那么此时我们用最小二乘法的方程来求此
x
∗
x^*
x∗:
∂
F
∂
x
j
=
0
j
=
1
,
2..
n
\frac{\partial F}{\partial x_j} = 0\quad j = 1,2..n
∂xj∂F=0j=1,2..n
因此可以求出来:
∂
F
∂
x
j
=
2
∫
a
b
R
(
t
;
x
)
(
−
ϕ
¨
j
(
t
)
+
q
(
t
)
ϕ
j
(
t
)
)
d
t
=
0
\frac{\partial F}{\partial x_j} = 2\int_a^bR(t;x)(-\ddot \phi_j(t)+q(t)\phi_j(t))dt = 0
∂xj∂F=2∫abR(t;x)(−ϕ¨j(t)+q(t)ϕj(t))dt=0
化简一下:
2
∫
a
b
(
−
ϕ
¨
j
(
t
)
+
q
(
t
)
ϕ
j
(
t
)
)
(
∑
i
=
1
n
(
−
ϕ
¨
i
(
t
)
+
q
(
t
)
ϕ
i
(
t
)
)
x
i
)
d
t
−
2
∫
a
b
f
(
t
)
(
−
ϕ
¨
j
(
t
)
+
q
(
t
)
ϕ
j
(
t
)
)
d
t
=
0
2\int_a^b(-\ddot \phi_j(t)+q(t)\phi_j(t))(\sum_{i = 1}^n(-\ddot\phi_i(t)+q(t)\phi_i(t))x_i)dt-2\int_a^bf(t)(-\ddot \phi_j(t)+q(t)\phi_j(t))dt =0
2∫ab(−ϕ¨j(t)+q(t)ϕj(t))(i=1∑n(−ϕ¨i(t)+q(t)ϕi(t))xi)dt−2∫abf(t)(−ϕ¨j(t)+q(t)ϕj(t))dt=0
因此上面的式子可以化成矩阵形式:
A
x
⃗
=
b
⃗
A\vec x = \vec b
Ax=b
其中
b
j
=
∫
a
b
f
(
t
)
(
−
ϕ
¨
j
(
t
)
+
q
(
t
)
ϕ
j
(
t
)
)
d
t
b_j = \int_a^bf(t)(-\ddot \phi_j(t)+q(t)\phi_j(t))dt
bj=∫abf(t)(−ϕ¨j(t)+q(t)ϕj(t))dt,
a
i
j
=
∫
a
b
(
−
ϕ
¨
i
+
q
ϕ
i
)
(
−
ϕ
¨
j
+
q
ϕ
j
)
d
t
a_{ij} = \int_a^b(-\ddot \phi_i+q\phi_i)(-\ddot\phi_j+q\phi_j)dt
aij=∫ab(−ϕ¨i+qϕi)(−ϕ¨j+qϕj)dt。求解上面的线性方程组便可以求出最合适的
x
∗
x^*
x∗。
如果我们采用Galerkin方法,我们就是要求满足以下条件的
x
∗
x^*
x∗:
∫
a
b
R
(
t
;
x
)
ϕ
i
(
t
)
d
t
=
0
i
=
1
,
2...
n
\int_a^bR(t;x)\phi_i(t)dt = 0\\ i = 1,2...n
∫abR(t;x)ϕi(t)dt=0i=1,2...n
同理也可以求出以下线性方程组:
B
y
⃗
=
c
⃗
B\vec y = \vec c
By=c
其中
B
i
j
=
∫
a
b
(
ϕ
˙
j
ϕ
˙
i
+
q
ϕ
j
ϕ
i
)
d
t
B_{ij} = \int_a^b(\dot \phi_j \dot \phi_i+q\phi_j\phi_i)dt
Bij=∫ab(ϕ˙jϕ˙i+qϕjϕi)dt,
c
j
=
∫
a
b
f
ϕ
j
d
t
,
c_j = \int_a^bf \phi_jdt,
cj=∫abfϕjdt,
j
=
1
,
2..
n
j = 1,2..n
j=1,2..n。
为了提高精度,我们可以不断的增加
n
n
n,也就是基函数的个数,由逼近定理:
u
(
t
)
=
l
i
m
n
→
∞
U
n
(
t
)
(
U
n
(
t
)
两
端
为
0
,
二
次
可
微
分
)
u(t ) = lim_{n \rightarrow \infty}U_n(t)(Un(t)两端为0,二次可微分)
u(t)=limn→∞Un(t)(Un(t)两端为0,二次可微分)
但是缺点也是明显的,就是说整体多项式的方法对于给定区域很难满足事先给定的边界条件。
3.方法的改进
那么我们可以尝试下面的改进方法
我们将
[
a
,
b
]
[a,b]
[a,b]划分成一些充分小的区间:
a
=
t
0
<
t
1
<
t
2
.
.
.
<
t
n
+
1
=
b
a = t_0<t_1<t_2...<t_{n+1} = b
a=t0<t1<t2...<tn+1=b
我们记有以下式子的成立:
h
i
=
t
i
+
1
−
t
i
h
=
m
a
x
0
≤
i
≤
n
h
i
h_i = t_{i+1} - t_i\\ h = max_{0\leq i \leq n}h_i
hi=ti+1−tih=max0≤i≤nhi
这时候我们的基函数
ϕ
i
(
t
)
\phi_i(t)
ϕi(t)可以记录为满足以下条件的
k
k
k此多项式:
(1)在 [ t i , t i + 1 ] [t_i,t_{i+1}] [ti,ti+1]上 ϕ i ( t ) \phi_i(t) ϕi(t)为 k k k此多项式。
(2)在整个区间 [ a , b ] [a,b] [a,b]上连续。
(3)满足边界条件 ϕ i ( a ) = ϕ i ( b ) = 0 \phi_i(a) = \phi_i(b) = 0 ϕi(a)=ϕi(b)=0。
我们取
k
=
1
k = 1
k=1,即是
ϕ
i
(
t
)
=
a
i
t
+
b
i
(
t
∈
[
t
i
,
t
i
+
1
]
)
\phi_i(t) = a_it+b_i(t\in[t_i,t_{i+1}])
ϕi(t)=ait+bi(t∈[ti,ti+1])
ϕ
i
(
t
)
=
{
t
−
t
i
−
1
h
i
−
1
,
t
∈
[
t
i
−
1
,
t
i
]
t
i
+
1
−
t
h
i
,
t
∈
[
t
i
,
t
i
+
1
]
0
,
e
l
s
e
\phi_i(t) = \begin{cases}\frac{t-t_{i-1}}{h_{i-1}},t\in[t_{i-1},t_i]\\ \frac{t_{i+1}-t}{h_i},t \in[t_{i},t_{i+1}]\\ 0,else\end{cases}
ϕi(t)=⎩⎪⎨⎪⎧hi−1t−ti−1,t∈[ti−1,ti]hiti+1−t,t∈[ti,ti+1]0,else
其导数可以描述为:
ϕ
˙
i
(
t
)
=
{
1
h
i
−
1
,
t
∈
[
t
i
−
1
,
t
i
]
−
1
h
i
,
t
∈
[
t
i
,
t
i
+
1
]
0
,
e
l
s
e
\dot \phi_i(t) = \begin{cases}\frac{1}{h_{i-1}},t\in[t_{i-1},t_i]\\ -\frac{1}{h_i},t \in[t_{i},t_{i+1}]\\ 0,else\end{cases}
ϕ˙i(t)=⎩⎪⎨⎪⎧hi−11,t∈[ti−1,ti]−hi1,t∈[ti,ti+1]0,else
而
ϕ
i
(
t
i
)
\phi_i(t_{i})
ϕi(ti)只有在
(
i
−
1
,
i
+
1
)
(i-1,i+1)
(i−1,i+1)上才不为0,所以有:
ϕ
i
(
t
)
ϕ
j
(
t
)
=
0
ϕ
˙
i
(
t
)
ϕ
˙
j
(
t
)
=
0
\phi_i(t)\phi_j(t) = 0\\ \dot \phi_i(t)\dot \phi_j(t) = 0
ϕi(t)ϕj(t)=0ϕ˙i(t)ϕ˙j(t)=0
因此有:
B
i
i
=
∫
a
b
(
(
ϕ
˙
i
(
t
)
)
2
+
q
(
ϕ
i
(
t
)
)
2
)
d
t
B
i
i
−
1
=
∫
a
b
(
ϕ
˙
i
(
t
)
ϕ
˙
i
−
1
(
t
)
+
q
ϕ
i
(
t
)
ϕ
i
−
1
(
t
)
)
d
t
B
i
i
+
1
=
∫
a
b
(
ϕ
˙
i
(
t
)
ϕ
˙
i
+
1
(
t
)
+
q
ϕ
i
(
t
)
ϕ
i
+
1
(
t
)
)
d
t
B_{ii} = \int_a^b((\dot \phi_i(t))^2+q(\phi_i(t))^2)dt\\ B_{ii-1} = \int_a^b(\dot \phi_i(t) \dot \phi_{i-1}(t)+q\phi_i(t)\phi_{i-1}(t))dt\\ B_{ii+1} = \int_a^b(\dot \phi_i(t) \dot \phi_{i+1}(t)+q\phi_i(t)\phi_{i+1}(t))dt\\
Bii=∫ab((ϕ˙i(t))2+q(ϕi(t))2)dtBii−1=∫ab(ϕ˙i(t)ϕ˙i−1(t)+qϕi(t)ϕi−1(t))dtBii+1=∫ab(ϕ˙i(t)ϕ˙i+1(t)+qϕi(t)ϕi+1(t))dt
我们对上面的式子化简可以得到:
Q
1
,
i
=
1
h
i
2
∫
t
i
t
i
+
1
(
t
i
+
1
−
t
)
(
t
−
t
i
)
q
(
t
)
d
t
,
i
=
1
,
2...
n
−
1
Q
2
,
i
=
1
h
i
−
1
2
∫
t
i
−
1
t
i
(
t
−
t
i
−
1
)
2
q
(
t
)
d
t
,
i
=
1
,
2...
n
Q
3
,
i
=
1
h
i
2
∫
t
i
t
i
+
1
(
t
−
t
i
+
1
)
2
q
(
t
)
d
t
,
i
=
1
,
2...
n
Q
4
,
i
=
1
h
i
−
1
∫
t
i
−
1
t
i
(
t
−
t
i
−
1
)
f
(
t
)
d
t
,
i
=
1
,
2...
n
Q
5
,
i
=
1
h
i
∫
t
i
t
i
+
1
(
t
i
+
1
−
t
)
f
(
t
)
d
t
,
i
=
1
,
2...
n
Q_{1,i} = \frac{1}{h_i^2}\int_{t_i}^{t_{i+1}}(t_{i+1}-t)(t - t_{i})q(t)dt,i = 1,2...n-1\\ Q_{2,i} = \frac{1}{h_{i-1}^2}\int_{t_{i-1}}^{t_{i}}(t - t_{i-1})^2q(t)dt,i = 1,2...n\\ Q_{3,i} = \frac{1}{h_{i}^2}\int_{t_{i}}^{t_{i+1}}(t - t_{i+1})^2q(t)dt,i = 1,2...n\\ Q_{4,i} = \frac{1}{h_{i-1}}\int_{t_{i-1}}^{t_{i}}(t - t_{i-1})f(t)dt,i = 1,2...n\\ Q_{5,i} = \frac{1}{h_{i}}\int_{t_{i}}^{t_{i+1}}(t_{i+1} - t)f(t)dt,i = 1,2...n\\
Q1,i=hi21∫titi+1(ti+1−t)(t−ti)q(t)dt,i=1,2...n−1Q2,i=hi−121∫ti−1ti(t−ti−1)2q(t)dt,i=1,2...nQ3,i=hi21∫titi+1(t−ti+1)2q(t)dt,i=1,2...nQ4,i=hi−11∫ti−1ti(t−ti−1)f(t)dt,i=1,2...nQ5,i=hi1∫titi+1(ti+1−t)f(t)dt,i=1,2...n
矩阵系数是:
B
i
i
=
Q
2
,
i
+
Q
3
,
i
+
1
h
i
−
1
+
1
h
i
B
i
,
i
+
1
=
Q
1
,
i
−
1
h
i
,
i
=
1
,
2..
n
−
1
B
i
,
i
−
1
=
Q
1
,
i
−
1
−
1
h
i
−
1
,
i
=
2
,
3..
n
c
i
=
Q
4
,
i
+
Q
5
,
i
B_{ii} = Q_{2,i}+Q_{3,i}+\frac{1}{h_{i-1}}+\frac{1}{h_i}\\ B_{i,i+1} = Q_{1,i}-\frac{1}{h_i},i = 1,2..n-1\\ B_{i,i-1} = Q_{1,i-1}-\frac{1}{h_{i-1}},i =2,3..n\\ c_i = Q_{4,i}+Q_{5,i}
Bii=Q2,i+Q3,i+hi−11+hi1Bi,i+1=Q1,i−hi1,i=1,2..n−1Bi,i−1=Q1,i−1−hi−11,i=2,3..nci=Q4,i+Q5,i
4.代码实现
例如用有限元法求解以下的方程组:
−
u
¨
+
π
2
u
=
2
π
2
s
i
n
π
t
0
<
t
<
1
u
(
0
)
=
0
u
(
1
)
=
0
-\ddot u + \pi^2u = 2\pi^2sin \pi t \quad 0<t<1\\ u(0) = 0\quad u(1) = 0
−u¨+π2u=2π2sinπt0<t<1u(0)=0u(1)=0
在这里取
h
i
=
h
=
0.1
,
t
i
=
0.1
i
,
i
=
0
,
1..9
h_i = h = 0.1,t_i = 0.1i,i = 0,1..9
hi=h=0.1,ti=0.1i,i=0,1..9。
clc,clear;
q = @(t)pi^2;
f = @(t)2*pi^2*sin(pi*t);
h = 0.1;
Q = zeros(5,10);
B = zeros(10,10);
c = zeros(10,1);
t = @(i)h*i;
for i = 1:10
if i< 10
Q(1,i) = 1/h^2*quad(@(x)(t(i+1)-x).*(x-t(i))*pi^2,t(i),t(i+1));
else
Q(1,i) = 1/h^2*quad(@(x)(t(10)-x).*(x-t(9))*pi^2,t(9),t(10));
end
Q(2,i) = 1/h^2*quad(@(x)(x - t(i-1)).^2*pi^2,t(i-1),t(i));
Q(3,i) = 1/h^2*quad(@(x)(x - t(i+1)).^2*pi^2,t(i),t(i+1));
Q(4,i) = 1/h*quad(@(x)(x - t(i-1))*2*pi^2.*sin(pi*x),t(i-1),t(i));
Q(5,i) = 1/h*quad(@(x)(t(i+1) - x)*2*pi^2.*sin(pi*x),t(i),t(i+1));
B(i,i) = Q(2,i)+Q(3,i)+1/h+1/h;
if i<=9
B(i,i+1) = Q(1,i) - 1/h;
end
if i>1
B(i,i-1) = Q(1,i-1)-1/h;
end
c(i) = Q(4,i)+Q(5,i);
end
x = inv(B)*c
结果是:
x =
0.3165
0.6032
0.8335
0.9864
1.0489
1.0177
0.8992
0.7100
0.4750
0.2261