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);

% 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 nsn’, 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;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

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 nsn’, 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; 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

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 nsn’, 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; pde, pdepe MATLAB Answers — New Questions

Cart

Cart