目錄
目錄
- 目錄
- 參考文獻
- 1 李榮華 微分方程的數值解法
- 零、 偏微分方程
- 偏微分方程形式
- 拋物線方程：擴散方程
- 雙曲線方程：對流方程
- 橢圓方程(1):拉普拉斯方程
- 橢圓方程(2):泊松方程
- 偏微分方程形式
- 一、常微分方程初值問題
- 常微分方程形式
- 1 常微分方程求解(1)：尤拉法
- 1.1 尤拉法與改進尤拉法matlab程式碼
- 1）尤拉法
- 2）尤拉法+改進的尤拉法
- 3）尤拉法
- 1.1 尤拉法與改進尤拉法matlab程式碼
- 2 常微分方程求解(2)：龍格庫塔法
- 2.1 龍格庫塔法Matlab程式碼
- 二、 橢圓方程的有限差分法
- 1 二階常微分方程邊值問題
- 1.1 差分方程的matlab程式碼
- 2 一維差分格式
- 2.1 直接差分法
- 2.2 有限體積法（1）
- 2.2.1 有限體積法（2）：matlab程式碼
- 2.3 待定係數法
- 3 矩陣網的差分格式
- 3.1 泊松方程的matlab程式碼
- 1 二階常微分方程邊值問題
- 三、拋物線方程的有限差分法
- 1 最簡單的差分格式
參考文獻
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)');
三、拋物線方程的有限差分法
1 最簡單的差分格式
\[\frac{\partial u }{\partial t}=a\frac{\partial^2 u }{\partial x^2}+f(X),
\quad 0<t<T,
\]
\[\int _{-\infty}^{+\infty} e^{(-x^2)}dx
\]