Defining mass matrix for solving transport equation

1 次查看（过去 30 天）

显示更早的评论

Deepa Maheshvare 2019-1-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/438547-defining-mass-matrix-for-solving-transport-equation

评论： Deepa Maheshvare 2019-9-23

在 MATLAB Online 中打开

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?

13 个评论
显示 11更早的评论隐藏 11更早的评论

Deepa Maheshvare 2019-1-8

在 MATLAB Online 中打开

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 2019-1-8

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 2019-9-23

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?

请先登录，再进行评论。

请先登录，再回答此问题。