temp.m

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

function out=PDE1Dd_CrankNicolson
% 使用Crank-Nicolson有限差分方法求解一维动态传热模型
%
% calling the function U = CrankNicolson(f,c1,c2,a,b,c,n,m) to solve the heat equation: 
%   du(x,t)/dt = α * 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
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/04/15 $
tic

% Increment of r and z
c1=8;
a =1;
b =5;
a1=10
b1=0.9533;
k1=0.3572;
n =150;
m =150;  
dx=a/(n-1);
dt=b/(m-1);
% input parameters
t_bi=12.5173;
Pe=30;
porosity=0.8;
F=(1-porosity)/porosity;
alpha=dt/(2*dx*dx*Pe);
beta=-dt/(4*dx);

     
U=CrankNicolson(@ic,c1,a,b,a1,b1,k1,alpha,beta,n,m,F);
toc
t=0:dt:b;
expe1=[10,20,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76;0,0.0015,0.003375,0.00725,0.021,0.06475,0.1285,0.231,0.3285,0.41225,0.506,0.596,0.68975,0.7335,0.8135,0.84225,0.86725,0.8885,0.921,0.93225,0.936,0.95725,0.946,0.96475,0.9685,1.001,0.961,0.9785,0.96475,0.97225,0.9685,0.9685,0.9935,0.9935,0.96475,0.9935,0.9785,0.99225,0.9935,0.9685,0.96475,0.9685,0.99975,0.98225,1.001,0.9685,0.9785,1.001,1.001;]';

plot(t'*t_bi,U(:,end),'',expe1(:,1),expe1(:,2),'mo');

% ------------------------------------------------------------------
function f = ic(x)
f = 0;