中心差分格式的数值试验（含MATLAB源码）
考虑二阶常微分方程边值问题：
（1）
其中q，f为[a,b]上的连续函数， 为常数
1、考虑问题：考虑二阶常微分方程边值问题
2、网格剖分与差分格式
3、截断误差
将方程（1）在节点离散化，由泰勒公式展开得截断误差
附件 ：程序源代码
### 中心差分格式的数值试验
#### 一、引言
本文主要介绍了一种针对二阶常微分方程边值问题的数值求解方法——中心差分格式，并通过MATLAB进行了数值试验验证。该方法是数值分析领域常用的一种手段，能够有效逼近微分方程的解。
#### 二、考虑的问题
考虑到一个特定形式的二阶常微分方程边值问题:
\[ -\left(p(x)u'(x)\right)' + q(x)u(x) = f(x), \quad x \in [a, b] \]
边界条件为:
\[ u(a) = \alpha, \quad u(b) = \beta \]
其中 \( p(x), q(x), f(x) \) 分别为定义在区间 \([a, b]\) 上的连续函数，\(\alpha\) 和 \(\beta\) 是常数。我们的目标是找到满足上述方程及边界条件的函数 \(u(x)\)。
#### 三、网格剖分与差分格式
为了数值求解上述问题，首先需要对区间 \([a, b]\) 进行网格剖分。具体做法是将 \([a, b]\) 均匀分割成 \(N\) 段，即每个子区间的长度为 \(h = (b - a)/N\)。这样，可以得到 \(N+1\) 个网格节点 \(x_i = a + ih, i = 0, 1, \ldots, N\)。
在每个网格节点上，采用中心差分近似二阶导数，得到差分格式:
\[ -\frac{p(x_i)}{h^2} \left(u_{i+1} - 2u_i + u_{i-1}\right) + q(x_i)u_i = f(x_i), \quad i = 1, 2, \ldots, N-1 \]
#### 四、截断误差分析
差分格式的精度可以通过截断误差来衡量。根据泰勒公式展开，我们可以得到:
\[ -\left(p(x_i)u'(x_i)\right)' + q(x_i)u(x_i) - f(x_i) = O(h^2) \]
这意味着，当 \(h\) 趋于零时，上述差分格式的解将逐渐接近真实解。
#### 五、数值例子与结果
假设具体的微分方程为:
\[ -\left(p(x)u'(x)\right)' + q(x)u(x) = f(x), \quad x \in [0, 1] \]
其中 \(p(x) = 1, q(x) = 1 + \sin(x), f(x) = e^x \cdot \sin(x)\)，边界条件为 \(u(0) = 1, u(1) = e\)。
根据上述差分格式，可以将向量形式的差分格式转换为矩阵形式:
\[ Au = F \]
其中，
\[ A = \begin{pmatrix} 2 + q_1 h^2 & -1 & 0 & \cdots & 0 \\ -1 & 2 + q_2 h^2 & -1 & \ddots & \vdots \\ 0 & -1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & 2 + q_{N-1} h^2 & -1 \\ 0 & \cdots & 0 & -1 & 2 + q_N h^2 \end{pmatrix}, \quad F = \begin{pmatrix} h^2 f_1 + \alpha \\ h^2 f_2 \\ \vdots \\ h^2 f_{N-1} \\ h^2 f_N + \beta \end{pmatrix} \]
其中 \(\alpha = 1, \beta = e\)。解该线性方程组即可得到 \(u(x)\) 的数值解。
#### 六、结果分析
使用MATLAB编程实现上述过程，并进行数值试验。选取不同的网格密度 \(N\) 进行对比，例如 \(N = 10\) 和 \(N = 100\)。结果表明，随着 \(N\) 的增加，数值解更加接近真实解，且收敛效果明显。
- 当 \(N = 10\) 时，可以看到数值解已经呈现出良好的逼近趋势。
- 当 \(N = 100\) 时，数值解更接近于真实解，曲线也更加平滑。
#### 七、结论
通过本试验可以看出，中心差分格式是一种有效的数值求解二阶常微分方程边值问题的方法。随着网格密度的增加，数值解的精度会显著提高。此外，MATLAB作为工具，不仅方便了编程实现，还使得数值模拟变得更加直观和高效。
#### 八、参考文献
1. 李荣华. 偏微分方程数值解法[M]. 高等教育出版社, 2018.
2. 李维国. 数值计算方法[M]. 中国石油大学出版社, 2015.
### 九、程序源代码
```matlab
%u(x)=exp(x),q(x)=1+sinx,f(x)=exp(x)*sinx0<x<1
a=0; % 输入条件a
b=1; % 输入条件b
N=10; % 输入等分数N
h=(b-a)/N; % 步长
for i=1:N-1
x(i)=a+i*h; % 网格节点
q(i)=1+sin(x(i)); % 连续函数q网格节点上各值
f(i)=exp(x(i))*sin(x(i)); % 连续函数f网格节点上各值
end
% 构造矩阵
c1=linspace(2+q(1)*h*h,2+q(N-1)*h*h,N-1); % 主对角线各值
c2=linspace(-1,-1,N-2); % 次对角线各值
d1=diag(c1,0);
d2=diag(c2,-1);
d3=diag(c2,1);
d=d1+d2+d3
D=inv(d); % 矩阵的逆
b1=linspace(h*h*f(1),h*h*f(N-1),N-1);
for i=1:N-1
if i==1
b2(i)=exp(a);
elseif i==N-1
b2(i)=exp(b);
else
b2(i)=0;
end
end
B1=b1+b2; % 构造向量
B=B1';
u=D*B; % 输出u(x)数值解
plot(a+h:h:b-h,u,'*')
ylabel('y');
xlabel('x');
title('拟合曲线');
```
通过以上步骤，可以清晰地看到中心差分格式应用于二阶常微分方程边值问题的整个过程及其结果。
