數學建模之差分演算法（求解偏微分方程）
差分演算法(求解偏微分方程)
定義
差分方法又稱為有限差分方法或網格法，是求偏微分方程定解問題的數值解中應用 最廣泛的方法之一。它的基本思想是：先對求解區域作網格剖分，將自變數的連續變化 區域用有限離散點（網格點）集代替；將問題中出現的連續變數的函式用定義在網格點 上離散變數的函式代替；通過用網格點上函式的差商代替導數，將含連續變數的偏微分 方程定解問題化成只含有限個未知數的代數方程組（稱為差分格式）。如果差分格式有 解，且當網格無限變小時其解收斂於原微分方程定解問題的解，則差分格式的解就作為 原問題的近似解（數值解）。因此，用差分方法求偏微分方程定解問題一般需要解決一下問題：
- 選取網路；
- 對微分方程及定解條件選擇差分近似，列出差分格式；
- 求解差分格式；
- 討論差分格式解對於微分方程解的收斂性及誤差估計
演算法詳解
因此，只要確定了步長，我們就可以將連續變化的自變數用有限離散點來表示
對於(9-3)的式子，為了方便計算，我們用差分來表示偏微分方程
\[\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y) \]由於式子中進行了兩次偏導，因此我們一步一步進行分析。
進行一次差分，我們用向前差分。
\[\frac{\partial u}{\partial x}=\frac{u(k+1,j)-u(k,j)}{h} \]由於向前差分有誤差，如果我們進行兩次向前差分的話，計算的誤差可能會增大，因此，第二次偏導我們選擇向後差分。即我們混合向前差分、向後差分來近似代替兩次偏導。
因此，第二次我們用向後差分
\[\begin{aligned} \frac{\partial^2u}{\partial x^2}&=\frac{\partial}{\partial x}(\frac{u(k+1,j)-u(k,j)}{h})\\ & = \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})-\frac{\partial}{\partial x}(\frac{u(k,j)}{h}) \end{aligned} \]\[\begin{aligned} \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})&=\frac{\frac{u(k+1,j)-u(k,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k+1,j)-u(k,j)}{h^2} \end{aligned} \]\[\begin{aligned} \frac{\partial}{\partial x}(\frac{u(k,j)}{h})&=\frac{\frac{u(k,j)-u(k-1,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k,j)-u(k-1,j)}{h^2} \end{aligned} \]綜上，
同理可得
\[\frac{\partial^2u}{\partial y^2}=\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2} \]所以原方程就變為
\[\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y)\\ \frac{u(k+1,j)-2u(k,j)+u(k-1,j)}{h^2}+\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2}=f(x,y) \]以上的進行通過混合向前差分、向後差分的方法求二階偏導的方法實際上就是二階中心差分
常見的差分
向前差分
函式的前向差分通常簡稱為函式的差分。對於函式f(x)，如果在等距節點：
\[x_k=x_0+kh(k=0,1,\cdots,n) \]\[\Delta f(x_k)=f(x_{k+1})-f(x_k) \]則稱Δf(x)，函式在每個小區間上的增量$$y(k+1)-y(k)$$為f(x)的一階前向差分。在微積分學中的有限差分（finite differences），前向差分通常是微分在離散的函式中的等效運算。差分方程的解法也與微分方程的解法相似。當是多項式時，前向差分為Delta運算元，一種線性運算元。前向差分會將多項式階數降低1
向後差分
對於函式\(f(x_k)\)，一階向後差分為：
\[\Delta f(x_k)=f(x_k)-f(x_{k-1}) \]中心差分
對於函式\(f(x_k)\)，一階中心差分為:
\[\Delta f(x_k)=\frac{1}{2}(f(x_{k+1})-f(x_k)) \]例題
用五點菱形個格式求解Laplace方程第一邊值問題
\[\begin{cases} \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=0,&(x,y)\in \Omega\\ u(x,y)|_{(x,y)\in \Gamma}=\lg [(1+x)^2+y^2],&\Gamma=\partial \Omega, \end{cases} \]其中\(\Omega=\left\{ (x,y)|0 \leq x,y \leq 1 \right\}\)，取\(h=\tau=\frac{1}{3}\)
解：
根據題意，我們畫出網格出來，正好構成四個五點菱形，即得到四個方程，我們將線性方程組寫出來
\[\begin{cases} \frac{1}{h^2}(u_{1,2}+u_{2,1}+u_{1,0}+u_{0,1}-4u(1,1))=0\\ \frac{1}{h^2}(u_{2,2}+u_{3,1}+u_{2,0}+u_{1,1}-4u(2,1))=0\\ \frac{1}{h^2}(u_{1,3}+u_{2,2}+u_{1,1}+u_{0,2}-4u(1,2))=0\\ \frac{1}{h^2}(u_{2,3}+u_{3,2}+u_{2,1}+u_{1,2}-4u(2,2))=0 \end{cases} \]線性方程組又可以化簡成
\[\begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} \]解非齊次線性方程組求得：
\[u_1=0.6348, u_2=1.06, u_3=0.7985, u_4=1.1698 \]計算的MATLAB程式如下
clc;
clear;
f1 = @(x)2 * log(1+x);
f2 = @(x)log((1+x).^2+1);
f3 = @(y)log(1+y.^2);
f4= @(y)log(4+y.^2);
u=zeros(4);
m=4;% 總列數
n=4;% 總行數
h=1/3;
u(1,1:m)=feval(f3,0:h:(m-1)*h)';
u(n,1:m)=feval(f4,0:h:(m-1)*h)';
u(1:n,1)=feval(f1,0:h:(n-1)*h);
u(1:n,m)=feval(f2,0:h:(n-1)*h);
b = -[u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3)];
a = [-4,1,1,0;1,-4,0,1;1,0,-4,1;0,1,1,-4];
x=a\b;
程式碼解讀
觀察線性方程：
\[\begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} \]這個形式為最小二乘法的形式
參考Matlab中的線性最小二乘的標準型：
\[\min _A \Vert RA-Y \Vert_2^2 \]所以我們只需要求出等式右邊的式子，那麼我們就可以解出等式中間的向量
現在我們來解等式右邊的向量，觀察可知，為定義域的邊界
注意這裡從1開始，與上式子有些不同
\[\begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\ u_{2,1} & u_{2,2} & u_{2,3} & u_{2,4}\\ u_{3,1} & u_{3,2} & u_{3,3} & u_{3,4}\\ u_{4,1} & u_{4,2} & u_{4,3} & u_{4,4} \end{bmatrix} \]因此我們只需要解出邊界的值，就可以把b表示出來