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

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

——计算物理实验作业七

陈万

● 题目：

编程求解一维扩散方程的解

取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 ）

有限差分解如下：

解析解如下：

分析数据可知误差量级为)(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列对角元素直接拿出来，因此在调用时应当注意各对角元素的位置，避免调用错误。 （4）