二阶常微分方程的数值解法（中心差分法和有限体积法）

这里我们介绍中心差分法和有限体积法求解方程。
题目：
用差分法的中心差分格式和有限体积法求解两点边值问题
u′′−α(2x−1)u′−2αu=0,0<x<1u(0)=u(1)=1,u^{\prime\prime}-\alpha\left(2x-1\right)u^\prime-2\alpha u=0,0<x<1 u\left(0\right)=u\left(1\right)=1,u′′−α(2x−1)u′−2αu=0,0<x<1u(0)=u(1)=1,
其中，参数α=10\alpha=10α=10，得到不同网格最大误差和收敛阶。

问题分析

Step 1：网格剖分：首先取 N + 1 个节点为： a=x0<x1<x2<⋯<xN=ba=x_0<x_1<x_2<\cdots<x_N=ba=x0<x1<x2<⋯<xN=b,
将区间 [a,b] 作等距剖分，划分为 N 个小区间,记 h=xi+1−xi,i=1,2,⋯,N−1h=x_{i+1}-x_i,i=1,2,\cdots,N-1h=xi+1−xi,i=1,2,⋯,N−1为网格步长.再进行对偶剖分：取相邻节点 xi+1,xix_{i+1} ,x_ixi+1,xi 的中点xi+12=12(xi+1−xi),i=1,2,⋯,N−1.x_{i+\frac{1}{2}}=\frac{1}{2}\left(x_{i+1}-x_i\right),i=1,2,\cdots,N-1.xi+21=21(xi+1−xi),i=1,2,⋯,N−1.由这些节点构成的剖分为对偶剖分.
(i)中心差分法：
1h2(ui+1−2ui+ui−1)−α(2xi−1)2h(ui+1−ui−1)−2αui=0\frac{1}{h^2}\left(u_{i+1}-2u_i+u_{i-1}\right)-\frac{\alpha\left(2x_i-1\right)}{2h}\left(u_{i+1}-u_{i-1}\right)-2\alpha u_i=0h21(ui+1−2ui+ui−1)−2hα(2xi−1)(ui+1−ui−1)−2αui=0
⇒\Rightarrow⇒
(1h2−ri)ui+1+(−2h2−2α)ui+(1h2+ri)ui−1=0\left(\frac{1}{h^2}-r_i\right)u_{i+1}+\left(-\frac{2}{h^2}-2\alpha\right)u_i+\left(\frac{1}{h^2}+r_i\right)u_{i-1}=0(h21−ri)ui+1+(−h22−2α)ui+(h21+ri)ui−1=0,

where ri=α(2xi−1)2h.r_i=\frac{\alpha\left(2x_i-1\right)}{2h}.ri=2hα(2xi−1).
Let A=[−2h2−2α1h2−r10⋯001h2+r2−2h2−2α1h2−r2⋯00⋯⋯⋯⋯⋯⋯000⋯1h2+rN−1−2h2−2α],A=\left[\begin{matrix}-\frac{2}{h^2}-2\alpha&\frac{1}{h^2}-r_1&0&\cdots&0&0\\\frac{1}{h^2}+r_2&-\frac{2}{h^2}-2\alpha&\frac{1}{h^2}-r_2&\cdots&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&\frac{1}{h^2}+r_{N-1}&-\frac{2}{h^2}-2\alpha\\\end{matrix}\right],A=−h22−2αh21+r2⋯0h21−r1−h22−2α⋯00h21−r2⋯0⋯⋯⋯⋯00⋯h21+rN−100⋯−h22−2α,
U=(u1,u2,⋯,uN−1)T,b=(−u0(1h2+r1),0,0,⋯,0,−uN(1h2+rN−1))TU=\left(u_1,u_2,\cdots,u_{N-1}\right)^T,\\ b=\left(-u_0\left(\frac{1}{h^2}+r_1\right),0,0,\cdots,0,-u_N\left(\frac{1}{h^2}+r_{N-1}\right)\right)^TU=(u1,u2,⋯,uN−1)T,b=(−u0(h21+r1),0,0,⋯,0,−uN(h21+rN−1))T.
Then AU=b.
接着解出U即可.

