美赛整理之偏微分方程的数值求解(一)

两点边值问题的GalerKin方法

1.问题的提出
2.求解方法(试函数法)
3.方法的改进
4.代码实现

1.问题的提出

考虑一个定义在[a,b][a,b][a,b]上的二阶常微分方程边值问题:
−u¨(t)+q(t)u(t)=f(t),a<t<bu(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.求解方法(试函数法)

我们定义以下nnn个线性无关的基函数：
ϕ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=1nxjϕ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)=sin(iπt−ab−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)ti,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=1nxjϕ¨j(t)+∑j=1nxjq(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=(x1,x2,..xn)Tx = (x_1,x_2,..x_n)^Tx=(x1,x2,..xn)T使得R(t;x)=‾0R(t;x)\overline = 0R(t;x)=0。但是这实际上是不可能找到这样的函数的，那么我们就只能退而求其次，用最小二乘法求x∗x^{*}x∗:
F(x∗)=minx∈Rn∫abR(t;x)2dxF(x^*) = min_{x\in R^n}\int_a^bR(t;x)^2dx F(x∗)=minx∈Rn∫abR(t;x)2dx
那么此时我们用最小二乘法的方程来求此x∗x^*x∗:
∂F∂xj=0j=1,2..n\frac{\partial F}{\partial x_j} = 0\quad j = 1,2..n ∂xj∂F=0j=1,2..n
因此可以求出来:
∂F∂xj=2∫abR(t;x)(−ϕ¨j(t)+q(t)ϕj(t))dt=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∫ab(−ϕ¨j(t)+q(t)ϕj(t))(∑i=1n(−ϕ¨i(t)+q(t)ϕi(t))xi)dt−2∫abf(t)(−ϕ¨j(t)+q(t)ϕj(t))dt=02\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
因此上面的式子可以化成矩阵形式:
Ax⃗=b⃗A\vec x = \vec b Ax=b
其中bj=∫abf(t)(−ϕ¨j(t)+q(t)ϕj(t))dtb_j = \int_a^bf(t)(-\ddot \phi_j(t)+q(t)\phi_j(t))dtbj=∫abf(t)(−ϕ¨j(t)+q(t)ϕj(t))dt，aij=∫ab(−ϕ¨i+qϕi)(−ϕ¨j+qϕj)dta_{ij} = \int_a^b(-\ddot \phi_i+q\phi_i)(-\ddot\phi_j+q\phi_j)dtaij=∫ab(−ϕ¨i+qϕi)(−ϕ¨j+qϕj)dt。求解上面的线性方程组便可以求出最合适的x∗x^*x∗。

如果我们采用Galerkin方法，我们就是要求满足以下条件的x∗x^*x∗：
∫abR(t;x)ϕi(t)dt=0i=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
同理也可以求出以下线性方程组:
By⃗=c⃗B\vec y = \vec c By=c
其中Bij=∫ab(ϕ˙jϕ˙i+qϕjϕi)dtB_{ij} = \int_a^b(\dot \phi_j \dot \phi_i+q\phi_j\phi_i)dtBij=∫ab(ϕ˙jϕ˙i+qϕjϕi)dt,cj=∫abfϕjdt,c_j = \int_a^bf \phi_jdt,cj=∫abfϕjdt,j=1,2..nj = 1,2..nj=1,2..n。

为了提高精度，我们可以不断的增加nnn，也就是基函数的个数，由逼近定理:
u(t)=limn→∞Un(t)(Un(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=t0<t1<t2...<tn+1=ba = t_0<t_1<t_2...<t_{n+1} = b a=t0<t1<t2...<tn+1=b
我们记有以下式子的成立:
hi=ti+1−tih=max0≤i≤nhih_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)可以记录为满足以下条件的kkk此多项式:

(1)在[ti,ti+1][t_i,t_{i+1}][ti,ti+1]上ϕi(t)\phi_i(t)ϕi(t)为kkk此多项式。

(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=1k = 1k=1,即是ϕi(t)=ait+bi(t∈[ti,ti+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−ti−1hi−1,t∈[ti−1,ti]ti+1−thi,t∈[ti,ti+1]0,else\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)={1hi−1,t∈[ti−1,ti]−1hi,t∈[ti,ti+1]0,else\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(ti)\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
因此有：
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))dtB_{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
我们对上面的式子化简可以得到:
Q1,i=1hi2∫titi+1(ti+1−t)(t−ti)q(t)dt,i=1,2...n−1Q2,i=1hi−12∫ti−1ti(t−ti−1)2q(t)dt,i=1,2...nQ3,i=1hi2∫titi+1(t−ti+1)2q(t)dt,i=1,2...nQ4,i=1hi−1∫ti−1ti(t−ti−1)f(t)dt,i=1,2...nQ5,i=1hi∫titi+1(ti+1−t)f(t)dt,i=1,2...nQ_{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
矩阵系数是:
Bii=Q2,i+Q3,i+1hi−1+1hiBi,i+1=Q1,i−1hi,i=1,2..n−1Bi,i−1=Q1,i−1−1hi−1,i=2,3..nci=Q4,i+Q5,iB_{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¨+π2u=2π2sinπt0<t<1u(0)=0u(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
在这里取hi=h=0.1,ti=0.1i,i=0,1..9h_i = h = 0.1,t_i = 0.1i,i = 0,1..9hi=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:10if i< 10Q(1,i) = 1/h^2*quad(@(x)(t(i+1)-x).*(x-t(i))*pi^2,t(i),t(i+1));elseQ(1,i) = 1/h^2*quad(@(x)(t(10)-x).*(x-t(9))*pi^2,t(9),t(10));endQ(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<=9B(i,i+1) = Q(1,i) - 1/h;endif i>1B(i,i-1) = Q(1,i-1)-1/h;endc(i) = Q(4,i)+Q(5,i);
end
x = inv(B)*c

结果是:


x =0.31650.60320.83350.98641.04891.01770.89920.71000.47500.2261