Chromatography PDE MOL

  • #1
13
0

Main Question or Discussion Point

Hello all

I am using the method of lines to solve the following PDE:

## \frac {\partial C} {\partial t} + F\frac {\partial q} {\partial t} + u \frac {dC} {dz} = D_{ax} \frac{\partial^2 C} {\partial z^2} ##

## \frac {\partial q} {\partial t} = k (q^{*}-q) ##

With these initial conditions:
## C(0,z) = 0 ##
## q(0,z) = 0 ##

and boundary conditions:
## C(t,0) = C_{in} (1) ##
## \frac {\partial C} {\partial z}(t,L) = 0 (2) ##

I discretized the equations above in space and obtained a system of ODE's which I solved with ode15s. To deal with the boundary condition (2) I said:

## \frac{C(t,L+1) - C(t,L-1)}{2\Delta z}= 0 ## or ## C(t,L+1) = C(t,L-1) ##

The problem is that my system gives a solution which tends to explode and makes no sense. My guess is that the way I am discretizing the system is not correct or consistent. I am using between 50-150 nodes.

My gradient matrix for the 1st derivative (GM1) looks like:

\begin{bmatrix}
0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & ... & 0\\
\frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & 0 & ... & 0\\
\frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & ... & 0 \\
0 & \frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & ... & 0 \\
... & ... & ... & ... & ... & ... & ... & ... & ... \\
0 & ... & 0 & 0 & 0 & \frac{1}{12} & \frac{-8}{12} & \frac{-1}{12} & \frac{8}{12} \\
0 & ... & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}

and b1 = \begin{bmatrix}
-\frac{C_{in}}{2} \\
\frac{C_{in}}{12} \\
0 \\
... \\
0 \\
-C_{in}
\end{bmatrix}

For the 2nd derivative (GM2):

\begin{bmatrix}
-2 & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0\\
1 & -2 & 1 & 0 & 0 & 0 & 0 & ... & 0\\
0 & 1 & -2 & 1 & 0 & 0 & 0 & ... & 0 \\
... & ... & ... & ... & ... & ... & ... & ... & ... \\
0 & ... & 0 & 0 & 0 & 0 & 0 & 2 & -2
\end{bmatrix}

and b2 = \begin{bmatrix}
C_{in} \\
0 \\
... \\
0
\end{bmatrix}

In the end I say that:

## \frac {\partial C} {\partial z} = \frac{1}{h}GM_{1} \cdot C + \frac{1}{h}b_{1}## and ## \frac {\partial^2 C} {\partial z^2} = \frac{1}{h^2}GM_{2} \cdot C + \frac{1}{h^2}b_{2}##
 
Last edited:

Answers and Replies

  • #2
13
0
My code:

FUNCTION 1:

Code:
function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
% Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
% using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
% Change the isotherm model by modifying the isotherm function
% Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
% Multicomponent
% Uses pdepe functin to solve the system of partial differential equations


global data Grad1 Grad2 res Cs_col

