一维热传导方程求数值解
一维热传到方程求数值解
本文主要利用泰勒展开将方程中的一阶还有二阶偏导数进行离散化，推导出一种可以用程序求解的形式
求解原理
一维热传导方程
\[\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
