Solving Chromatography PDE with MOL and ode15s

In summary: Q = 1e-10; Cfeed = 1e-5; KLDF = 1e-4; Dax = 1e-6; tpulse = 1e-4; tfinal = 1e-8%----------------------------------------------------------------------% Solve the system of partial differential equations%----------------------------------------------------------------------% Set up the boundary conditions%----------------------------------------------------------------------% Solve the system of differential equations%----------------------------------------------------------------------% Print the resultsIn summary,The TDM method was used to solve the PDE for C. The isotherm type was changed from linear to linear-langmuir, and the
  • #1
msanx2
13
0
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:
Physics news on Phys.org
  • #2
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 equationsglobal 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
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
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 equationsglobal 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
 

1. What is chromatography PDE?

Chromatography PDE (partial differential equation) refers to a mathematical model used to describe the behavior of chromatography systems, which are commonly used in chemistry and biochemistry to separate and analyze mixtures of substances.

2. What is MOL in relation to chromatography PDE?

MOL (method of lines) is a numerical method used to solve differential equations, including chromatography PDEs. It involves discretizing the PDE into a system of ODEs (ordinary differential equations), which can then be solved using standard numerical techniques.

3. What is ode15s and how does it relate to solving chromatography PDEs?

ode15s is a numerical solver specifically designed for stiff systems of ODEs, which often occur when solving chromatography PDEs. It uses a variable-order, variable-step method to efficiently and accurately solve the equations. It is frequently used in conjunction with MOL to solve chromatography PDEs.

4. What are the benefits of using MOL and ode15s to solve chromatography PDEs?

Using MOL and ode15s allows for efficient and accurate solutions to chromatography PDEs, which can be difficult to solve analytically. The combination of these methods also allows for a more flexible and customizable approach to solving PDEs, as different numerical techniques can be applied to different parts of the problem.

5. Are there any limitations to using MOL and ode15s for chromatography PDEs?

While MOL and ode15s are powerful tools for solving chromatography PDEs, they do have limitations. For example, they may not be suitable for very complex or large systems, and the accuracy of the solutions can be affected by the specific parameters and assumptions used in the model. It is always important to carefully validate and verify the results obtained from numerical methods.

Similar threads

MATLAB Sliding mode observer in Matlab or Simulink
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
LaTeX Need help -- Having trouble with eqnarray*/gathered/aligned
    • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
Struggling With Plasma Physics Question
    • Advanced Physics Homework Help
Replies
7
Views
241
Majorana Fermions: Lagrangean and equations of motion
    • Advanced Physics Homework Help
Replies
1
Views
876
MATLAB Solving transportation problems (LPP) using linprog in MATLAB
    • MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
3K
Solving Klein Gordon’s equation
    • Calculus and Beyond Homework Help
Replies
1
Views
444
Perform suitable gauge transformations
    • Advanced Physics Homework Help
Replies
3
Views
393
A Newmann-Penrose Spin coefficients for Schwarschild metric
    • Differential Geometry
Replies
1
Views
987
A Gamma matrices and Gell-Mann's I - Y categorization
    • Quantum Physics
Replies
22
Views
1K
A Solve PDE w/ Comsol 5.3: Numerical Solution & Time Evolution
    • Differential Equations
Replies
5
Views
1K

Hot Threads

Recent Insights

Back
Top