简单体会数值求解一维扩散方程

借用这个题目
然后用比较简单的matlab代码来求解

%空间0到1，切了5份，时间0到1，切了十份。
clc,clear
v = 1;
dx = 0.2;
x = 0:dx:1;
dt = 0.1;
nx = 6;
nt = 11;
t = 0:dt:1;
b = v/(2*dx*dx);
c = b;
a = 1/dt+b+c;

Uold = exp(x);
Unew = ones(1,nx);    
Unew(1) = exp(t(2));
Unew(nx) = exp(1+t(2));

%其实最关键的就是得到A E F这些矩阵，本案例用的是crank-nicolson差分格式，不同格式得到的A等矩阵不同，在openfoam中有特定的设定方法
A = [a -b 0 0;-c a -b 0;0 -c a -b;0 0 -c a];
E = [c (1/dt-b-c) b 0 0 0;0 c (1/dt-b-c) b 0 0;0 0 c (1/dt-b-c) b 0;0 0 0 c (1/dt-b-c) b];
F = [c 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 b];

%时间推进，大规模矩阵的话就不能简单地用A\D来求了，需要其他迭代法，比如gauss-seidel之类的，还有更高效的算法，注意一个点就是对于第一类边界条件的情况，U1和Unx其实不用解矩阵求解，边界条件就已经设定了，第二第三类边界条件暂时还不清楚，后续再补充。 
for k = 2:nt
    D = E*Uold' + F*Unew';
    Uin = A\D;
    Unew(2:nx-1) = Uin';
    Unew(1) = exp(t(k));
    Unew(nx) = exp(1+t(k));
    Uold = Unew;
end

plot(x,Unew)