pdepe: Unable to meet integration tolerances without reducing the step size below the smallest value allowed

Hi, I'm trying to solve a system of transient convection diffusion reaction equations composed of 5 identical format PDEs. The function is below:

where d is the diffusivity, v is the velocity, G is the reaction term. The system is in cylindrical coorcinates. The problem is when if I remove the convection term in my code, it worked well. But if I added the convection term back in, the error came out: Warning: Failure at t=3.656244e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.

I don't know why this would happen. Could someone please help me? Thank you very much! Below is my code:

%main calling function

function model

m = 1;

x = linspace(0, 2.8, 29);

t = linspace(0, 200, 201);

sol = pdepe(m, @model_pde, @model_ic, @model_bc, x, t);

u_Ia = sol(:,:,1);

figure

surf(x, t, u_Ia)

title ('Numerical solution of Ia')

xlabel('Distance x')

ylabel ('Time t')

figure

[C, h] = contour(t, x, u_Ia', [350, 350], 'b')

xlabel('Time (min)')

ylabel('Clot Radius (mm)')

end

%boundary conditions

function [pl, ql, pr, qr] = model_bc(xl, ul, xr, ur, t)

pl = [0; 0; 0; 0; 0];

ql = [1; 1; 1; 1; 1];

pr = [0; 0; 0; 0; 0];

qr = [1; 1; 1; 1; 1];

end

%initial conditions

function u0 = model_ic(x)

u0 = [7000*(1-heaviside(x-2.5)); 2.18*heaviside(x-2.5); 2180*heaviside(x-2.5); 447.7*heaviside(x-2.5); 105*heaviside(x-2.5)];

end

%pde function

function [c, f, s] = model_pde(x, t, u, DuDx)

h_1 = 1500;

H_1m = 250000;

k_pla_upa = 60;

K_plam_upa = 50000;

h_pla = 0.096;

flow_rate = 10*10^3;

v = flow_rate / (pi*(2.4^2-(x/2)^2));

t_initial = 20;

a = 5*10^-2;

b = 5*10^-3;

m = 9*10^-6;

d = 5*10^-8;

e = 1;

rate = 1+a*v^1+b*v^2+m*v^3;

logi = 1/(1+exp(-(t - t_initial/(1+d*v^e))));

d_ia = 1.482*10^-3*logi;

d_pla = 2.958*10^-3*logi;

d_pls = 2.886*10^-3*logi;

d_l2ap = 3.15*10^-3*logi;

d_upa = 4.446*10^-3*logi;

G_Ia =-h_1*u(2)*u(1)/(H_1m+u(1));

G_PLA_generate = k_pla_upa*u(4)*u(3)/(K_plam_upa+u(3));

G_PLA_deplete = -h_pla*u(2)*u(5);

G_PLS = -k_pla_upa*u(4)*u(3)/(K_plam_upa+u(3));

G_uPA = 0;

G_L2AP = -h_pla*u(2)*u(5);

c = [1; 1; 1; 1; 1];

s = [G_Ia*rate*logi; G_PLA_generate*rate*logi + G_PLA_deplete*logi; G_PLS*rate*logi; G_uPA*logi; G_L2AP*logi];

f = [d_ia; d_pla; d_pls; d_upa; d_l2ap] .* DuDx - [v*u(1); v*u(2); v*u(3); v*u(4); v*u(5)];

% f = [d_ia; d_pla; d_pls; d_upa; d_l2ap] .* DuDx; %if there were no convection term, the code worked well

end

Answer 1

[Bill Greene]

Hi,

pdepe is designed for PDEs where the diffusion term is relatively large compared with the convection term. In your equations the convection terms are much larger than the diffusion terms. A key parameter is the cell Peclet number which is c*h/(2*d) where c is the convection coefficient, d is the diffusion coefficient, and h is the size of each cell (element) in the mesh. For the algorithm to be stable this number should be less than one. You can find discussion of this in books on numerical methods for convection-diffusion equations. (e.g. page 508 in http://www.amazon.com/Computational-Science-Engineering-Gilbert-Strang/dp/0961408812).

Some approaches for solving this class of problem add some "artificial" diffusion-- i.e. a larger diffusion coefficient than actually exists. You might try this and see how that affects the solution. Also, you have a relatively coarse mesh with only 29 elements along the 2.8 unit length. In general, the mesh density affects accuracy and, in this case, also the stability of the method. I recommend experimenting with a single, simple convection diffusion equation to get a feel for how these parameters affect the solution.

Finally, I noticed one thing about your problem that seems odd. The length is 2.8 units, the velocity is on the order of 550, and you are solving the equations for a time span of 200. It seems like whatever is being convected will be transported well off the model long before you reach the final time?

Hope this helps.

Bill

[Zhen]

Hi Bill,

Thank you for your reminder. I also thought of it might be the problem of dimension scale, but didn't take action. I just tried to reduce the velocity to 10*10^-3. It worked, but some odd thing came out aroung x=0 and t=0. Anyway the problem is, the units of length and velocity are in mm scale, and the units of those property parameters are in nM scale. The original flow rate was 10 ml/min. I have to transfer it to mm^3/min before proceeding. That's why I put 10*10^3 in the flow rate, and it crashed.

I don't quite get your point of adding "artificial" diffusion coefficient. The diffusion coefficients I currently use come from literatures and have been transferred to mm/min scale. If I add some other large diffusion coefficient, how can I transfer back to the real one? I tried to reduce the mesh size to 10 times finer than before (281 elements), but it still didn't work at the velocity of 10*10^3. It worked at 10*10^-3, but the odd thing became more serious.

The time span was just a safety region. I expected that by varying velocity, which was a controlled paramter, the final time of u_Ia reached 0 also the model length would also change. Without the convection term, the final time would be ~120 min.

Your suggestion really helps me a lot. Thank you so much! I would be very appreciated if you can further explain to me the questions I have in your answer.

Zhen

[Zhen]

The diffusion coefficient I transfered was in mm^2/min. Sorry for the typo above.

[Bill Greene]

>If I add some other large diffusion coefficient, how can I transfer back to the real one? Introducing artificial diffusion definitely introduces some error in the solution so the idea is to introduce as little as possible. That is the major defect to this approach.

If you want to understand more about the difficulties of solving this class of problem and the appropriate numerical techniques, I suggest this book:

http://www.amazon.com/Methods-Hyperbolic-Problems-Cambridge-Mathematics/dp/0521009243

Bill

Answer 2

[tom_brot]

If you still have stability issues, check my answer at

http://www.mathworks.com/matlabcentral/answers/56037-solving-one-dimensional-pde-s-using-the-pde-toolbox

BR, Tom