一維熱傳到方程求數值解
本文主要利用泰勒展開將方程中的一階還有二階偏導數進行離散化，推導出一種可以用程序求解的形式
求解原理
一維熱傳導方程
\[\begin{align} \begin{cases} \frac{\partial u}{\partial x} \left ( x,t \right ) &=a^{2} \frac{\partial^{2}u}{\partial x^2}u(x,t)+f(x,t)\\ u(x,0)&=\varphi({x})\\ u(a,t)&=\gamma_{1}(t)\\ u(b,t)&=\gamma_2(t) \end{cases} \end{align} \]
由於熱傳導方程較為復雜，只能將方程中的一階和二階偏導進行離散化。和歐拉法采用相同的思路，下面進行推導：
-
將\(x\)與\(t\)分別在橫坐標與縱坐標上進行划分
\(x\)步長: \(\Delta{x}= \frac{b-a}{N}\)，得到關於\(x_j\)與\(t_n\)的表達式：
\[\begin{aligned} x_j &= a + (j-1)\Delta{x} \\ t_n &= 0 + (n-1)\Delta{t} \\ \end{aligned} \]
將函數進行近似替換\(u_j^n\approx u(x_j,t_n)\)
-
根據泰勒展開將公式進行代換
對於任意一個\(x_j\)對\(t\)進行展開：
\[u(x_j,t_n+\Delta{t})=u(x_j,t_n)+\frac{\partial u}{\partial t} (x_j,t_n)\Delta{t}+··· \]由於很難求出函數的偏導，所以需要將其所有偏導形式轉換成容易求解出來的離散形式首先用一維熱傳導方程進行替換
\[\frac{\partial u}{\partial t} (x_j,t_n) = a^2 \frac{\partial^{2}u}{\partial x^2}(x_j,t_n)+f(x_j,t_n) \]利用上式聯立下面兩個式子
\[\begin{aligned} u(x_j+\Delta{x},t_n) &= u(x_j,t_n)+\frac{\partial u}{\partial x} (x_j,t_n)\Delta{x}+1/2\frac{\partial^{2}u}{\partial x^2}(x_j,t_n)\Delta{x}^2+···\\ u(x_j-\Delta{x},t_n) &= u(x_j,t_n)-\frac{\partial u}{\partial x} (x_j,t_n)\Delta{x}+1/2\frac{\partial^{2}u}{\partial x^2}(x_j,t_n)\Delta{x}^2-··· \\ \frac{\partial^{2}u}{\partial x^2}(x_j,t_n) &\approx \frac{u_{j+1}^n+u_{j-1}^n-2u_j^n}{\Delta{x}^2} \end{aligned} \]最后得到遞推關系式
\[u_j^{n+1}=u_j^n+[a^2\frac{u_{j+1}^n+u_{j-1}^n-2u_j^n}{\Delta{x}^2}+f_j^n]\Delta{t} \]
化成易於用程序求解的形式
-
在時間維度上進行遞推
首先設置兩個時間向量，將所有的位置包括其中
\[u^n= \begin{pmatrix}u_1^n \\\vdots \\\vdots \\u_{N+1}^n \end{pmatrix}\qquad u^{n+1}= \begin{pmatrix}u_1^{n+1} \\\vdots \\\vdots \\u_{N+1}^{n+1} \end{pmatrix} \] -
建立系數矩陣
\[\begin{pmatrix} \phi \\ 第一取值 \\ \vdots \\ 第N取值 \\ \phi \end{pmatrix}= \begin{pmatrix}-2\quad1\quad0\quad0\quad··· \\1\quad-2\quad1\quad0\quad··· \\\vdots \\···\quad0\quad1\quad-2\quad1 \\···\quad0\quad0\quad1\quad-2 \end{pmatrix} \begin{pmatrix}u_1^n \\u_2^n \\\vdots \\\vdots \\u_{N+1}^n \end{pmatrix} \]為何矩陣要這么建立，系數矩陣\(A\)的第二行為例，與右邊的列向量相乘得到結果\(u_1^n - 2u_2^n + u_3^n\) 將結果表示成以下列向量。由於初值和末值原方程會給出，所以開始兩個值可能不正確，所以先不用管。
\[\begin{pmatrix} \phi \\ u_1^n - 2u_2^n + u_3^n \\ u_2^n - 2u_3^n + u_4^n \\ \vdots \\ u_{N - 1}^n - 2u_N^n + u_{N + 1}^n \\ \phi \\ \end{pmatrix} \] -
擴散源或熱源表示
\[f^n= \begin{pmatrix}f_1^n \\\vdots \\\vdots \\f_{N+1}^n \end{pmatrix} \] -
得到最終表達式
\[u^{n+1}=u^n+(a^2\frac{1}{\Delta{x}^{2}}Au^n+f^n)\Delta{t}\\ \]以下為初始條件：
\[\begin{aligned} u_1^{n+1} &= \gamma_1^{n+1} = \gamma(t_{n+1}) \\ u_{N+1}^{n+1} &= \gamma_2^{n+1} \\ u^1 &= \gamma \\ \gamma_j &= \gamma(x_j) \\ \end{aligned} \]通過程序可以設定步長\(\Delta x\) 和 \(\Delta t\) 其它量就可以通過邊界條件推出。
示例與代碼
設方程如下：
\[\begin{align} \begin{cases} \frac{\partial u}{\partial x} \left ( x,t \right ) &=a^{2} \frac{\partial^{2}u}{\partial x^2}u(x,t)+f(x,t)\\ u(x,0) &= 0\\ u(a,t) &= 0 + 0.0 \times sin(t)\\ u(b,t) &= 0 - 0.0 \times sin(t)\\ \end{cases} \end{align} \]
求解\(matlab\)代碼如下：
clc; clear
a = 1;
dx = 0.01;
x = 0 : dx : 1;
dt = 0.00005;
t = 0 : dt : 1;
u = zeros(length(x), length(t));
u(:,1) = 0;
f = 5 * exp(-20 * (x - 1 / 2) .^ 2);
m1 = 0 + 0.0 * sin(t);
m2 = 0 - 0.0 * sin(t);
A = -2 * eye(length(x)) + diag(ones(1,length(x)-1),1) + diag(ones(1,length(x)-1),-1)
for n = 1 : length(t)-1
u(:, n+1) = u(:,n) + (a^2 / dx^2 * A * u(:,n) + f') * dt;
u(1, n+1) = m1(n+1);
u(end, n+1) = m2(n+1);
plot(x, u(:,end))
axis([x(1) x(end) 0 1])
getframe;
end