MATLAB: Error using pdepe (line 293) Spatial discretization has failed.

MATLABpdepe

I keep getting following error for my code. not sure where I went wrong.
Error using pdepe (line 293) Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.
Error in Kronstadt_HW4_problem3_pdepe (line 28) sol=pdepe(m,pdepb,@pbIC,pbBC,x,t,options); %solution
function Kronstadt_HW4_problem3_pdepe
 %This is a dimensional formulation 
  close all
Nx=201;    %number of points in x direction
Nt=101;    %number of points in time
m=0;       %cartesian
D=2.4E-9;  %m^2/s, diffusivity
Na=0.025;  %mol/m^3-s, sink
Ci=0.04;   %mol/m^3, initial condition
Cw = 0.04;  %mol/m^3, this is the set concentration at the right boundary
Xmax=.1;   %m
tmax=1000000;  %sec
dx=Xmax/(Nx-1);
dt=tmax/(Nt-1);
t=linspace(0,tmax,Nt);
x=linspace(0,Xmax,Nx);
options=odeset; %('MaxStep',100);
 %anonymous functions to pass parameters
  pbBC=@(yl,Cl,yr,Cr,t)pbBC1(yl,Cl,yr,Cr,t,Cw);  
  pdepb=@(x,t,C,DCDx)pdepb1(x,t,C,DCDx,Na,D);
  sol=pdepe(m,pdepb,@pbIC,pbBC,x,t,options);  %solution      
  C=sol(:,:,1);                               %extract solution 
x_value1=6;
x_value2=16;   
x_value3=36; 
x_value4=51;
plot(t,C(:,x_value1),t,C(:,x_value2),t,C(:,x_value3),t,C(:,x_value4))   %Plot vs t at specified x locations
xlabel('Time (s)')
ylabel('Concentration (mol/m^3)')
legend([num2str(x_value1*dx-dx) ' m'],[num2str(x_value2*dx-dx) ' m'],[num2str(x_value3*dx-dx) ' m'],[num2str(x_value4*dx-dx) ' m'])
figure(2)   %Plot vs x at specified times
t1=2;
t2=11;
t3=31;
t4=101;
plot(x,C(t1,:),x,C(t2,:),x,C(t3,:),x,C(t4,:))
xlabel('x (m)')
ylabel('Concentration (mol/m^3)')
ylim([0 Cw])
xlim([0 Xmax])
legend([num2str(t1*dt-dt) ' sec'],[num2str(t2*dt-dt) ' sec'],[num2str(t3*dt-dt) ' sec'], [num2str(t4*dt-dt) ' sec'])
figure(3)   %Flux at both wall surfaces
flux_l=-D*(C(:,2)-C(:,1))/dx;    %mol/s-m

flux_r=-D*(C(:,Nx)-C(:,Nx-1))/dx;    %mol/s-m
plot(t,flux_l,t,flux_r)
legend('Flux (left wall)','Flux (right wall)')
xlabel('Time (sec)')
ylabel('Mass Flux (mol/s-m^2)')
ylim([0 .5E-8])
figure(4)   %Flux at both left surface
semilogy(t,flux_l)
legend('Flux (left wall)')
xlabel('Time (sec)')
ylabel('Mass Flux (mol/s-m^2)')
grid on
figure(5)
Ns=5;
ta=linspace(0,tmax,(Nt-1)/Ns);
xa=linspace(0,Xmax,(Nx-1)/Ns);
[X,T]=meshgrid(xa,ta);
for i=1:length(ta)
    for j=1:length(xa)
        C1(i,j)=C(1+i*Ns,1+j*Ns);
    end
end
surface(X,T,C1)
ylabel('Time (s)')
xlabel('Position (m)')
zlabel('Concentration (mol/m^3)')
function [c,f,s]=pdepb1(x,t,C,DCDx,Na,D)  %define pde system – see Matlab help for pdepe
c=1;           %coefficient on time derivative
f=D*DCDx;      %general flux   mol/s-m^3 (negative of diffusive flux)
s=Na;         %sink   mol/s-m^3
function C0=pbIC(x)   %initial condition
C0=0.04;   %mol/m^3
function [pl,ql,pr,qr]=pbBC1(yl,Cl,yr,Cr,t,Cw)  %boundary conditions
pl=0;   %left side bc. This one sets flux to zero since the left side is impermeable
ql=1;       %Form is pl + ql*f = 0 where f is flux from definition of pde system in function pdepb
pr=Cl-Cw;       %right side bc. This one sets concentration to Cw
qr=0;

Best Answer

