目录
目录参考文献
1 李荣华 微分方程的数值解法
零、 偏微分方程
偏微分方程形式
抛物线方程：扩散方程
\[\frac{\partial\mu}{\partial t} =\frac{\partial^2\mu}{\partial x^2}+ \frac{\partial^2\mu}{\partial y^2}+ \frac{\partial^2\mu}{\partial z^2} \]双曲线方程：对流方程
\[\frac{\partial\mu}{\partial t}= \frac{\partial\mu}{\partial x} \]椭圆方程(1):拉普拉斯方程
\[\frac{\partial^2\mu}{\partial x^2}+ \frac{\partial^2\mu}{\partial y^2}+ \frac{\partial^2\mu}{\partial z^2} =0 \]椭圆方程(2):泊松方程
\[\frac{\partial^2\mu}{\partial x^2}+2\frac{\partial^2\mu}{\partial y^2}+3\frac{\partial^2\mu}{\partial z^2}4=f(x,y,z) \]一、常微分方程初值问题
常微分方程形式
\[\frac{du}{dt}=f(t,u), \qquad 0<t<T, \\ u(0)=u_0 \]1 常微分方程求解(1)：欧拉法
泰勒展开
\[u(t_1)=u(t_0+h)=u(t_0)+hu'(t_0)+\frac {h^2}{2}u''(t_0)+\frac{h^3}{6}u'''(\zeta)\\ = u(t_0)+hu'(t_0,u_0)+R_0 \]其中 $ \zeta \in (t_0,t_1) \(,并略去二阶小量\)R_0$
\[u_1=u_0+hf(t_0,t_1) \]然后得到一个递推公式
\[u_{n+1}=u_n+hf(t_n,t_n) \]欧拉法的几何意义
微分方程的解释\(t,u\)平面上的积分曲线族，过一点恰好有一积分曲线通过，按欧拉法，过初始点\((t_0,u_0)\)作经过此点的切线，斜率为\(f(t_0,u_0)\) , 沿切线取点\((t_1,u_1)\), 然后继续做切线。**
1.1 欧拉法与改进欧拉法matlab代码
1）欧拉法
% 无需导入特定库，MATLAB内置了所需的数学和绘图函数
% 定义导函数
dy = @(y) y;
h = 0.1; % 求解步长
% 求解区间
a = 0;
b = 3;
% 在求解区间内间隔 h 取点作为 x 数组
x = a:h:b; % MATLAB的冒号操作符用于创建等差数列
% 定义和 x 形状相同的数组存贮数值解的值
w = zeros(size(x));
% 计算解析解
y = exp(x);
% 设置初始值
w(1) = 1;
% 开始计算
for i = 1:length(x) - 1
w(i+1) = w(i) + h * dy(w(i));
end
% 绘图
figure; % 创建一个新的图形窗口
plot(x, y, 'k-', 'LineWidth', 2, 'DisplayName', "Analytic"); % 使用DisplayName属性来设置图例标签
hold on; % 保持当前图形，以便在同一个图上绘制更多的线
plot(x, w, 'g--', 'LineWidth', 2, 'DisplayName', "Numerical");
legend('show'); % 显示图例
% saveas(gcf, './1.jpg'); % 如果你想保存图片，可以使用这个命令
% 注意：MATLAB中没有直接对应的plt.show()，因为图形窗口在脚本执行时会自动显示
2）欧拉法+改进的欧拉法
% MATLAB 代码
h = 0.1; % 步长
a = 0; % 起始点
b = 3; % 终止点
% 创建x值的数组（从a到b，包含端点b但不包括a+h）
x = a:h:b;
% 初始化权重数组（与x具有相同的形状）
w1 = zeros(size(x));
w2 = zeros(size(x));
% 解析解
y = exp(x);
% 欧拉方法
w1(1) = 1; % 初始条件
for i = 1:length(x)-1
w1(i+1) = w1(i) + h * dy_euler(x(i), w1(i)); % dy_euler是欧拉方法的导数函数
end
% 改进的欧拉方法（预测-校正）
w2(1) = 1; % 初始条件
for i = 1:length(x)-1
% 预测
w2_pred = w2(i) + h * dy_euler(x(i), w2(i));
% 校正
w2(i+1) = w2(i) + (h/2) * (dy_euler(x(i), w2(i)) + dy_euler(x(i+1), w2_pred));
end
% 绘图
figure;
plot(x, y, 'k-', 'LineWidth', 2, 'DisplayName', 'Analytic');
hold on;
plot(x, w1, 'g--', 'LineWidth', 1.5, 'DisplayName', 'Euler Method');
plot(x, w2, 'mo', 'MarkerSize', 8, 'DisplayName', 'Advanced Euler');
legend('Location', 'upperleft');
grid on;
% saveas(gcf, './3.jpg'); % 如果需要保存图像，取消注释此行
function dy = dy_euler(x, y)
dy = exp(x); % 对于指数函数 exp(x)，其导数是 exp(x)
end
3）欧拉法
function [t, y] = eulerMethod(f, y0, t0, tf, h)
% 初始化时间向量和近似解向量
t = t0:h:tf;
y = zeros(size(t));
% 初始条件
y(1) = y0;
% 欧拉法迭代
for i = 1:length(t)-1
y(i+1) = y(i) + h * f(t(i), y(i));
end
end
% 示例：求解 dy/dt = y - t^2, y(0) = 1, 在区间 [0, 2] 上，使用步长 h = 0.1
f = @(t, y) y - t.^2; % 微分方程 dy/dt = y - t^2
y0 = 1; % 初始条件 y(0) = 1
t0 = 0; % 初始时间 t0 = 0
tf = 2; % 最终时间 tf = 2
h = 0.1; % 时间步长 h = 0.1
% 调用欧拉法函数求解
[t, y] = eulerMethod(f, y0, t0, tf, h);
% 绘制结果
plot(t, y, 'o-');
xlabel('t');
ylabel('y(t)');
title('欧拉法求解 dy/dt = y - t^2');
grid on;
2 常微分方程求解(2)：龙格库塔法
泰勒展开法，用\(f\) 在同一点\((t_n,u_n)\)的高阶导数表示\(\varphi(t_n,u_n,h)\) ,这不方面数值计算，龙格库塔是用\(f\)在一些点上的值表示\(\varphi(t_n,u_n,h)\)，使得单步局部截断误差的阶和泰勒展开法相等。
把初值问题写成积分形式：
\[u(t+h)=u(t)+\int_t^{t+h}f(\tau,u(\tau) )d\tau \]在 \([t,t+h]\)取\(m\)个点，$t_1=t\le t_2\le t_3\le \cdot \cdot \cdot \le t_m \le t+h $.
若知道 \(k_i=f(t_i,u(t_i),i=1,\cdot\cdot\cdot,m,\)则用他们的一次组合去近似\(f\):
\[\sum_{i=1}^m c_i k_i =f \]怎么计算\(k_i\) 因为\(u(t_i)\) 未知，可以用欧拉法 \(u(t_2) \approx (t_1)+(t_2-t_1)f(t_1,u(t_1)k_1)\)
于是
\[k_2 \approx f(t_2,u(t_1)+(t_2-t_1)k_1) \\ \\ k_3 \approx f(t_3,u(t_1)+(t_2-t_1)k_1+(t_3-t_2)k_2) \]四阶龙格库塔法：
\[\left\{ \begin{array}{l} u_{n+1}=u_n+\frac h6(k_k+2k_2+2+k_3+k4) \\ k_1=f(t_n,u_n)\\ k_2=f(t_n + \frac 12h,u_n+\frac12 hk_1)\\ k_3=f(t_n + \frac 12h,u_n+\frac12 hk_2)\\ k_2=f(t_n + h,u_n+ hk_3) \end{array} \right. \]2.1 龙格库塔法Matlab代码
function rk4_example()
% 定义时间范围和步长
tspan = [0 2]; % 时间范围从0到2
h = 0.1; % 步长
t = tspan(1):h:tspan(2); % 生成时间向量
% 初始条件
y0 = 1;
y = zeros(size(t));
y(1) = y0; % 初始值
% 四阶龙格-库塔方法
for i = 1:(length(t)-1)
k1 = f(t(i), y(i));
k2 = f(t(i) + h/2, y(i) + h/2*k1);
k3 = f(t(i) + h/2, y(i) + h/2*k2);
k4 = f(t(i) + h, y(i) + h*k3);
y(i+1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
% 绘制结果
plot(t, y);
xlabel('Time (t)');
ylabel('y(t)');
title('Solution of dy/dt = t^2 + y using Runge-Kutta 4th order method');
grid on;
end
% 定义微分方程的函数
function dy = f(t, y)
dy = t^2 + y;
end
% 调用主函数
rk4_example();
二、 椭圆方程的有限差分法
1 二阶常微分方程边值问题
考虑最简单的椭圆方程第一边值问题
\[\begin{array}{l} Lu=-\frac{d^2u}{dx^2}+qu=f \\ u(a)=\alpha， \quad u(b)=\beta \end{array} \]其中，\(q,f\)为\([a,b]\)上的连续函数，\(q \ge 0; \quad \alpha,\beta\) 为给定常数，将区间\([a,b]\) 分成\(N\)等分，分点为
\[x_i=a+ih \quad i=0,1,\cdot\cdot\cdot,N, \]\(h=(b-a)/N\), 于是我们得到区间\(I=[a,b]\)的一个网络划分， \(x_i\)称为节点，\(h\)为步长 .
现将方程在节点\(x_i\)离散化，对充分光滑的解 \(u\), 由泰勒展开可得
\[\frac{u(x_{i+1})-2u(x_i)+u(x_{i-1})}{h^2}\\ =\bigg[\frac {d^2 u(x)}{dx^2} \bigg]_i +\frac {h^2}{12} \bigg[\frac {d^4 u(x)}{dx^4} \bigg]_i + O(h^3) \]其中\([\quad]_i\) 表示函数在\(x_i\)取值
于是在\(x_i\)可将方程写成
\[-\frac{u(x_{i+1})-2u(x_i)+u(x_{i-1})}{h^2} +q(x_i)u(x_i)=f(x_i)+R_I(u) \]其中
\[R_i(u)=-\frac {h^2}{12} \bigg[\frac {d^4 u(x)}{dx^4} \bigg]_i + O(h^3) \]当\(h→0\)，\(R_i(u)\)是 \(h\) 的二阶无穷小量，若舍去，得到逼近的差分方程
\[L_hu_i=-\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+q_iu_i=f_i, \]式中\(q_i=q(x_i),f_i=f(x_i).\) 称\(R_i(u)\)为差分方程的截断误差 . 利用差分算子\(L_h\)
$$L_hu(x_i)=f(x_i)+R_i(u).$$
而在节点\(x_i\)处，微分方程 为
$$[Lu]_i=f(x_i),$$
相减得
\[R_i=L_hu(x_i)-[Lu]_i \]\(R_i(u)\) 是用差分算子\(L_h\)逼近微分算子\(L\)引起的截断误差.
差分方程当\(i=1,2,\cdots,N-1\) 时成立，加上边值条件\(u_0=\alpha,u_N=\beta,\)
就得到关于\(u_i\)的线性方程组：
\[\begin{aligned} &L_hu_i=-\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+q_iu_i=f_i,\quad i=1,2,\cdots,N-1,\\ &u_{0}=\alpha,\quad u_{N}=\beta. \end{aligned} \]它的解\(u_i\)是\(u(x)\)在\(x=x_i\)的近似. 上述为成为差分方程，或者差分格式.
由于是二阶中心差商代替二阶微商得到，所以也教中心差分格式.
注意方程的个数等于网格内点\(x_1,x_2,\cdots,x_{N-1}\)的个数 \(N-1\),因此它是\(N-1\)阶方程组. 这是高阶方程组，但是每个方程未知数最多三个，因为系数矩阵大量元素为零. 如何把方程和未知数由左向右排列，则A是堆成三对角矩阵，例如N=5，则
\[\left. \boldsymbol{A}= \left[\begin{array}{cccc} \frac{2}{h^{2}}+q_{1}&-\frac{1}{h^{2}}&0&0\\ -\frac{1}{h^{2}}&\frac{2}{h^{2}}+q_{2}&-\frac{1}{h^{2}}&0\\ 0&-\frac{1}{h^{2}}&\frac{2}{h^{2}}+q_{3}&-\frac{1}{h^{2}}\\ 0&0&-\frac{1}{h^{2}}&\frac{2}{h^{2}}+q_{4} \end{array}\right.\right]. \]1.1 差分方程的matlab代码
% 定义方程
function f = ode_func(x, u)
% 这里q是已知参数，需要事先定义
% 假设 q = 1, 但你可以根据需要更改它
q = 1;
f = [u(2); -q*u(1) + f_given(x)];
% 假设f_given(x)是f(x)的表达式，你需要定义它
% 例如，如果f(x) = sin(x)，则定义f_given为
% function y = f_given(x)
% y = sin(x);
% end
end
% 定义边界条件
function res = bc_func(ya, yb)
% 这里alpha和beta是边界值
alpha = 0; % 假设的alpha值
beta = 1; % 假设的beta值
res = [ya(1) - alpha; yb(1) - beta];
end
% 假设的f_given函数（根据具体f(x)替换）
function y = f_given(x)
y = sin(x); % 示例：f(x) = sin(x)
end
% 主程序
% 定义区间
a = 0;
b = pi;
% 初始猜测解（需要足够好的猜测以便收敛）
solinit = bvpinit(linspace(a,b,5), [0; 0]); % 初始猜测可以是零或任何其他合理的值
% 求解BVP
sol = bvp4c(@ode_func, @bc_func, [a, b], solinit);
% 绘制解
x = linspace(a,b,100);
u = deval(sol, x);
plot(x, u(1,:));
xlabel('x');
ylabel('u(x)');
title('Solution of the Boundary Value Problem');
grid on;
2 一维差分格式
考虑两点边值问题
\[Lu= -\frac{\mathrm{d}}{\mathrm{d}x}\left(p\frac{\mathrm{d}u}{\mathrm{d}x}\right)+ r\frac{\mathrm{d}u}{\mathrm{d}x}+qu=f,\quad a<x<b,\\ u(a)=\alpha,\quad u(b)=\beta. \]假定\(p\in C^1[a,b],p(x)\geqslant p_{\min}>0,r,q,f\in C[a,b],\alpha,\beta\) 是给定的常数 .
本节构造差分格式的三种方法：
- 直接差分法
- 有限体积法
- 特定系数法
2.1 直接差分法
首先去\(N+1\)个节点：
\[a=x_0<x_1<\cdots<x_i<\cdots<x_N=b, \]将区间\(I=[a,b]\)分成\(N\)个小区间
\[I_{i}:x_{i-1}\leqslant x\leqslant x_{i},\quad i=1,2,\cdots,N. \]记\(h_i=x_i-x_i\)，称\(h=\underset{i}\max h_{i}\)为最大步长.
用\(I_h\)表示网格内点\(x_{1},x_{2},\cdots,x_{N-1}\)的集合，$\bar{I}_h \(表示内点和界点\)x_0=a,x_N=b$的集合.
取相邻节点\(x_{i-1},x_i\)的中点\(x_{i-\frac{1}{2}}=\frac{1}{2}(x_{i-1}+x_{i})(i=1,2,\cdots,N),\)称为半整点，则由节点
\[a=x_{0}<x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{i-\frac{1}{2}}<\cdots<x_{N-\frac{1}{2}}<x_{N}=b \]又构成了一个\([a,b]\)的对偶剖分
其次由差商代替微商，将方程离散化（推导是为了导出截断误差公式）
\[\begin{aligned} \frac{u(x_{i+1})-u(x_{i-1})}{h_{i}+h_{i+1}}& =\left[\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i}+\frac{h_{i+1}-h_{i}}{2}\left[\frac{\mathrm{d}^{2}u}{\mathrm{d}x^{2}}\right]_{i}+O(h^{2}), \\ \end{aligned} （1） \]\[\begin{aligned} p\left(x_{i-\frac{1}{2}}\right)\frac{u(x_{i})-u(x_{i-1})}{h_{i}}& =\left[p\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i-\frac12}+\frac{h_{i}^{2}}{24}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i-\frac12}+O(h^{3}) \\ &=\left[p\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i-\frac{1}{2}}+\frac{h_{i}^{2}}{24}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i}+O(h^{3}), \\ \end{aligned} （2） \]\[\begin{aligned} p\left(x_{i+\frac12}\right)\frac{u(x_{i+1})-u(x_{i})}{h_{i+1}}& =\left[p\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i+\frac{1}{2}}+\frac{h_{i+1}^{2}}{24}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i}+O(h^{3}). \end{aligned} （3） \]将（3）-（2）并除以 \(\frac{(h_{i}+h_{i+1})}2\)得到
\[\begin{gathered} \frac{2}{h_{i}+h_{i+1}}\left\lfloor p\left(x_{i+\frac{1}{2}}\right)\frac{u(x_{i+1})-u(x_{i})}{h_{i+1}}-p\left(x_{i-\frac{1}{2}}\right)\frac{u(x_{i})-u(x_{i-1})}{h_{i}}\right\rfloor \\ =\frac{2}{h_{i}+h_{i+1}}\left(\left[p\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i+\frac{1}{2}}-\left[p\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i-\frac{1}{2}}\right)+\frac{h_{i+1}-h_{i}}{12}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i}+O(h^{2}) \\ =\left[\frac{\mathrm{d}}{\mathrm{d}x}\left(p\frac{\mathrm{d}u}{\mathrm{d}x}\right)\right]_{i}+\frac{h_{i+1}-h_{i}}{4}\left[\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\left(p\frac{\mathrm{d}u}{\mathrm{d}x}\right)\right]_{i}+\frac{h_{i+1}-h_{i}}{12}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i}+O(h^{2}) \end{gathered} \]令\(p_{i-\frac{1}{2}}=p\left(x_{i-\frac{1}{2}}\right),r_{i}=r(x_{i}),q_{i}=q(x_{i}),f_{i}=f(x_{i}),\)
则边值问题的解\(u(x)\))满足方程
\[\begin{aligned} L_hu(x_i)& \equiv-\frac{2}{h_{i}+h_{i+1}}\left[p_{i+\frac{1}{2}}\frac{u(x_{i+1})-u(x_{i})}{h_{i+1}}-p_{i-\frac{1}{2}}\frac{u(x_{i})-u(x_{i-1})}{h_{i}}\right]+ \\ &\frac{r_i}{h_i+h_{i+1}}\left[u(x_{i+1})-u(x_{i-1})\right]+q_iu(x_i) \\ &=f_i+R_i(u), \end{aligned} \]其中
\[R_{i}(u)=-(h_{i+1}-h_{i})\left(\frac{1}{4}\left[\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\left(p\frac{\mathrm{d}u}{\mathrm{d}x}\right)\right]_{i}+\frac{1}{12}\left[p\frac{\mathrm{d}^{3}u}{\mathrm{d}x^{3}}\right]_{i}-\frac{1}{2}\left[r\frac{\mathrm{d}^{2}u}{\mathrm{d}x^{2}}\right]_{i}\right)+O(h^{2}) \]\(R_{i}(u)\)为差分算子的截断误差，舍去\(R_{i}(u)\),则边值问题的差分方程为
\[\begin{aligned} L_{h}u_{i}& \equiv-\frac{2}{h_{i}+h_{i+1}}\left[p_{i+\frac{1}{2}}\frac{u(x_{i+1})-u(x_{i})}{h_{i+1}}-p_{i-\frac{1}{2}}\frac{u(x_{i})-u(x_{i-1})}{h_{i}}\right]+ \\ &\frac{r_i}{h_i+h_{i+1}}(u_{i+1}-u_{i-1})+q_iu_i=f_i,\quad i=1,\cdots,N-1, \\ &u_0=\alpha,\quad u_N=\beta. \end{aligned} \]差分方程也可以数值微商公式
\[\begin{aligned} \left[\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i}& \approx\frac{u_{i+1}-u_{i-1}}{h_{i}+h_{i+1}}, \\ \left[\frac{\mathrm{d}}{\mathrm{d}x}\left(p\frac{\mathrm{d}u}{\mathrm{d}x}\right)\right]_{i}& \approx\left(p_{i+\frac{1}{2}}\left[\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i+\frac{1}{2}}-p_{i-\frac{1}{2}}\left[\frac{\mathrm{d}u}{\mathrm{d}x}\right]_{i-\frac{1}{2}}\right)/\frac{h_{i}+h_{i+1}}{2} \\ &\approx\frac{2}{h_{i}+h_{i+1}}\left(p_{i+\frac{1}{2}}\frac{u_{i+1}-u_{i}}{h_{i+1}}-p_{i-\frac{1}{2}}\frac{u_{i}-u_{i-1}}{h_{i}}\right) \end{aligned} \]前面的推导是为了导出截断误差公式。
截断误差可表示为
\[\begin{aligned} L_hu_i& =-\frac1{h^2}\left[p_{i+\frac12}u_{i+1}-(p_{i+\frac12}+p_{i-\frac12})u_i+p_{i-\frac12}u_{i-1}\right]+ \\ &r_i\frac{u_{i+1}-u_{i-1}}{2h}+q_iu_i=f_i. \end{aligned} \]差分方程可以化简为
\[\begin{aligned} L_hu_i& =-\frac1{h^2}\left[p_{i+\frac12}u_{i+1}-(p_{i+\frac12}+p_{i-\frac12})u_i+p_{i-\frac12}u_{i-1}\right]+ \\ &r_i\frac{u_{i+1}-u_{i-1}}{2h}+q_iu_i=f_i. \end{aligned} \]2.2 有限体积法（1）
考虑守恒性微分方程
\[Lu=-\frac{\mathrm{d}}{\mathrm{d}x}\left(p(x)\frac{\mathrm{d}u}{\mathrm{d}x}\right)+q(x)u=f(x). \]把他看成分布在一根杆上的稳定温度场方程，则在 \([a,b]\)内人一小区间 \([x^{(1)},x^{(2)}]\)上的热量守恒定律具有形式
\[-\int_{x^{(1)}}^{x^{(2)}}\frac{\mathrm{d}}{\mathrm{d}x}\left(p(x)\frac{\mathrm{d}u}{\mathrm{d}x}\right)\mathrm{d}x+\int_{x^{(1)}}^{x^{(2)}}qu\mathrm{d}x=\int_{x^{(1)}}^{x^{(2)}}f\mathrm{d}x, \]或
\[W(x^{(1)})-W(x^{(2)})+\int_{x^{(1)}}^{x^{(2)}}q(x)u\mathrm{d}x=\int_{x^{(1)}}^{x^{(2)}}f\mathrm{d}x, \]其中
\[W(x)=p(x)\frac{\mathrm{d}u}{\mathrm{d}x}. \]把微分方程写成积分守恒形式后，最高阶微商由二阶降为一阶，从而减弱了对\(p、u\)光滑性要求，从积分守恒型方程出发构造差分格式，便于推广到任意网格和处理第二边值条件.
既然具有守恒形式的微分方程放映了物理、力学某些守恒定律，那么构造的差分格式也反应了这一基本形式. 现在来构造差分格式.
取\([x^{(1)},x^{(2)}]\)为对偶单元 $$[x_{i-\frac12},x_{i+\frac12}]$$ ,则
\[W\left(x_{i-\frac{1}{2}}\right)-W\left(x_{i+\frac{1}{2}}\right)+\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}qu\mathrm{d}x=\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f\mathrm{d}x. \]考虑到\(p(x)\) 可能有间断点，进一步差分是不合适，但是“热流量”\(W(X)\) 恒连续，因而
\[\frac{\mathrm{d}u}{\mathrm{d}x}=\frac{W(x)}{p(x)}, \]再沿\([x^{(1)},x^{(2)}]\)积分，得到
\[u_i-u_{i-1}=\int_{x_{i-1}}^{x_i}\frac{W(x)}{p(x)}\mathrm{d}x, \]利用中矩形公式
\[W_{i-\frac{1}{2}}\approx a_{i}\frac{u_{i}-u_{i-1}}{h_{i}},\\a_{i}=\left[\frac1{h_i}\int_{x_{i-1}}^{x_i}\frac{\mathrm{d}x}{p(x)}\right]^{-1}. \]又
\[\begin{gathered} \int_{x_{i-\frac12}}^{x_{i+\frac12}}qu\mathrm{d}x\approx\frac{h_{i}+h_{i+1}}2d_{i}u_{i}, \\ d_{i}=\frac{2}{h_{i}+h_{i+1}}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}q(x)\mathrm{d}x. \end{gathered} \]分别代入，得到守恒型差分方程：
\[-\left[a_{i+1}\frac{u_{i+1}-u_{i}}{h_{i+1}}-a_{i}\frac{u_{i}-u_{i-1}}{h_{i}}\right]+\frac{1}{2}(h_{i}+h_{i+1})d_{i}u_{i}=\frac{1}{2}(h_{i}+h_{i+1})\phi_{i},\\\phi_{i}=\frac{2}{h_{i}+h_{i+1}}\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}f(x)\mathrm{d}x. \]如果系数 \(p,q\)及右端\(f\)光滑，则利用中矩形公式计算
\[\begin{cases}a_i=p_{i-\frac12}=p\left(x_{i-\frac12}\right),\\d_i=q_i=q(x_i),\\\phi_i=f_i=f(x_i),\end{cases} \]2.2.1 有限体积法（2）：matlab代码
% 参数定义
u_inf = 1; % 无穷远处速度
L = 1; % 求解区域长度
Nx = 100; % 网格数
dx = L / (Nx - 1); % 网格间距
t_final = 1; % 模拟时间
CFL = 0.5; % CFL数，用于控制时间步长
u0 = zeros(1, Nx); % 初始条件，这里假设全为0
u0(Nx//2) = u_inf; % 假设在中心位置有一个初始扰动
% 时间步长计算
dt = CFL * dx / u_inf; % 使用CFL条件确保数值稳定性
nt = ceil(t_final / dt); % 总时间步数
% 有限体积法求解
u = u0; % 初始化u
for n = 1:nt
un = u; % 保存当前时间步的解用于计算下一时间步
for i = 2:Nx-1 % 忽略边界点（需要特殊处理）
u(i) = un(i) - dt/dx * (un(i) * u_inf - un(i-1) * u_inf); % 一维对流方程的离散化形式
end
% （这里省略了边界条件的处理，通常需要根据具体问题设置）
% ...
% 更新时间
t = n * dt;
% （可选）可视化或其他后处理步骤
% ...
end
% 可视化结果
x = linspace(0, L, Nx);
plot(x, u);
xlabel('x');
ylabel('u(x, t_final)');
title(['1D Convection Equation Solution at t = ' num2str(t_final)]);
2.3 待定系数法
略
3 矩阵网的差分格式
考虑泊松方程
\[-\Delta u=f(x,y),\quad(x,y)\in G, \]G是x,y平面的有界区域，其边界\(\Gamma\)为分段光滑曲线，在\(\Gamma\)上\(u\)满足下列边界条件之一：
\[\begin{aligned} &u\Big|_{\Gamma}=\alpha(x,y) \qquad (\text{第一边值条件}), \\ &\left.\frac{\partial u}{\partial n}\right|_{\Gamma}=\beta(x,y) \quad(\text{第二边值条件}), \\ &\frac{\partial u}{\partial n}+ku\Big|_{\Gamma}=\gamma(x,y)(\text{第三边值条件}), \end{aligned} \]3.1 泊松方程的matlab代码
% 参数设置
N = 21; % 网格大小（例如，21x21的网格）
h = 1 / (N - 1); % 网格步长
x = 0:h:1; % x坐标
y = 0:h:1; % y坐标
[X, Y] = meshgrid(x, y); % 创建网格
% 初始化u，包括边界条件
u = zeros(N); % 初始化一个N x N的矩阵，但这里我们先处理一维边界
u_boundary = @(x, y) sin(pi*x).*sin(pi*y); % 假设的边界条件函数
u(1, :) = u_boundary(0, y); % x=0边界
u(end, :) = u_boundary(1, y); % x=1边界
u(:, 1) = u_boundary(x, 0); % y=0边界
u(:, end) = u_boundary(x, 1); % y=1边界
u = u(:); % 将二维矩阵转换为一维向量，以便进行迭代
% 创建一个用于存储中间结果的向量
u_old = u;
% 迭代参数
max_iter = 10000; % 最大迭代次数
tol = 1e-6; % 收敛容差
% 迭代求解
for iter = 1:max_iter
u_old = u; % 保存上一次的解
% 有限差分格式（五点中心差分）
for i = 2:N-1
for j = 2:N-1
i_idx = (i-1)*N + j; % 一维索引
u(i_idx) = 0.25 * (u_old((i-1)*N + j-1) + u_old((i-1)*N + j+1) + ...
u_old((i+1)*N + j-1) + u_old((i+1)*N + j+1) - h^2 * f(X(i,j), Y(i,j)));
end
end
% 应用边界条件
u(1:N) = u_boundary(0, y); % x=0边界
u((N-1)*N+1:end) = u_boundary(1, y); % x=1边界
u(1:N:end) = u_boundary(x, 0); % y=0边界
u(N:N:end) = u_boundary(x, 1); % y=1边界
% 检查收敛性
if norm(u - u_old, inf) < tol
fprintf('Converged in %d iterations\n', iter);
break;
end
end
% 如果需要，将解转换回二维矩阵形式
u_matrix = reshape(u, N, N);
% 可视化解
surf(X, Y, u_matrix);
colorbar;
title('Solution of the 2D Poisson Equation');
xlabel('x');
ylabel('y');
zlabel('u(x, y)');