Kinks/discontinuities in pdepe solution

Hello everyone,

I've been working with Matlab pdepe to solve numerically a set of coupled partial differential equations. I'm fairly new and I've encountered a problem with the numerical solution that I get from my code and can't seem to get rid of it. The equations model the diffusion and decay of two species in a material over 100s of ns. The problem I'm getting is that when I look at the solution in a log time scale, several kinks appear. Here are different things I've tried and haven't helped:

• Make tspan finer. (30-1000 points)
• Make xmesh finer. (50-1000 points)
• Increase or decrease relative tolerance (1E0 to 1E-6).
• Make tspan time separation linear after t = 0.
• Make tspan time separation increase logarithmically after t = 0.1 s

Any insights into how I could solve the problem would be useful.

% Physical constants
%params = [DR TR TT A B];
DR_guess = 1E6; % in nm^2 /ns
DT_guess = 1E1; % in nm^2 /ns
TR_guess = 14; % in ns 
TT_guess = 130; % in ns
T0 = 295; % in K
R0 = 0.016;% in e/nm^3
A  = 1;
B  = 200;
Y  = 0.3056;
CR = -80000;
CT = 220;

params = [DR_guess DT_guess TR_guess TT_guess T0 R0 A B Y CR CT];

% Setting up the mesh for space and time
xmesh = linspace(0,2000,55); %in nm
tspan = [-5, -2, -1, -0.5, 0, logspace(-1, 3, 80)];% in ns


% Solving the differential equation
options = odeset('RelTol',1e-0,'OutputFcn',@myOut,'Stats','on');
pde = @(x, t, u, dudx) AgSepde(x, t, u, dudx, params);
soln = pdepe(1, pde, @AgSeic, @AgSebc, xmesh, tspan, options);

15 successful steps
0 failed attempts
37 function evaluations
1 partial derivatives
5 LU decompositions
30 solutions of linear systems

% Partial  differential equation function
function [c,f,s] = AgSepde(x,t,u,dudx,params)  
  % Parameters for PDEs
  DR = params(1);
  DT = params(2);
  TR = params(3);
  TT = params(4);
  T0 = params(5);
  R0 = params(6);
  A  = params(7);
  B  = params(8);
  Y  = params(9);

  
  % Generation function
  sigmar = 180; % in nm
  sigmat = 0.2; %in ns
  t0 = 0.00; % in ns
  G = (1 / (sigmar * sqrt(2*pi))) * exp(-x.^2 / (2 * sigmar^2)) ...
    * (1 / (sigmat * sqrt(2*pi))) * exp(-(t - t0).^2 / (2 * sigmat^2));

  c = [1; 1];
  f = [DR*dudx(1); DT*dudx(2)];
  s = [-(u(1) - R0)/TR  + A*G; Y*(u(1)- R0) - (u(2) - T0)/TT  + B*G];
end

% Intial conditions for the PDEs
function u0 = AgSeic(x)
    T0 = 295;
    R0 = 0.016;
    u0 = [R0; T0];
end


