PDEPE: Unable to meet integration tolerances

Hi,

I am trying to solve a system of 2 PDEs using Matlab's built-in PDE solver, pdepe. I am being returned an "Unable to meet integration tolerances" warning during the ODE15s routine -- see the bottom of this message to see the specific warning. I read that this warning my occur if the equation has singularities, but I didn't notice any such singularities in my equations.

Below is my code, in a self-contained file. I thought that the problem would be either with the p.h parameter being to large, or that it may be due to the presence of the M^2 term. However, I still get the same warning message even after setting p.h = 1, as well as after removing the M^2 term.

Needless to say, any help is greatly appreciated!

Matlab code:

function unable_to_meet_tolerance()
      clear all;

    % -------- Model Parameters ----------------------
        p.x_m = 100;
        p.D_M = 10;
*        p.D_E = 1;
        p.a_M_1 = 1E-1;
        p.a_M_2 = 1;
        p.a_E = 1E-5;
        p.b_E = 1E-2;
        p.T_p2 = 1E-3;
        p.h = 2; 
        p.j = 10; % flux

    % -------- Model ---------------------------------
        % PDE function
        function [c,f,s] = pde(~,~,u,dudx)
            M = u(1);
            E = u(2);

            % differential equations
            % see Matlab PDE solver help file for definitions of c, f, and s
            c = [1; 1];
            f = [p.D_M; p.D_E] .* dudx;
            s = [ - ( (1+E) * p.a_M_1 * M ) - ( (1+E) * p.a_M_2 * M^2 )
                - ( p.a_E * E ) + ( p.b_E/(1+(M/p.T_p2)) )];
        end

        % boundary conditions
        function [pl,ql,pr,qr] = pde_bc(~,~,~,~,~)
            pl = [-p.j; 0];
            ql = [1; 1/p.D_E]; 
            pr = [0; 0];
            qr = [1/p.D_M; 1/p.D_E]; 
        end

        % initial conditions
        function u0 = pde_ic(~)
            u0 = [0; 0]; % zero initial conditions
        end

    % -------- Simulation ----------------------------
        nsteps_x = 1000;
        nsteps_t = 5;
        t_fin = 1;

        x_mesh = linspace(0,p.x_m,nsteps_x+1);
        t_mesh = linspace(0,t_fin,nsteps_t+1);

        sol = pdepe(0,@pde,@pde_ic,@pde_bc,x_mesh,t_mesh);
end

This is the warning message:

Warning: Failure at t=4.060671e-03.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.387779e-17) at time t. 
> In ode15s at 819
  In pdepe at 320
  In unable_to_meet_tolerance at 53
Warning: Time integration has failed. Solution is available at requested time points
up to t=0.000000e+00. 
> In pdepe at 326
  In unable_to_meet_tolerance at 53

[Abed Alnaif]

I forgot to mention, as it may or may not be clear from the code above, I am going for zero initial conditions and no-flux boundary conditions (with the exception of M, which has a flux of p.j at x=0).

Answer 1

[Grzegorz Knor]

Are you sure that boundary conditions are correct?

Maybe one minus is unnecessary?

        % boundary conditions
        function [pl,ql,pr,qr] = pde_bc(~,~,~,~,~)
            pl = [p.j; 0];
            ql = [1; 1/p.D_E]; 
            pr = [0; 0];
            qr = [1/p.D_M; 1/p.D_E]; 
        end

[Abed Alnaif]

Grzegorz, you are correct (an embarrassing mistake on my part!). Thanks for your help; it works now after I made the change that you suggested!