0%

一维热传导方程的差分格式

一维热传导方程差分方法的 MATLAB 编程实现

一维热传导方程

其中 是正常数, 为已知函数。 为初始条件, 为边界条件。

向前差分格式

以空间步长 、时间步长 分别将 轴上区间 轴上区间 分成 等分,可得

一维热传导方程的向前差分格式

其中 . 以 表示网比

以上格式可改写差分格式

先取 ,利用 和边值 算出第一层值 ,再取 ,利用 和边值便可算出 。如此下去,便可求出所有

差分格式写成矩阵的形式:

其中矩阵 、向量 的定义如下,注意向量 的首尾元素已包含了边界条件。

数值例子

方程的真解:.

计算数值解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
% fdm_heat.m
% forward difference scheme for heat equation
% u_t=u_{xx}, (x,t) in (0,1)x(0,1],
% u(x,0)=exp(x), x in [0,1],
% u(0,t)=exp(t), u(1,t)=exp(1+t), t in (0,1]
% exact solution: u(x,t)=exp(x+t)
clear all; close all;
a=1;
h=0.05; x=[0:h:1];
tau=0.00125; t=[0:tau:1];
r=a*tau/h^2;
M=length(x)-1; N=length(t)-1;
% constructing the coefficient matrix
e=r*ones(M-1,1);
A=spdiags([e 1-2*e e],[-1 0 1],M-1,M-1);
% setting initial and boundary conditions
u=zeros(M+1,N+1);
u(:,1)=exp(x);
u(1,:)=exp(t);
u(end,:)=exp(1+t);
for n=1:N
u(2:M,n+1)=A*u(2:M,n);
u(2,n+1)=u(2,n+1)+r*u(1,n);
u(M,n+1)=u(M,n+1)+r*u(end,n);
end
% plot the figure
mesh(t(1:20:end),x,u(:,1:20:end))
set(gca,'fontsize',12)
xlabel('t','fontsize', 14)
ylabel('x','fontsize',14)
zlabel('u','fontsize',14)

% calculating maximum error
[T X]=meshgrid(t,x);
ue=exp(X+T);
Error=max(max(abs(ue-u)))

% print -dpng -r600 fdm_heat.png
% print -depsc2 fdm_heat.eps

输出结果

1
2
Error =
2.1748e-04

fdm_heat.png