cranknicolson.m

有限差分法（crank_nicholson finite difference ) 求解色谱模型

function U = CrankNicolson(f,c1,theta,fai,yita,keli,lamda,a,b,n,m)
% 使用Crank-Nicolson有限差分方法求解一维动态传热模型
% Crank-Nicolson Method for the heat equation: du(x,t)/dt = alpha * d2u(x,t)/dx2 over R={(x,t):0<=x<=a,
%   0<=t<=b} with u(x,0)=f(x), for 0<=x<=a, and u(0,t)=c1,u(a,t)=c2,for 0<=t<=b.
%
% Input     - f=u(x,0) as a string 'f'
%           - c1=u(0,t) and c2=u(a,t)
%           - a and b right end points of [0,a] and [0,b]
%           - c the constant in the heat equation
%           - n and m number of grid points over [0,a] and [0,b]
% Output    - U solution matrix

% Initialize parameters and U
h = a/(n-1);
k = b/(m-1);

U = zeros(n,m);
q=zeros(n,m);
% Boundary conditions
U(1,1:m) =c1;
a1=72.973;
b1=0.4840;
k1=0.3273;


% Generate first row (that is, the first column in matrix U)
U(1:n-1,1) = feval(f,h:h:(n-2)*h)';
q(1:n-1,1)=0;
% Form the diagonal and off-diagonal elements of A and
% the constant vector B and solve tridiagonal system AX=B
Vd(1:n) = ones(1,n);
Vd(n) = 1-theta;
Va = -theta*ones(1,n-1);
Vc = -fai*ones(1,n-1);
Vc(1)=0;

Vb(1)=c1;
g=ones(n,m);

for j=2:m
    for i=2:n-1
        q(i,j)=lamda*a1*U(i,j-1)/(1+b1*U(i,j-1))+lamda*k1*U(n,j-1)-(lamda-1)*q(i,j-1);
        g(i,j-1)=a1*U(i,j-1)/(1+b1*U(i,j-1))+k1*U(i,j-1);
        Vb(i) = theta*U(i-1,j-1)+fai*U(i+1,j-1)+yita*(g(i,j-1)-q(i,j-1))-keli*U(i,j-1);
    end
    g(n,j-1)=a1*U(n,j-1)/(1+b1*U(n,j-1))+k1*U(n,j-1);

    Vb(n)=theta*U(n-1,j-1)+fai*(2*U(n,j-1)-U(n-1,j-1))+yita*(g(n,j-1)-q(n,j-1))-keli*U(n,j-1);
    X = TDMA(Va,Vd,Vc,Vb);
    U(1:n,j)=X';
end
U = U';