function [pl, ql, pr, qr] = AgSebc(xl, ul, xr, ur, t)
    % Left boundary: u(0,t) = 1
    pl = [0; 0];   
    ql = [1; 1];   

    % Right boundary: u(x',t) = 0
    pr = [0; 0];
    qr = [1; 1];
end


function status = myOut(t,y,flag)
    persistent last
    if strcmp(flag,'init'), last = tic; end
    if isempty(flag) && toc(last) > 1
        fprintf('t = %.3g ns\n', t(end));
        last = tic;
    end
    status = 0;
end

%%  Figures
tiledlayout(2,1)
nexttile(1,[1,1])
hold on;
plot(tspan,soln(:,1,1),'Marker','o');
set(gca,'XScale','log');
xlabel('Time / ns');
ylabel('Species 1');
hold off;

nexttile(2,[1,1])
hold on;
plot(tspan,soln(:,1,2),'Marker','square')
set(gca,'XScale','log');
xlabel('Time / ns');
ylabel('Species 2');
hold off;

Warning: Negative data ignored

Warning: Negative data ignored

[Walter Roberson]

You should expect kinks near 0 since you are using negative time.

% Physical constants
%params = [DR TR TT A B];
DR_guess = 1E6; % in nm^2 /ns
DT_guess = 1E1; % in nm^2 /ns
TR_guess = 14; % in ns 
TT_guess = 130; % in ns
T0 = 295; % in K
R0 = 0.016;% in e/nm^3
A  = 1;
B  = 200;
Y  = 0.3056;
CR = -80000;
CT = 220;

params = [DR_guess DT_guess TR_guess TT_guess T0 R0 A B Y CR CT];

% Setting up the mesh for space and time
xmesh = linspace(0,2000,55); %in nm
tspan = [-5, -2, -1, -0.5, 0, logspace(-1, 3, 80)];% in ns


% Solving the differential equation
options = odeset('RelTol',1e-0,'OutputFcn',@myOut,'Stats','on');
pde = @(x, t, u, dudx) AgSepde(x, t, u, dudx, params);
soln = pdepe(1, pde, @AgSeic, @AgSebc, xmesh, tspan, options);

15 successful steps
0 failed attempts
37 function evaluations
1 partial derivatives
5 LU decompositions
30 solutions of linear systems

% Partial  differential equation function
function [c,f,s] = AgSepde(x,t,u,dudx,params)  
  % Parameters for PDEs
  DR = params(1);
  DT = params(2);
  TR = params(3);
  TT = params(4);
  T0 = params(5);
  R0 = params(6);
  A  = params(7);
  B  = params(8);
  Y  = params(9);

  
  % Generation function
  sigmar = 180; % in nm
  sigmat = 0.2; %in ns
  t0 = 0.00; % in ns
  G = (1 / (sigmar * sqrt(2*pi))) * exp(-x.^2 / (2 * sigmar^2)) ...
    * (1 / (sigmat * sqrt(2*pi))) * exp(-(t - t0).^2 / (2 * sigmat^2));

  c = [1; 1];
  f = [DR*dudx(1); DT*dudx(2)];
  s = [-(u(1) - R0)/TR  + A*G; Y*(u(1)- R0) - (u(2) - T0)/TT  + B*G];
end

% Intial conditions for the PDEs
function u0 = AgSeic(x)
    T0 = 295;
    R0 = 0.016;
    u0 = [R0; T0];
end


function [pl, ql, pr, qr] = AgSebc(xl, ul, xr, ur, t)
    % Left boundary: u(0,t) = 1
    pl = [0; 0];   
    ql = [1; 1];   

    % Right boundary: u(x',t) = 0
    pr = [0; 0];
    qr = [1; 1];
end


function status = myOut(t,y,flag)
    persistent last
    if strcmp(flag,'init'), last = tic; end
    if isempty(flag) && toc(last) > 1
        fprintf('t = %.3g ns\n', t(end));
        last = tic;
    end
    status = 0;
end

%%  Figures
tiledlayout(2,1)
nexttile(1,[1,1])
hold on;
plot(tspan,soln(:,1,1),'Marker','o');
%set(gca,'XScale','log');
xlim([-10 10])
xlabel('Time / ns');
ylabel('Species 1');
hold off;

nexttile(2,[1,1])
hold on;
plot(tspan,soln(:,1,2),'Marker','square')
%set(gca,'XScale','log');
xlim([-10 10])
xlabel('Time / ns');
ylabel('Species 2');
hold off;

Answer 1

[Torsten]

Are you sure you want to solve your problem in cylindrical coordinates (since you set m = 1 in the call to "pdepe") ?

Tighten your tolerances as done in the below code:

% Physical constants
%params = [DR TR TT A B];
DR_guess = 1E6; % in nm^2 /ns
DT_guess = 1E1; % in nm^2 /ns
TR_guess = 14; % in ns 
TT_guess = 130; % in ns
T0 = 295; % in K
R0 = 0.016;% in e/nm^3
A  = 1;
B  = 200;
Y  = 0.3056;
CR = -80000;
CT = 220;

params = [DR_guess DT_guess TR_guess TT_guess T0 R0 A B Y CR CT];

% Setting up the mesh for space and time
xmesh = linspace(0,2000,55); %in nm
tspan = [-5, -2, -1, -0.5, 0, logspace(-1, 3, 80)];% in ns

% Solving the differential equation
%options = odeset('RelTol',1e-0,'OutputFcn',@myOut,'Stats','on');
options = odeset('RelTol',1e-12,'AbsTol',1e-12,'OutputFcn',@myOut,'Stats','on');
pde = @(x, t, u, dudx) AgSepde(x, t, u, dudx, params);
soln = pdepe(1, pde, @AgSeic, @AgSebc, xmesh, tspan, options);

566 successful steps
15 failed attempts
1163 function evaluations
1 partial derivatives
84 LU decompositions
1156 solutions of linear systems

% Partial  differential equation function
function [c,f,s] = AgSepde(x,t,u,dudx,params)  
  % Parameters for PDEs
  DR = params(1);
  DT = params(2);
  TR = params(3);
  TT = params(4);
  T0 = params(5);
  R0 = params(6);
  A  = params(7);
  B  = params(8);
  Y  = params(9);

  
  % Generation function
  sigmar = 180; % in nm
  sigmat = 0.2; %in ns
  t0 = 0.00; % in ns
  G = (1 / (sigmar * sqrt(2*pi))) * exp(-x.^2 / (2 * sigmar^2)) ...
    * (1 / (sigmat * sqrt(2*pi))) * exp(-(t - t0).^2 / (2 * sigmat^2));

  c = [1; 1];
  f = [DR*dudx(1); DT*dudx(2)];
  s = [-(u(1) - R0)/TR  + A*G; Y*(u(1)- R0) - (u(2) - T0)/TT  + B*G];
end

% Intial conditions for the PDEs
function u0 = AgSeic(x)
    T0 = 295;
    R0 = 0.016;
    u0 = [R0; T0];
end


function [pl, ql, pr, qr] = AgSebc(xl, ul, xr, ur, t)
    % Left boundary: u(0,t) = 1
    pl = [0; 0];   
    ql = [1; 1];   

    % Right boundary: u(x',t) = 0
    pr = [0; 0];
    qr = [1; 1];
end


function status = myOut(t,y,flag)
    persistent last
    if strcmp(flag,'init'), last = tic; end
    if isempty(flag) && toc(last) > 1
        fprintf('t = %.3g ns\n', t(end));
        last = tic;
    end
    status = 0;
end

%%  Figures
tiledlayout(2,1)
nexttile(1,[1,1])
hold on;
plot(tspan,soln(:,1,1),'Marker','o');
set(gca,'XScale','log');
xlabel('Time / ns');
ylabel('Species 1');
hold off;

nexttile(2,[1,1])
hold on;
plot(tspan,soln(:,1,2),'Marker','square')
set(gca,'XScale','log');
xlabel('Time / ns');
ylabel('Species 2');
hold off;

Warning: Negative data ignored

Warning: Negative data ignored

[Jose]

Hi Torsten thanks for the help. I didn't realize I had to tighten the tolerances that far.

As per your question. I think I have to use cylindrical coordinates as the system I'm trying to model is set in cylindrical coordinates (i.e., 2D Gaussian, azimuthally isotropic).