一维扩散方程的有限差分法matlab.pdf

. 一维扩散方程的有限差分法 ——计算物理实验作业七 陈万 物理学 2013 级 题目： 编程求解一维扩散方程的解 2 u u D (0 x a ,0 t t ) 2 0 max t x x u(x,t ) | e t 0 u a u b c (x 0) 1 1 1 n u a u b c (x a ) 2 2 2 0 n 取 a1 1,b1 1,c1 0, a2 1, b2 1,c2 0,a0 1.0,t max 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： 1 / 6 . 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); 2 / 6 . 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