MATLAB Solving Chromatography PDE with MOL and ode15s

AI Thread Summary
The discussion focuses on using the method of lines to solve a system of partial differential equations (PDEs) related to chromatography. The author encounters issues with the solution becoming unstable, potentially due to incorrect discretization or boundary condition handling. They provide details on their gradient matrices for first and second derivatives and the initial and boundary conditions applied. Despite attempts to refine the model, noise appears in the solution as the wave reaches the boundary. The conversation highlights the challenges of numerical stability in solving PDEs in chromatography applications.
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 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
 
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 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
 
Back
Top