一维扩散方程的有限差分法
——计算物理实验作业七
陈万 物理学2013级 13020011006
? 题目：
编程求解一维扩散方程的解
??u?2u??D2(0?x?a0,0?t?tmax)?x??t?u(x,t)|t?0?ex???u au?b?c(x?0)?111?n??u?au?b?c2(x?a0)22??n?取a1?1,b1?1,c1?0,a2?1,b2??1,c2?0,a0?1.0,tmax?10,D?0.1,h?0.1,??0.1。输出t=1,2,...,10时刻的x和u(x),并与解析解u=exp(x+0.1t)作比较。
? 主程序：
% 一维扩散方程的有限差分法 clear,clc;
%定义初始常量
a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0; a0 = 1.0; t_max = 10; D = 0.1; h = 0.1; tao = 0.1;
%调用扩散方程子函数求解
u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);
? 子程序1：
function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)
% 一维扩散方程的有限差分法，采用隐式六点差分格式(Crank-Nicolson) % a0: x的最大值 % t:_max: t的最大值 % h: 空间步长 % tao: 时间步长 % D：扩散系数
% a1,b1,c1是（x=0）边界条件的系数；a2,b2,c2是（x=a0）边界条件的系数
x = 0:h:a0;
1
n = length(x); t = 0:tao:t_max; k = length(t);
P = tao * D/h^2; P1 = 1/P + 1; P2 = 1/P - 1;
u = zeros(k,n); %初始条件 u(1,:) = exp(x);
%求A矩阵的对角元素d d = zeros(1,n);
d(1,1) = b1*P1+h*a1; d(2:(n-1),1) = 2*P1; d(n,1) = b2*P1+h*a2;
%求A矩阵的对角元素下面一行元素e e = -ones(1,n-1); e(1,n-1) = -b2;
%求A矩阵的对角元素上面一行元素f f = -ones(1,n-1); f(1,1) = -b1;
R = zeros(k,n);%求R %追赶法求解 for i = 2:k
R(i,1) = (b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1; for j = 2:n-1
R(i,j) = u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1); end
R(i,n) = b2*u(i-1,n-1)+( b2*P2-h*a2)*u(i-1,n)+2*h*c2; M = chase(e,d,f,R(i,:)); u(i,:) = M';
plot(x,u(i,:)); axis([0 a0 0 t_max]); pause(0.1) end
output = u;
% 绘图比较解析解和有限差分解 [X,T] = meshgrid(x,t); Z = exp(X+0.1*T);
surf(X,T,Z),xlabel('x'),ylabel('t'),zlabel('u'),title('解析解'); figure
2
surf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title('有限差分解');
? 子程序2：
function M = chase(a,b,c,f)
% 追赶法求解三对角矩阵方程，Ax=f % a是对角线下边一行的元素 % b是对角线元素
% c是对角线上边一行的元素
% M是求得的结果，以列向量形式保存 n = length(b);
beta = ones(1,n-1); y = ones(1,n); M = ones(n,1);
for i = (n-1):(-1):1 a(i+1) = a(i); end
% 将a矩阵和n对应
beta(1) = c(1)/b(1); for i = 2:(n-1)
beta(i) = c(i)/( b(i)-a(i)*beta(i-1) ); end
y(1) = f(1)/b(1); for i = 2:n
y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end
M(n) = y(n);
for i = (n-1):(-1):1
M(i) = y(i)-beta(i)*M(i+1); end end
? 结果：
3