MATLAB Answers

Error using pdepe - Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

13 views (last 30 days)
Jennifer Yang
Jennifer Yang on 25 Jul 2018
Answered: Bill Greene on 28 Jul 2018
Hello,
I'm using pdepe to solve for a reaction-diffusion equation. When I set the initial concentration to anything between 0-250, it runs fine, but the moments, I get this error saying
"Error using pdepe Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative."
Am I setting my boundary conditions incorrectly? I want the left hand side (at x = 0, C_L = C_L0) and on the right hand side I want the flux to equal to 0 (at x = x_f, dC_L/dx = 0).
function time_dependent_pdepe
clear all; close all; clc;
D_ij = 30;
L0 = 1000; %c0 [nM]
x_f =40; %Length of domain [um]
maxt = 2; %Max simulation time [s]
m = 0;
x = linspace(0,x_f,1000); %xmesh
t = linspace(0,maxt,50) %tspan
sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
u = sol;
% Plotting
hold all
for n = linspace(1,length(t),50)
plot(x,sol(n,:),'LineWidth',2)
end
title('Time Dependent (Diffusion Only)')
xlabel('Distance (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])
function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)
D = D_ij;
%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;
%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;
%Total Recptors
N_T = 1.7;
N_r = -(K_3.*K_4.*K_i.*(K_2 + u - ((8.*K_2.*u.^2.*N_T + 8.*K_2.*K_3.*u.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*u.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*u + 8.*K_2.*K_3.*K_i.*u.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*u.^2 + 4.*K_3.*u + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*u);
N_R = (u./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (u./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (u./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((u.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;
R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*u.*N_r))+...
((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*u.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
(d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*u.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));
% PDE
c = 1;
f = D_ij.*dudx;
s = R_L;
end
function u0 = DiffusionICfun(x)
u0 = 0;
end
function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)
c0 = L0;
% BCs: No flux boundary at the right boundary and constant concentration on
% the left boundary
% At x = 0, c0 = L0
pl = ul-c0;
ql = 0;
% At x = L, dc/dx = 0
pr = 0;
qr = 1;
end
end
Where am I going wrong?
Thank you.

Accepted Answer

Bill Greene
Bill Greene on 28 Jul 2018
The fundamental cause of your problem is that your initial conditions are inconsistent with your boundary condition at x=0. Mathematically this is referred to as an ill-posed problem and the PDE has no solution.
Sometimes, if the inconsistency is not too large, the numerical algorithm in pdepe can adjust the initial conditions to make them consistent with the boundary conditions. But, in your case, for large values of L0, it is not able to do this.
One simple way to fix this problem is to modify your initial condition function slightly:
function u0 = DiffusionICfun(x)
if(x==0)
u0 = L0;
else
u0 = 0;
end
end

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!