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

 function[phi,potential]=pdex4EEpC(kVal,lVal)
 m = 0;

 x = linspace(0,4,100);
 t = linspace(0,2,500);

 nt=length(t);
 potential=zeros(nt,1);

 for i=1:nt
    potential(i)=appPot(t(i));

   end
  global K L tspan
  K=kVal;
  L=lVal;
  tspan=t;

options = odeset('RelTol',1e-9,'AbsTol',1e-9);

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t,options);
u1 = sol(:,:,1);
phi = zeros(1,nt);
for j = 1:nt
   % At time t(j), compute Du/Dt at x = 0.
   [~,phi(j)] = pdeval(m,t,u1(j,:),0);
end

figure
plot(potential,phi)
title('CV Plot')
xlabel('Potential')
ylabel('Current')

% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
global K L
c = [1; 1; 1]; 
f = [1; 1; 1] .* DuDx;
s = [u(2)-K*u(3)*u(1); -u(2)+K*u(1)*u(3); u(2)-K*u(1)*u(3)-L*u(3)]; 
% --------------------------------------------------------------
function u0 = pdex4ic(x)
u0 = [1; 0; 0]; 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
zeta = overPot(t);
conP = 1/(1+exp(zeta));
conQ = 1/(1+exp(-zeta));
pl = [ul(1)-conP;ul(2)-conQ;0];
ql = [0;0;1]; 
pr = [ur(1)-1; ur(2);0]; 
qr = [0; 0; 0];

% Defined as functions in the directory
function [e]=appPot(t)
t0=1;
Ei=-2;
Ef=2;
v=Ef-Ei;

if t<t0
    e=Ei+(v*t);
else
    e=Ef-(v*(t-t0));
end
end
% Defined as function in the directory
function [zeta]=overPot(t)

E0=1.2;
E=appPot(t);
zeta = (25.693*0.001)*(E-E0);

Answer 1

[Torsten]

You forgot to define the boundary condition of the third equation in xr.

Best wishes

Torsten.

[Kiran Vaddi]

Hi,I changed it to the following but I still get the same error

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
zeta = overPot(t);
conP = 1/(1+exp(zeta));
conQ = 1/(1+exp(-zeta));
pl = [ul(1)-conP;ul(2)-conQ;0];
ql = [0;0;1]; 
pr = [ur(1)-1; ur(2);ur(3)]; 
qr = [0; 0; 0];

I am trying to run the solution using the following

[phi,potential]=pdex4EEpC(1e6,1e-4)

[Torsten]

Use this function to supply the initial conditions:

function u0 = pdex4ic(x)

if x==0 
   zeta = overPot(0.0);
   conP = 1/(1+exp(zeta));
   conQ = 1/(1+exp(-zeta));
   u0 = [conP;conQ;0];
else
   u0 = [1; 0; 0];
end

end

[Kiran Vaddi]

Hi, Thanks for the help I was also trying to solve a fairly similar problem with just two variables this time (u1,u2). But When I try to solve it with the parameters

[phi,potential]=pdex4_EEp(1e-2,1e3);

I get the following error:

Warning: Failure at t=1.296555e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (3.552714e-15) at time t. 
> In ode15s (line 668)
  In pdepe (line 289)
  In pdex4_EEp (line 20) 
Warning: Time integration has failed. Solution is available at requested time points up to t=1.294589e+00. 
> In pdepe (line 303)
  In pdex4_EEp (line 20)

I would think my mesh is fairly okay but not sure. Any help in this regard will be appreciated

New solving code us here for two variables

function[phi,potential]=pdex4_EEp(kVal,lVal)
m = 0;

x = linspace(0,4,500);
t = linspace(0,2,500); %set t0 in the appPot.m accordingly

nt=length(t);
potential=zeros(nt,1);

for i=1:nt
    potential(i)=appPot(t(i));
end

global K L 
K=kVal;
L=lVal;

options = odeset('RelTol',1e-5,'AbsTol',1e-7);

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t,options);
u1 = sol(:,:,1);
phi = zeros(1,nt);
for j = 1:nt
   % At time t(j), compute Du/Dx at x = 0.
   [~,phi(j)] = pdeval(m,x,u1(j,:),0);
end

figure
plot(potential,phi)
title('CV Plot')
xlabel('Potential')
ylabel('Current')

% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
global K L
c = [1/L; 1/L]; 
f = [K; K] .* DuDx;
s = [-u(2); u(2)]; 
% --------------------------------------------------------------
function u0 = pdex4ic(x)
if x==0 
   zeta = overPot(0.0);
   conP = 1/(1+exp(zeta));
   conQ = 1/(1+exp(-zeta));
   u0 = [conP;conQ];
else
   u0 = [1; 0];
end

% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
zeta = overPot(t);
conP = 1/(1+exp(zeta));
conQ = 1/(1+exp(-zeta));
pl = [ul(1)-conP; ul(2)-conQ];
ql = [0; 0]; 
pr = [ur(1)-1; ur(2)]; 
qr = [0; 0];

[Bill Greene]

What are overPot and appPot?

[Torsten]

function [e]=appPot(t)
t0=1;
Ei=-2;
Ef=2;
v=Ef-Ei;

if t<t0
    e=Ei+(v*t);
else
    e=Ef-(v*(t-t0));
end
end
% Defined as function in the directory
function [zeta]=overPot(t)

E0=1.2;
E=appPot(t);
zeta = (25.693*0.001)*(E-E0);

[Bill Greene]

The solution is becoming very large as t increases. If you solve to .02 seconds and then plot the first solution variable, the problem will be obvious.

[Kiran Vaddi]

So what would be a better way to solve this system instead? Do you have any suggestions?

[Bill Greene]

When the solution goes to plus or minus infinity, usually the equations are mathematically ill-posed. Do you have a reference for the equations you are trying to solve?