Try
    function Kronstadt_HW4_problem3_pdepe
     %This is a dimensional formulation 
      close all
    Nx=201;    %number of points in x direction
    Nt=101;    %number of points in time
    m=0;       %cartesian
    D=2.4E-9;  %m^2/s, diffusivity
    Na=0.025;  %mol/m^3-s, sink
    Ci=0.04;   %mol/m^3, initial condition
    Cw = 0.04;  %mol/m^3, this is the set concentration at the right boundary
    Xmax=.1;   %m
    tmax=1000000;  %sec
    dx=Xmax/(Nx-1);
    dt=tmax/(Nt-1);
    t=linspace(0,tmax,Nt);
    x=linspace(0,Xmax,Nx);
    options=odeset; %('MaxStep',100);
     %anonymous functions to pass parameters
      pbBC=@(yl,Cl,yr,Cr,t)pbBC1(yl,Cl,yr,Cr,t,Cw);
      pbIC = @(x)pbIC1(x,Xmax,Cw);  
      pdepb=@(x,t,C,DCDx)pdepb1(x,t,C,DCDx,Na,D);
      sol=pdepe(m,pdepb,pbIC,pbBC,x,t,options);  %solution      
      C=sol(:,:,1);                               %extract solution 
    x_value1=6;
    x_value2=16;   
    x_value3=36; 
    x_value4=51;
    plot(t,C(:,x_value1),t,C(:,x_value2),t,C(:,x_value3),t,C(:,x_value4))   %Plot vs t at specified x locations
    xlabel('Time (s)')
    ylabel('Concentration (mol/m^3)')
    legend([num2str(x_value1*dx-dx) ' m'],[num2str(x_value2*dx-dx) ' m'],[num2str(x_value3*dx-dx) ' m'],[num2str(x_value4*dx-dx) ' m'])
    figure(2)   %Plot vs x at specified times
    t1=2;
    t2=11;
    t3=31;
    t4=101;
    plot(x,C(t1,:),x,C(t2,:),x,C(t3,:),x,C(t4,:))
    xlabel('x (m)')
    ylabel('Concentration (mol/m^3)')
    ylim([0 Cw])
    xlim([0 Xmax])
    legend([num2str(t1*dt-dt) ' sec'],[num2str(t2*dt-dt) ' sec'],[num2str(t3*dt-dt) ' sec'], [num2str(t4*dt-dt) ' sec'])
    figure(3)   %Flux at both wall surfaces
    flux_l=-D*(C(:,2)-C(:,1))/dx;    %mol/s-m

    flux_r=-D*(C(:,Nx)-C(:,Nx-1))/dx;    %mol/s-m
    plot(t,flux_l,t,flux_r)
    legend('Flux (left wall)','Flux (right wall)')
    xlabel('Time (sec)')
    ylabel('Mass Flux (mol/s-m^2)')
    ylim([0 .5E-8])
    figure(4)   %Flux at both left surface
    semilogy(t,flux_l)
    legend('Flux (left wall)')
    xlabel('Time (sec)')
    ylabel('Mass Flux (mol/s-m^2)')
    grid on
    figure(5)
    Ns=5;
    ta=linspace(0,tmax,(Nt-1)/Ns);
    xa=linspace(0,Xmax,(Nx-1)/Ns);
    [X,T]=meshgrid(xa,ta);
    for i=1:length(ta)
        for j=1:length(xa)
            C1(i,j)=C(1+i*Ns,1+j*Ns);
        end
    end
    surface(X,T,C1)
    ylabel('Time (s)')
    xlabel('Position (m)')
    zlabel('Concentration (mol/m^3)')
    function [c,f,s]=pdepb1(x,t,C,DCDx,Na,D)  %define pde system – see Matlab help for pdepe
    c=1;           %coefficient on time derivative
    f=D*DCDx;      %general flux   mol/s-m^3 (negative of diffusive flux)
    s=Na;         %sink   mol/s-m^3
    function C0=pbIC1(x,Xmax,Cw)   %initial condition
    if x==Xmax
      C0=Cw;
    else  
      C0=0.04;   %mol/m^3
    end  
    function [pl,ql,pr,qr]=pbBC1(yl,Cl,yr,Cr,t,Cw)  %boundary conditions
    pl=0;   %left side bc. This one sets flux to zero since the left side is impermeable
    ql=1;       %Form is pl + ql*f = 0 where f is flux from definition of pde system in function pdepb
    pr=Cr-Cw;       %right side bc. This one sets concentration to Cw
    qr=0;
Related Question