Defining mass matrix for solving transport equation

I want to solve the following 1-D pde which is discretized in spacial direction to solve using ode15s.

Dc/Dt = D*D^c/Dx^2 - v*Dc/Dx.

With initial condition,

c(x,0) = Co

Boundary conditions ,

Dirichlet boundary condition at left boundary,

c(0,t) = cl;

Neumann boundary condition at right boundary,

dc/dx = 0;

I'm not very clear about defining the mass matrix for the input options in the call to ode15s

options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@fun,tspan,y0,options);

I'm trying to solve the above system as a index-1 DAE. The first row and second row of the mass matrix will contain one's in the diagonal entries.

I am not sure how the left boundary condition has to be specified in the mass matrix.

Could someone explain?

[Torsten]

Why don't you use "pdepe" ?

[Deepa Maheshvare]

I have been able to successfully use "pdepe" for solving the above equation.I understand pdepe discretizes the pde in space to convert to an ode and the resulting equation is solved using ode15s.

My actual system has other odes's. So I want to understand how to solve the above system as a DAE.This will help me in implementaing the same procedure for my actual system.

Many thanks.

[Torsten]

For your PDE from above, the mass matrix will be the identity matrix and your system will be a pure system of ODEs. Why do you expect something different ?

[Deepa Maheshvare]

Excuse me for the naive questions,

I'm confused about the boundary conditions.

If the space is discretized into 'N' nodes, there will be N ordinary differential equations.

The left boundary condition cl(algebraic variable) = (say) 100 , can be differentiated once to form an ode which gives cl' = 0

Similarly the initial condition co(algebraic variable) = (say) 10 will form an ode co' = 0

I'm confused about the left boundary condition.What happens to this? Is this directly treated as an ode.

I'm also not really clear about the size of the mass matrix.

Could you please clarify?

[Torsten]

At the right boundary (x=x_n)

D^2c/dx^2 ~

(dc_n/dx - dc_(n-1/2)/dx) /(x_n-x_(n-1/2)) ~

(dc_n/dx - (c_n-c_(n-1))/(x_n-x_(n-1)))/(x_n-x_(n-1/2))

Thus with your boundary condition dc_n/dx = 0 you get the ODE

dc_n/dt = D*(- (c_n-c_(n-1))/(x_n-x_(n-1)))/(x_n-x_(n-1/2))

with x_(n-1/2) = (x_n+x_(n-1))/2

[Deepa Maheshvare]

Thanks a lot for the clarification.

M = diag([zeros(2,1); ones(N-2,1);ones(1,1)]);

2 zeros corresponding to the 2 algebraic equations (one for initial and other for left boundary). N-2 1's for N-2 odes at the interior nodes and 1 ode pertaining to the node at the right boundary from what was explained in the above comment.

Is this above mass matrix right?

I would like to know if there is a way to define the order of equations and variables.

[Deepa Maheshvare]

Following a few examples , I have created the following functions for the above-mentioned system

clear all
global N x D
N = 11;
D = 900;
% Discretizing space
x = linspace(0,62,11);
% Initial condition
y0 = 4000*ones(N,1)
options = odeset('Mass','M','Jpattern','on');
[t,y] = ode15s('executeode',[0 0.5 1],y0,options);

function out = executeode(t,y,flag)
global N x D
if nargin < 3 | isempty(flag)
  out = fun(t,y,N);
else
  switch(flag)           
  case 'mass'
    out = diag([zeros(2,1); ones(N,1)]);
  case 'jpattern'
    out = Jacpat(N); % to be created
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end
%--------------------------------------------------------------
function yp = fun(t,y,P,N)
%  Evaluate the differential equations.
c = y;
f = zeros(N,1);
%  The are N equally spaced mesh points
h = 6.2;
%  Equations at interior mesh points:
for i = 2:(N-1)
  f(i) = D*(c(i+1) - 2*c(i) + c(i-1))/h^2 ;              
end
%  Boundary conditions at x = 0:
f(1) = 4100;
%  Boundary conditions at x end:  
f(N) = -D*((c(N)-c(N-1))/h)/(x(N)- (x(N)+x(N-1))/2);
yp = f;
%--------------------------------------------------------------

%% Jacobian /JPattern

Could you please advise on whether Jacbian matrix or the JPattern matrix has to be given as input when ode15s is called?

How should the JPattern or Jacobian be created?

[Torsten]

As I already wrote: A mass matrix does not need to be defined. At the left boundary, you can solve dc(1)/dt = 0 since c=cI for all times. At the right boundary, I provided an ODE for dc(n)/dt.

Before you care about providing a Jacobian pattern (which might be helpful to accelarate computation), you should try to make your code work in dense mode.

Best wishes

Torsten.

[Deepa Maheshvare]

Thanks a lot for the support.The following code works in dense mode. The results are in accordance with what is obtained using pdepe solver. I would like to know how the Jacobian pattern can be implemented.

clear all
global N x D
N = 11;
D = 900;
tspan = 0:0.1:20
% Discretizing space
x = linspace(0,62,11);
% Initial condition
y0 = [3500;4000*ones(N-1,1)];
%options = odeset('Mass','M');%,'Jpattern','on');
[t,yDAE] = ode15s('fun',tspan,y0)

function dc = fun(t,y)
global N x D
%  Evaluate the differential equations.
c  = y;
dc = zeros(N,1);
%  The are N equally spaced mesh points
h  = 6.2;
%  Equations at interior mesh points:
for i = 2:(N-1)
  dc(i) = D*(c(i+1) - 2*c(i) + c(i-1))/h^2 ;              
end
%  Boundary conditions at x = 0:
dc(1) = 0;
%  Boundary conditions at x end:  
dc(N) = -D*((c(N)-c(N-1))/h)/(x(N)- (x(N)+x(N-1))/2);
end

[Torsten]

Sorry, I have not yet worked with the option of "Jacobian Pattern". Thus I had to start at the same knowledge level as you do.

[Deepa Maheshvare]

No problem. Thanks a lot for your support.

There is an example available here.

In the file named 'ex12ode(t,y,flag)' in the examples provided in the above link,

the Jacobian Pattern matrix and the mas matrix are defined.From what I understand , the diagonal entries corresponding to the algebraic equations are zero .I suppose, for my case, since the algebraic equation of left boundary is written as dc(1)/dt = 0 there is no algebraic equation. Thus the mass matrix is identity matrix, as you already mentioned. (Hope it is right to understand that the mass matrix doesn't include the algebraic equation corresponding to the initial condition.)

However,I couldn't understand how Jacpat function is defined in 'ex12ode(t,y,flag)'.

Hope someone could help me with this.

Any help will be much appreciated.

[Torsten]

function S = Jacpat(N)
%  Define the sparsity structure of the Jacobian by indicating which 
%  variables appear in which equations.  This is done by inspection of
%  the function 'f' above for evaluating the equations.

S = sparse(N,N);
for i = 2:(N-1)
  S(i,i-1) = 1;      %  c(i-1)
  S(i,i) = 1;        %  c(i)
  S(i,i+1) = 1;      %  c(i+1) 
end
S(N,N-1) = 1;        %  c(N-1)  
S(N,N) = 1;          %  c(N)

[Deepa Maheshvare]

Hello, Could you please explain why half-step size is considered for discretizaition at the boundary node that was explained here . Is there any name for this type of scheme?