Solving Chromatography PDE with MOL and ode15s

Click For Summary

Discussion Overview

The discussion revolves around the application of the method of lines (MOL) to solve a system of partial differential equations (PDEs) related to chromatography. Participants explore the numerical implementation of the method, including the discretization of the equations, the formulation of boundary conditions, and the behavior of the resulting system of ordinary differential equations (ODEs).

Discussion Character

  • Technical explanation
  • Mathematical reasoning
  • Debate/contested

Main Points Raised

  • One participant describes the PDEs they are attempting to solve, including the initial and boundary conditions, and expresses concern about the stability of their numerical solution, which tends to explode.
  • The same participant provides details about their discretization approach, including the gradient matrices for the first and second derivatives, and questions the correctness of their method.
  • Another participant shares a code snippet that implements the method of lines, detailing the setup of parameters and the structure of the function used to solve the PDEs.
  • The code includes a function for calculating the concentration profile over time and position, but does not clarify how it addresses the stability issues mentioned earlier.
  • There is an implicit challenge regarding the effectiveness of the discretization method and the handling of boundary conditions, as the first participant's results are not as expected.

Areas of Agreement / Disagreement

Participants do not reach a consensus on the correctness of the discretization method or the stability of the numerical solution. Multiple competing views regarding the implementation and potential issues remain unresolved.

Contextual Notes

The discussion highlights potential limitations in the discretization approach, including assumptions made in the formulation of the gradient matrices and the handling of boundary conditions. The stability of the numerical solution is also a concern, but specific causes are not identified.

msanx2
Messages
13
Reaction score
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
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 resistance 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
 
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
 
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 resistance 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
 

Similar threads

MATLAB Unable to produce desired results using state-space model with observer for inverted pendulum
  • · Replies 0 ·
Replies
0
Views
2K
Mathematica Solving a complicated equation for approximate analytical Solution using Mathematica
  • · Replies 10 ·
Replies
10
Views
3K
Mathematica Solving equation and plotting issue with Mathematica
  • · Replies 2 ·
Replies
2
Views
2K
Undergrad Problem with calculating projections of curl using rotation of contour
  • · Replies 9 ·
Replies
9
Views
1K
MATLAB Solving transportation problems (LPP) using linprog in MATLAB
  • · Replies 3 ·
Replies
3
Views
5K
Undergrad Convention for writing $\binom{1}{1}$ tensors in matrix form?
  • · Replies 8 ·
Replies
8
Views
2K
Rotation of the stress tensor
  • · Replies 8 ·
Replies
8
Views
2K
Undergrad Stuck on solving a non homogenous (linear) PDE
  • · Replies 36 ·
2
Replies
36
Views
5K
LaTeX Need help -- Having trouble with eqnarray*/gathered/aligned
  • · Replies 1 ·
Replies
1
Views
2K
MATLAB Solving 2nd Order PDE System with Crank-Nicholson
  • · Replies 4 ·
Replies
4
Views
2K