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

1.问题的提出

​ 考虑一个定义在$[a,b]$上的二阶常微分方程边值问题:
$$\begin{gathered} -\ddot{u}(t)+q(t)u(t)=f(t),a<t<b \\ u(a)=0,u(b)=0 \end{gathered}$$
其中$q,f\in C[a,b]$为已知的函数，$\forall t\in [a,b],q(t)\ge 0$。

2.求解方法(试函数法)

我们定义以下$n$个线性无关的基函数：
$$\begin{gathered} {\varphi }_1(t),{\varphi }_2(t),...{\varphi }_n(t) \\ {\varphi }_i(a)={\varphi }_i(b)=0 \end{gathered}$$
我们将上述问题中的$u(t)$表述为上述基函数的线性组合。
$$u(t)\sim U(t)=\sum _{j=1}^n{x}_j{\varphi }_j(t)$$
现在的关键问题就是要先设定一个这样的基函数${\varphi }_i(t)$：

例如一般可以选择以下两个函数

${\varphi }_i(t)=sin(i\pi \frac{t-a}{b-a}),i=1,2...n$

${\varphi }_i(t)=(t-a)(b-t)t^i,i=1,2...n$

我们将$U(t)$带入上述的二阶常微分方程，求出残余量：
$$R(t;x)=-\sum _{j=1}^n{x}_j{\ddot{\varphi }}_j(t)+\sum _{j=1}^n{x}_{j}q(t){\varphi }_j(t)-f(t)$$
那么如果我们可以找到满足条件$x=(x_1,x_2,..x_n{)}^T$使得$R(t;x)\overline{=}0$。但是这实际上是不可能找到这样的函数的，那么我们就只能退而求其次，用最小二乘法求$x^{\ast }$:
$$F(x^{\ast })=mi{n}_{x\in R^n}{\int }_a^{b}R(t;x{)}^{2}dx$$
那么此时我们用最小二乘法的方程来求此$x^{\ast }$:
$$\frac{\partial F}{\partial x_j}=0j=1,2..n$$
因此可以求出来:
$$\frac{\partial F}{\partial x_j}=2{\int }_a^{b}R(t;x)(-{\ddot{\varphi }}_j(t)+q(t){\varphi }_j(t))dt=0$$
化简一下:
$$2{\int }_a^b(-{\ddot{\varphi }}_j(t)+q(t){\varphi }_j(t))(\sum _{i=1}^n(-{\ddot{\varphi }}_i(t)+q(t){\varphi }_i(t))x_i)dt-2{\int }_a^{b}f(t)(-{\ddot{\varphi }}_j(t)+q(t){\varphi }_j(t))dt=0$$
因此上面的式子可以化成矩阵形式:
$$A\vec{x}=\vec{b}$$
其中$b_j={\int }_a^{b}f(t)(-{\ddot{\varphi }}_j(t)+q(t){\varphi }_j(t))dt$，$a_{ij}={\int }_a^b(-{\ddot{\varphi }}_i+q{\varphi }_i)(-{\ddot{\varphi }}_j+q{\varphi }_j)dt$。求解上面的线性方程组便可以求出最合适的$x^{\ast }$。

如果我们采用Galerkin方法，我们就是要求满足以下条件的$x^{\ast }$：
$$\begin{gathered} {\int }_a^{b}R(t;x){\varphi }_i(t)dt=0 \\ i=1,2...n \end{gathered}$$
同理也可以求出以下线性方程组:
$$B\vec{y}=\vec{c}$$
其中$B_{ij}={\int }_a^b({\dot{\varphi }}_j{\dot{\varphi }}_i+q{\varphi }_j{\varphi }_i)dt$,$c_j={\int }_a^{b}f{\varphi }_{j}dt,$$j=1,2..n$。

为了提高精度，我们可以不断的增加$n$，也就是基函数的个数，由逼近定理:
$$u(t)=li{m}_{n\to \infty }{U}_n(t)(Un(t)两端为0，二次可微分)$$
但是缺点也是明显的，就是说整体多项式的方法对于给定区域很难满足事先给定的边界条件。

3.方法的改进

那么我们可以尝试下面的改进方法

我们将$[a,b]$划分成一些充分小的区间:
$$a=t_0<t_1<t_2...<t_{n+1}=b$$
我们记有以下式子的成立:
$$\begin{gathered} h_i=t_{i+1}-t_i \\ h=ma{x}_{0\le i\le n}{h}_i \end{gathered}$$
这时候我们的基函数${\varphi }_i(t)$可以记录为满足以下条件的$k$此多项式:

(1)在$[t_i,t_{i+1}]$上${\varphi }_i(t)$为$k$此多项式。

(2)在整个区间$[a,b]$上连续。

(3)满足边界条件${\varphi }_i(a)={\varphi }_i(b)=0$。

我们取$k=1$,即是${\varphi }_i(t)=a_{i}t+b_i(t\in [t_i,t_{i+1}])$
$${\varphi }_i(t)=\{\begin{array}{l}\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{array}$$

其导数可以描述为:
$${\dot{\varphi }}_i(t)=\{\begin{array}{l}\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{array}$$
而${\varphi }_i(t_i)$只有在$(i-1,i+1)$上才不为0，所以有:
$$\begin{gathered} {\varphi }_i(t){\varphi }_j(t)=0 \\ {\dot{\varphi }}_i(t){\dot{\varphi }}_j(t)=0 \end{gathered}$$
因此有：
$$\begin{gathered} B_{ii}={\int }_a^b(({\dot{\varphi }}_i(t){)}^2+q({\varphi }_i(t){)}^2)dt \\ {B}_{ii-1}={\int }_a^b({\dot{\varphi }}_i(t){\dot{\varphi }}_{i-1}(t)+q{\varphi }_i(t){\varphi }_{i-1}(t))dt \\ {B}_{ii+1}={\int }_a^b({\dot{\varphi }}_i(t){\dot{\varphi }}_{i+1}(t)+q{\varphi }_i(t){\varphi }_{i+1}(t))dt \end{gathered}$$
我们对上面的式子化简可以得到:
$$\begin{gathered} 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}{)}^{2}q(t)dt,i=1,2...n \\ {Q}_{3,i}=\frac{1}{h_i^2}{\int }_{t_i}^{t_{i+1}}(t-t_{i+1}{)}^{2}q(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 \end{gathered}$$
矩阵系数是:
$$\begin{gathered} 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} \end{gathered}$$

4.代码实现

例如用有限元法求解以下的方程组:
$$\begin{gathered} -\ddot{u}+{\pi }^{2}u=2{\pi }^{2}sin\pi t0<t<1 \\ u(0)=0u(1)=0 \end{gathered}$$
在这里取$h_i=h=0.1,t_i=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