(ii)有限体积法：
Let W(x)=dudxW\left(x\right)=\frac{du}{dx}W(x)=dxdu, then dWdx−α(2x−1)W(x)−2αu=0.\frac{dW}{dx}-\alpha\left(2x-1\right)W\left(x\right)-2\alpha u=0.dxdW−α(2x−1)W(x)−2αu=0.
两边同时积分得：
W(xi+12)−W(xi−12)−∫xi−12xi+12α(2x−1)W(x)dx−∫xi−12xi+122αudx=0.W(x_{i+\frac{1}{2}})-W(x_{i-\frac{1}{2}})-\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}α(2x-1)W(x)dx-\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}2αudx=0.W(xi+21)−W(xi−21)−∫xi−21xi+21α(2x−1)W(x)dx−∫xi−21xi+212αudx=0.
运用积分中值定理可得：
W(xi+12)−W(xi−12)−W(xi)∫xi−12xi+12α(2x−1)dx−ui∫xi−12xi+122αdx=0.W(x_{i+\frac{1}{2}})-W(x_{i-\frac{1}{2}})-W(x_i)\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}α(2x-1)dx-u_i\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}2αdx=0.W(xi+21)−W(xi−21)−W(xi)∫xi−21xi+21α(2x−1)dx−ui∫xi−21xi+212αdx=0.
⇒\Rightarrow⇒
W(xi+12)−W(xi−12)−W(xi)(xi+12+xi−12−1)hα−2αuih=0W(x_{i+\frac{1}{2}})-W(x_{i-\frac{1}{2}})-W(x_i)(x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}}-1)hα-2α u_ih=0W(xi+21)−W(xi−21)−W(xi)(xi+21+xi−21−1)hα−2αuih=0
⇒\Rightarrow⇒
1h(ui+1−2ui+ui−1)−α4(xi+1+2xi+xi−1)(ui+1−ui−1)+α2(ui+1−ui−1)−2αuih=0\frac{1}{h}\left(u_{i+1}-2u_i+u_{i-1}\right)-\frac{\alpha}{4}\left(x_{i+1}+2x_i+x_{i-1}\right)\left(u_{i+1}-u_{i-1}\right)+\frac{\alpha}{2}\left(u_{i+1}-u_{i-1}\right)-{2\alpha\ u}_ih=0h1(ui+1−2ui+ui−1)−4α(xi+1+2xi+xi−1)(ui+1−ui−1)+2α(ui+1−ui−1)−2α uih=0
⇒\Rightarrow⇒
(1h−li+α2)ui+1+(−2h−2αh)ui+(1h+li−α2)ui−1=0,\left(\frac{1}{h}-l_i+\frac{\alpha}{2}\right)u_{i+1}+\left(-\frac{2}{h}-2\alpha h\right)u_i+\left(\frac{1}{h}+l_i-\frac{\alpha}{2}\right)u_{i-1}=0,(h1−li+2α)ui+1+(−h2−2αh)ui+(h1+li−2α)ui−1=0,
whereli=α4(xi+1+2xi+xi−1){\ l}_i=\frac{\alpha}{4}\left(x_{i+1}+2x_i+x_{i-1}\right) li=4α(xi+1+2xi+xi−1).
Let A=[−2h−2αh1h−l1+α20⋯001h+l2−α2−2h−2αh1h−l2+α2⋯00⋯⋯⋯⋯⋯⋯000⋯1h+lN−1−α2−2h−2αh],A=\left[\begin{matrix}-\frac{2}{h}-2\alpha h&\frac{1}{h}-l_1+\frac{\alpha}{2}&0&\cdots&0&0\\\frac{1}{h}+l_2-\frac{\alpha}{2}&-\frac{2}{h}-2\alpha h&\frac{1}{h}-l_2+\frac{\alpha}{2}&\cdots&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&\frac{1}{h}+l_{N-1}-\frac{\alpha}{2}&-\frac{2}{h}-2\alpha h\\\end{matrix}\right],A=−h2−2αhh1+l2−2α⋯0h1−l1+2α−h2−2αh⋯00h1−l2+2α⋯0⋯⋯⋯⋯00⋯h1+lN−1−2α00⋯−h2−2αh,

