热传导方程差分格式

R = zeros(N-1 , N-1);
fori = 2 : N-2
R(i ,i) = 1 + 2 * r;
R(i , i+1) = -r;
R(i, i-1) = -r;
end
R(1 , 1) = 1 + 2 * r;
R(1 , 2) = -r;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
向前差分格式：
：
:
:
:
:
:
:
向后差分格式：
六点差分格式：
误差：
t=1/1600 t=1/3200 t=1/6400
向前差分格式
t=0.05
8.888396579076e+21
0.0014
0.00034640856587426
t=0.1
8.83785296707723e+59
ቤተ መጻሕፍቲ ባይዱ0.0017
0.000422936658240724
t=0.1
0.0443229427820120
0.0442555796634600
0.0442226561119014
t=0.2
0.0265375944429049
0.0264980178171388
0.0264788945621501
3方法总结及分析
本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。根据matlab作，特别明显的是， 时，图像解析解波动特别大，与真解差距很大。这是因为 差分格式不稳定。根据误差比较，解本题时向后差分格式误差最小。衡量一个差分格式是否经济实用，有点多方面因素决定，应考虑不同的情况决定。
附件程序
向前差分格式：
function[ u , err ] = xq( t , t1 )
% t 是时刻值；
% t1 是时间步长；
N = 40;
h = 1 / N;
x=[h : h : (1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
R = zeros(N-1 , N-1);
fori = 2 : N-2
R(i , i) = 1 - 2 * r;
R(i , i+1) = r;
R(i, i-1) = r;
end
R(1 , 1) = 1 -2 * r;
R(1 , 2) = r;
R(N-1 , N-1) = 1 - 2 * r;
R(N-1, N-2) = r;
k = t / t1;
legend('数值解','精确解');
err=norm(u-up(:,k));
end
fori = 1 : k
u(:,i+1) =inv(R) *R1* u(:,i);
end
u=[0;u(:,i+1);0];
fori=1:k
forj=1:N+1
up(j,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u,'b.',x,up(:,k),'r--');
已知 在相应区域光滑，并且在 与边值相容，使问题有唯一充分光滑的解。
取空间步长 ，和时间步长 ，其中 都是正整数。用两族平行直线 和 将矩形域 分割成矩形网络，网络格节点为 。以 表示网格内点集合，即位于矩形 的网点集合； 表示闭矩形 的网格集合； 是网格界点的集合。
向前差分格式，即
（1）
，
其中， 以 表示网比。（1）式可改写成如下：
此格式为显格式。
其矩阵表达式如下：
（2）向后差分格式
向后差分格式，即
（2）
其中 （2）式可改写成
此种差分格式被称为隐格式。
其矩阵表达式如下：
（3）六点对称格式
六点差分格式：
（3）
将（3）式改写成
其矩阵表达式如下:
2利用MATLAB求解问题的过程
对每种差分格式依次取 ， ， ， ，用MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的 误差。
fori = 1 : k
u(: , i+1) =R * u(: , i);
end
u=[0 ; u(: , i+1) ; 0];
fori=1 : k
forj=1 : N+1
up(j , i)=exp(-pi*pi*t) * sin(pi*(j-1)*h);
end
end
x=0 : h : 1;
plot(x , u ,'b.', x , up(: , k) ,'r--');
legend('数值解','精确解');
err=norm(u - up(: , k));
end
向后差分格式：
function[ u , err ] = xh( t , t1 )
N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
R = zeros(N-1 , N-1);
fori = 2 : N-2
R(i ,i) = 1 + r;
R(i , i+1) = -r/2;
R(i, i-1) = -r/2;
end
R(1 , 1) = 1 + r;
R(1 , 2) = -r/2;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
legend('数值解','精确解');
err=norm(u-up(:,k));
end
六点差分格式：
function[ u , err ] =ld( t , t1 )
N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
0.00337797222842041
0.00211264169361075
t=0.2
0.00440853027126062
0.00252054496720081
0.00157579425393478
六点差分格式
t=0.05
0.0526000555873324
0.0525169033137740
0.0524758656498487
一维抛物方程的初边值问题
分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：
在 时刻的数值解，并与解析解 。
1差分格式形式
设空间步长 , 时间步长 ， ，网比 .
（1）向前差分格式
该问题是第二类初边值问题（混合问题），我们要求出所需次数的偏微商的函数 ，满足方程 和初始条件 ，及边值条件 。
t=0.2
1.16572934157246e+136
0.00126149089942459
0.000315223617970853
向后差分格式
t=0.05
0.00483090554361983
0.00276517069125328
0.00172971295195573
t=0.1
0.00590373504606959
fori = 2 : N-2
R1(i ,i) = 1 - r;
R1(i , i+1) = r/2;
R1(i, i-1) = r/2;
end
R1(1 , 1) = 1 - r;
R1(1 , 2) = r/2;
R1(N-1 , N-1) = 1 -2 * r;
R1(N-1, N-2) = r/2;
k = t / t1;
k = t / t1;
fori = 1 : k
u(:,i+1) =inv(R) * u(:,i);
end
u=[0;u(:,i+1);0];
fori=1:k
forj=1:N+1
up(j,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u,'b.',x,up(:,k),'r--');