一维扩散方程的差分法matlab
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
一维扩散方程的有限差分法
——计算物理实验作业七
陈万
题目:
编程求解一维扩散方程的解
取1.0,1.0,1.0,10,0.1,0,1,1,0,1,1max 0222111======-=====τh D t a c b a c b a 。输出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;
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
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
结果:
对比分析两图,结果令人满意。
取出t_max 时刻的u 值分析:(111~u u )
有限差分解如下:
3.004166024
解析解如下:
分析数据可知误差量级为)(O )(O 22h +τ=0.12+0.12=0.02
总结:
(1)隐式六点差分格式(Crank-Nicolson)基本思想是用前一时刻的三个点表示后一时刻的三个点。因为不是直接表示u(k+1,i),故称为隐式差分。
(2)x 由等步长被分割为N 个点,列出N 个方程。采用追赶法求解得到结果。原理很简单,关键是求解AU=R 中的A 和R 。
(3)子函数2的功能是实现追赶法,该程序中没有直接用A 来表示三对角矩阵,而是把3列对角元素直接拿出来,因此在调用时应当注意各对角元素的位置,避免调用错误。