U=(u1,u1,⋯,uN−1)T,F=(−u0(1h+l1−α2),0,0,⋯,0,−uN(1h−lN−1+α2))TU=\left(u_1,u_1,\cdots,u_{N-1}\right)^T,F=\left(-u_0\left(\frac{1}{h}+l_1-\frac{\alpha}{2}\right),0,0,\cdots,0,-u_N\left(\frac{1}{h}-l_{N-1}+\frac{\alpha}{2}\right)\right)^TU=(u1,u1,⋯,uN−1)T,F=(−u0(h1+l1−2α),0,0,⋯,0,−uN(h1−lN−1+2α))T
Then AU=F.
接着解出U即可.

全部代码如下(matlab)

function E=Centered(h)%中心差分法
u0=1;
uN=1;
N=1/h;
x=zeros(1,N+1);
for i=1:(N)x(i+1)=x(i)+h;
end
R=zeros(1,N+1);
for i=1:(N+1)R(i)=10*(2*x(i)-1)/(2*h);
end
A=zeros(N-1,N-1);
A(1,1)=-2/(h*h)-20;
A(1,2)=1/(h*h)-R(2);
A(N-1,N-1)=-2/(h*h)-20;
A(N-1,N-2)=1/(h*h)+R(N)
for i=2:(N-2)A(i,i)=-2/(h*h)-20;A(i,i+1)=1/(h*h)-R(i+1);A(i,i-1)=1/(h*h)+R(i+1);
end
b=zeros(1,N-1);
b(1)=-u0*(1/(h*h)+R(2));
b(N-1)=-uN*(1/(h*h)-R(N));
y=A\b';
Y=zeros(N+1,1);
Y(1,1)=1;
Y(N+1,1)=1;
for i=2:NY(i,1)=y(i-1,1);
end
plot(x,Y ,'r*-');
hold on;
exa = dsolve('D2u = 10*(2*x-1)*Du+20*u', 'u(0)=1','u(1)=1', 'x');       %求出解析解
ezplot(exa, [0,1]);       %画出解析解的图像
legend('中心差分法数值解','解析解' );
u=zeros(1,N+1);
exa=inline(exa);
for i=1:(N+1)u(i)=feval(exa,x(i));e(i)=abs(u(i)-Y(i));
end
E=max(e);function E=V(h)%有限体积法
u0=1;
uN=1;
N=1/h;
x=zeros(1,N+1);
for i=1:(N)x(i+1)=x(i)+h;
end
R=zeros(1,N+1);
for i=2:NR(i)=2.5*(x(i+1)+2*x(i)+x(i-1));
end
A=zeros(N-1,N-1);
A(1,1)=-2/h-20*h;
A(1,2)=1/h+5-R(2);
A(N-1,N-1)=-2/h-20*h;
A(N-1,N-2)=1/h-5+R(N);
for i=2:(N-2)A(i,i)=-2/h-20*h;A(i,i+1)=1/h+5-R(i+1);A(i,i-1)=1/h-5+R(i+1);
end
b=zeros(1,N-1);
b(1)=-u0*(1/h-5+R(2));
b(N-1)=-uN*(1/h+5-R(N));
y=A^(-1)*b';
Y=zeros(N+1,1);
Y(1,1)=1;
Y(N+1,1)=1;
for i=2:NY(i,1)=y(i-1,1);
end
plot(x,Y ,'r*-');
hold on;
exa = dsolve('D2u = 10*(2*x-1)*Du+20*u', 'u(0)=1','u(1)=1', 'x');       %求出解析解
ezplot(exa, [0,1]);       %画出解析解的图像
legend('有限体积法数值解','解析解' );
u=zeros(1,N+1);
exa=inline(exa);
for i=1:(N+1)u(i)=feval(exa,x(i));e(i)=abs(u(i)-Y(i));
end
E=max(e);
%调试
clear all;clc;close  all;
h=[0.1,0.2,0.01,0.02,0.05];
for i=1:5E(i)=V(h(i));F(i)=Centered(h(i));
end

问题结果

我们发现中心差分法的最大误差与有限体积法的最大误差相差较小，阶数相同，因此我们可以认为中心差分法与有限体积法两者的收敛阶相同，而我们已知中心差分法的收敛阶为2阶，因此有限体积法的收敛阶也是二阶。

（以h=0.1为例，我们画出了两种方法数值解与解析解的图像）