%% Default arguments (Example)
if nargin == 0
    isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
    feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
    parameter = [1e3 0.3/1e3 0.05; 1e3 0.1/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
    L =         10;                                                                 % cm, column length
    Di =        1;                                                                  % cm, column internal diameter
    epsb =      0.4;                                                                % column bulk porosity
    Q =         4;                                                                  % mL/min, flow rate
    Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
    KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
    Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
    tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
    tfinal =    5;                                                                  % min, final time for calculation
    opt.npz =   50;                                                                % number of discretization points in z
    opt.npt =   100;                                                                % number of discretization points in t
    opt.fig =   true;                                                               % true - show figures; false - do not show figures

    Cs_wash = 0; % Salt concentration in mM
    twash = 1;
  
    t0_elu = tpulse+twash; telu = 8; Cf = 80;
    Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

    Cs_reg = 80;
    treg = 3;
  
    Cs_eq = 0;
    teq = 3;
end

%% Calculations
addpath('../')
data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
nc = length(Cfeed);
data.nc = nc;
A = pi()*Di^2/4; % cm2
data.F = (1-epsb)/epsb;
data.ui = Q/A; % cm/min

% Run MOL
t = linspace(0,tfinal,opt.npt);
z = linspace(0,L,opt.npz);

[Grad1,Grad2] = GradientMatrix(opt.npz);
u0 = zeros(2*nc*opt.npz,1);

[t,uout] = ode15s(@odefun, t, u0);
for CompNo = 1:data.nc
    idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;
  
    C{CompNo} = uout(:,idx1:idx2);
end
res.t = t'; res.z = z'; res.C = C;

Cs_col = zeros(opt.npt,opt.npz);

Cs_col(:,1) = Cs_profile(t);
for i = 1:eek:pt.npt
    for j = 2:eek:pt.npz
        if t(i)-z(j)/data.ui >= 0
            Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
        else
            Cs_col(i,j) = 0;
        end
    end
end

ymat = [];
for CompNo = 1:data.nc
    ymat = [ymat C{CompNo}(:,end)];
end  
[ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
axis([0 tfinal  0 inf]) % fix the axes
set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
xlabel('t')
ylabel('C')

% Concentration over time with slider figure
    figure
    hold all
    slidermin = 1; % t = 0
    slidermax = opt.npt; % t = tfinal
    ymax_aux = [];
    for CompNo = 1:data.nc-1
        res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
        ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
    end
    ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
    res.ymax = str2double(num2str(max(ymax_aux),1));
    axis([0 L 0 2]);
  
    [ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
    set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
    hold off
    res.hslider(data.nc) = h1;
    res.hslider(data.nc+1) = h2;
  
    xlabel('Position, z')
    ylabel('Concentration, C')
    anno_text = sprintf('t = %.1f', 0);
    res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
    sld = uicontrol('Style', 'slider',...
        'Min',slidermin,'Max',slidermax,'Value',slidermin,...
        'SliderStep',[1 1]./(slidermax-slidermin),...
        'Position', [390 390 120 20],...
        'Callback', @sliderplot);

%% Auxiliary functions
function sliderplot(source,event)
% Handles the UI slider in the Concentration vs Position plot
global data res Cs_col

val = round(get(source,'Value'));
for CompNo = 1:data.nc
    set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
end
set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

anno_text = sprintf('t = %.2f',res.t(val));
set(res.time_box,'String',anno_text);

function udot = odefun(t,u)
    global data Grad1 Grad2
  
    nc = data.nc; npz = data.npz; L = data.L;
    z = linspace(0,L,npz); h = L/(npz-1);
  
    Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);
  
    C = u(1:nc*npz);
    for CompNo = 1:nc
        idx1 = 1+(CompNo-1)*npz;
        C(idx1) = Cin(CompNo);
    end
  
    q = u(nc*npz+1:end);
    Cs = zeros(npz,1);
  
    udot = zeros(2*nc*npz,1);  
    for j = 1:npz
        if t-z(j)/data.ui >= 0
            Cs(j) = Cs_profile(t-z(j)/data.ui);
        else
            Cs(j) = 0;
        end
    end
  
    qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

    dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
    for CompNo = 1:nc
        idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
        idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;
      
        b1 = zeros(npz,1); b1(2) = -Cin(CompNo)/2; b1(3) = Cin(CompNo)/12;
        b1(end) = -C(idx2);
        b2 = zeros(npz,1); b2(2) = Cin(CompNo);

        dCdz(idx1:idx2) = Grad1*C(idx1:idx2)./h + b1./h;
        d2Cdz2(idx1:idx2) = Grad2*C(idx1:idx2)./h^2 + b2./h^2;
      
        udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
        udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

        udot(idx1)=0; udot(idx3) = 0;
    end
  
%% Feed profile
function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
% Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

if strcmp(feedProf,'pulse')
    if t<=tpulse
        Cinj = Cfeed;
    else
        Cinj = zeros(1,length(Cfeed));
    end
  
elseif strcmp(feedProf,'step')
    Cinj = Cfeed;
  
else
    error('Invalid feed profile. feedProf must be "step" or "pulse"')
end

function Cs_in = Cs_profile(t)
global data

Cs_in = zeros(length(t),1);

for i = 1:length(t)
  
    if t(i) <= data.tpulse
        Cs_in(i) = 0;
    elseif t(i) <= data.tpulse+data.twash
        Cs_in(i) = data.Cs_wash;
    elseif t(i) <= data.tpulse+data.twash+data.telu
        Cs_in(i) = data.Cs_elu(t(i));
    elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
        Cs_in(i) = data.Cs_reg;
    else
        Cs_in(i) = data.Cs_eq;
    end
end

function [Grad1,Grad2] = GradientMatrix(N)

Grad1 = zeros(N,N);
Grad1(2,2) = 1/2; Grad1(3,1:4) = [-8/12, 0, 8/12, -1/12];
Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;

for i = 4:N-2
    Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
end

Grad2 = zeros(N,N); Grad2(2,1:2) = [-2, 1];
Grad2(end,end-1:end) = [2,-2];

for i = 3:N-1
    Grad2(i,i-2:i) = [1, -2, 1];
end
 
  • #3
13
0
FUNCTION 2

Code:
%% Isotherm Library

function q=isotherm1(isoType, C, Cs, nc, parameter)
% Defines the isotherm
parameter = transpose(parameter); N = length(Cs);

q = zeros(N,nc);
Cmat = zeros(N,nc);
for i = 1:nc
    idx1 = 1+(i-1)*N; idx2 = i*N;
    Cmat(:,i) = C(idx1:idx2);
end

if strcmp(isoType,'mod-langmuir')
% Linear-Langmuir isotherm:   q1(C1,C2) = a1*b1*exp(-m1*Cs)*C1/(1+b1*exp(-m1*Cs)*C1) , parameter = [a1 ; b1 ; m1]
    for i = 1:N
        langmuirDenominator = 1;
        for j = 1:nc
            langmuirDenominator = langmuirDenominator +  parameter(2,j)*exp(-parameter(3,j)*Cs(i))*Cmat(i,j);
        end
        q(i,:) = parameter(1,:).*parameter(2,:).*exp(-parameter(3,:)*Cs(i)).*Cmat(i,:)/langmuirDenominator;
    end
 
else
    error('Invalid isotherm type. isoType must be "linear" or "langmuir" or "linear-langmuir"')
end
 
  • #4
13
0
I think now I kind of solved it. However when the wave reaches the right boundary there is some noise in the wave. Do you know what could be the reason?

Code:
function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
% Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
% using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
% Change the isotherm model by modifying the isotherm function
% Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
% Multicomponent
% Uses pdepe functin to solve the system of partial differential equations


global data Grad1 Grad2 res Cs_col

%% Default arguments (Example)
if nargin == 0
    isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
    feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
    parameter = [1e3 1.5/1e3 0.15; 1e3 2/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
    L =         10;                                                                 % cm, column length
    Di =        1;                                                                  % cm, column internal diameter
    epsb =      0.4;                                                                % column bulk porosity
    Q =         4;                                                                  % mL/min, flow rate
    Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
    KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
    Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
    tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
    tfinal =    8;                                                                  % min, final time for calculation
    opt.npz =   300;                                                                % number of discretization points in z
    opt.npt =   100;                                                                % number of discretization points in t
    opt.fig =   true;                                                               % true - show figures; false - do not show figures

    Cs_wash = 0; % Salt concentration in mM
    twash = 1;
   
    t0_elu = tpulse+twash; telu = 8; Cf = 80;
    Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

    Cs_reg = 80;
    treg = 3;
   
    Cs_eq = 0;
    teq = 3;
end

%% Calculations
addpath('../')
data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
nc = length(Cfeed);
data.nc = nc;
A = pi()*Di^2/4; % cm2
data.F = (1-epsb)/epsb;
data.ui = Q/A; % cm/min

% Run MOL
t = linspace(0,tfinal,opt.npt);
z = linspace(0,L,opt.npz);

[Grad1,Grad2] = GradientMatrix(opt.npz);
u0 = zeros(2*nc*opt.npz,1);

[t,uout] = ode15s(@odefun, t, u0);
for CompNo = 1:data.nc
    idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;
   
    C{CompNo} = uout(:,idx1:idx2);
end
res.t = t'; res.z = z'; res.C = C;

Cs_col = zeros(opt.npt,opt.npz);

Cs_col(:,1) = Cs_profile(t);
for i = 1:opt.npt
    for j = 2:opt.npz
        if t(i)-z(j)/data.ui >= 0
            Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
        else
            Cs_col(i,j) = 0;
        end
    end
end

ymat = [];
for CompNo = 1:data.nc
    ymat = [ymat C{CompNo}(:,end)];
end   
[ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
axis([0 tfinal  0 inf]) % fix the axes
set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
xlabel('t')
ylabel('C')

% Concentration over time with slider figure
    figure
    hold all
    slidermin = 1; % t = 0
    slidermax = opt.npt; % t = tfinal
    ymax_aux = [];
    for CompNo = 1:data.nc-1
        res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
        ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
    end
    ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
    res.ymax = str2double(num2str(max(ymax_aux),1));
    axis([0 L 0 2]);
   
    [ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
    set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
    hold off
    res.hslider(data.nc) = h1;
    res.hslider(data.nc+1) = h2;
   
    xlabel('Position, z')
    ylabel('Concentration, C')
    anno_text = sprintf('t = %.1f', 0);
    res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
    sld = uicontrol('Style', 'slider',...
        'Min',slidermin,'Max',slidermax,'Value',slidermin,...
        'SliderStep',[1 1]./(slidermax-slidermin),...
        'Position', [390 390 120 20],...
        'Callback', @sliderplot);

%% Auxiliary functions
function sliderplot(source,event)
% Handles the UI slider in the Concentration vs Position plot
global data res Cs_col

val = round(get(source,'Value'));
for CompNo = 1:data.nc
    set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
end
set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

anno_text = sprintf('t = %.2f',res.t(val));
set(res.time_box,'String',anno_text);

function udot = odefun(t,u)
    global data Grad1 Grad2
   
    nc = data.nc; npz = data.npz; L = data.L;
    z = linspace(0,L,npz); h = L/(npz-1);
   
    Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);
   
    C = u(1:nc*npz);
    for CompNo = 1:nc
        idx1 = 1+(CompNo-1)*npz;
        C(idx1) = Cin(CompNo);
    end
   
    q = u(nc*npz+1:end);
    Cs = zeros(npz,1);
   
    udot = zeros(2*nc*npz,1);   
    for j = 1:npz
        if t-z(j)/data.ui >= 0
            Cs(j) = Cs_profile(t-z(j)/data.ui);
        else
            Cs(j) = 0;
        end
    end
   
    qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

    dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
    for CompNo = 1:nc
        idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
        idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;
       
        b1 = zeros(npz-1,1); b1(1) = -Cin(CompNo)/2; b1(2) = Cin(CompNo)/12;
        b1(end) = -C(idx2);
        b2 = zeros(npz-1,1); b2(1) = Cin(CompNo);
       
        dCdz(idx1) = 3.14; % Value does not matter
        d2Cdz2(idx1) = 3.14; % Value does not matter
       
        dCdz(idx1+1:idx2) = Grad1*C(idx1+1:idx2)./h + b1./h;
        d2Cdz2(idx1+1:idx2) = Grad2*C(idx1+1:idx2)./h^2 + b2./h^2;
       
        udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
        udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

        udot(idx1)=0; udot(idx3) = 0;
    end
   
%% Feed profile
function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
% Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

if strcmp(feedProf,'pulse')
    if t<=tpulse
        Cinj = Cfeed;
    else
        Cinj = zeros(1,length(Cfeed));
    end
   
elseif strcmp(feedProf,'step')
    Cinj = Cfeed;
   
else
    error('Invalid feed profile. feedProf must be "step" or "pulse"')
end

function Cs_in = Cs_profile(t)
global data

Cs_in = zeros(length(t),1);

for i = 1:length(t)
   
    if t(i) <= data.tpulse
        Cs_in(i) = 0;
    elseif t(i) <= data.tpulse+data.twash
        Cs_in(i) = data.Cs_wash;
    elseif t(i) <= data.tpulse+data.twash+data.telu
        Cs_in(i) = data.Cs_elu(t(i));
    elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
        Cs_in(i) = data.Cs_reg;
    else
        Cs_in(i) = data.Cs_eq;
    end
end

function [Grad1,Grad2] = GradientMatrix(N)

% Grad1 = zeros(N,N);
% Grad1(2,2) = 1/2; Grad1(3,1:4) = [-8/12, 0, 8/12, -1/12];
% Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;
%
% for i = 4:N-2
%     Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
% end
Grad1 = zeros(N-1,N-1);
Grad1(1,2) = 1/2; Grad1(2,1:4) = [-8/12, 0, 8/12, -1/12];
Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;

for i = 3:N-3
    Grad1(i,i-2:i+2) = [1/12, -8/12, 0, 8/12, -1/12];
end

Grad2 = zeros(N-1,N-1); Grad2(1,1:2) = [-2, 1];
Grad2(end,end-1:end) = [2,-2];

for i = 2:N-2
    Grad2(i,i-1:i+1) = [1, -2, 1];
end
 

Related Threads on Chromatography PDE MOL

Replies
2
Views
2K
  • Last Post
Replies
3
Views
566
Replies
4
Views
4K
Replies
4
Views
680
Replies
4
Views
3K
Replies
4
Views
1K
Replies
1
Views
895
Replies
3
Views
759
  • Last Post
Replies
9
Views
1K
Top