Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Unable to meet tolerance using bvp4c
Replies: 0

 Abed Alnaif Posts: 13 Registered: 10/21/08
Unable to meet tolerance using bvp4c
Posted: Nov 2, 2010 1:56 PM

Hello,

I am trying to solve a boundary-value problem - a system of 3 ODEs - using the bvp4c command. I am returned the warning below regarding the tolerance. The warning states that it needs to use 1666 mesh points to meet my specified tolerance (1E-9). However, as the warning states itself, the simulation consists of more mesh points (2004 mesh points) than the required 1666. I also tried the simulation with many more mesh points (20005 mesh points), and I still got a similar warning message. Why, then, is it not able to meet the tolerance? I included the exact warning message, and my code (in case anyone wanted to try this for themselves - takes about 20sec to run on my laptop), below.

Thanks,

Abed

% ------------------------------------------------------------------------------
WARNING MESSAGE:

Warning: Unable to meet the tolerance without using more
than 1666 mesh points.
The last mesh of 2004 points and the solution are
available in the output argument.
The maximum residual is 0.305596, while requested
accuracy is 1e-09.
> In bvp4c at 326
In test_model at 123

% ------------------------------------------------------------------------------
% test_model CODE:
% ------------------------------------------------------------------------------

function test_model
clear all;
tStart=tic; % start timer

% --------------------- Parameters ---------------------------
% Compartment boundaries
X_max = 100E-6; % [m] (100 microns)
X_min = 20E-6; % [m] (20 microns)

x_m = X_min/X_max;

% Diffusion constants
D_L = 1E-11; % [m^2/s]
D_E = 1E-11; % [m^2/s]
D_EL = 1E-11; % [m^2/s]

k_D_E = D_E/D_L;
k_D_EL = D_EL/D_L;

% k_D_E_sink = 1; % use different parameters in sink region
% k_D_EL_sink = 1;

A_X = 3.3E-4; % [1/s]
A_XX = 3.3E-3; % [1/s]
A_L = A_X; % [1/s]
A_R = A_X; % [1/s]
A_E = A_X; % [1/s]
A_LR = A_XX; % [1/s]
A_EL = A_XX; % [1/s]
A_ER = A_XX; % [1/s]

m_a = (X_max^2)/D_L; % nondimensionalizing factor
a_L = m_a*A_L;
a_R = m_a*A_R;
a_E = m_a*A_E;
a_LR = m_a*A_LR;
a_EL = m_a*A_EL;
a_ER = m_a*A_ER;

% a_L_sink = 0.5; % use different parameters in sink region
% a_E_sink = 0.5;
% a_EL_sink = 0.5;

% Production parameters
V_L_max = 1E-4; % [mol/s]
V_E_max = 1E-4; % [mol/s]
V_R_max = 1E-4; % [mol/s]
K_E = 1;
K_R = 1;
n_E = 2;
n_R = 2;
p = 0.2;

R_0_bar = V_R_max/A_R; % initial receptor constant

m_v = (X_max^2)/(D_L*R_0_bar); % nondimensionalizing factor
v_L_max = m_v*V_L_max;
v_E_max = m_v*V_E_max;
v_R_max = m_v*V_R_max;
k_E = K_E*R_0_bar^n_E;
k_R = K_R*R_0_bar^n_R;

% Rate constants
K_p = 1E-2; % [1/(s*M)]
K_m = 1E-5; % [1/s]
K_LR_p = K_p; % [1/(s*M)]
K_LR_m = K_m; % [1/s]
K_EL_p = K_p; % [1/(s*M)]
K_EL_m = K_m; % [1/s]
K_ER_p = K_p; % [1/(s*M)]
K_ER_m = K_m; % [1/s]

m_k = (X_max^2)/D_L;
k_LR_p = m_k*R_0_bar*K_LR_p;
k_LR_m = m_k*K_LR_m;
k_EL_p = m_k*R_0_bar*K_EL_p;
k_EL_m = m_k*K_EL_m;
k_ER_p = m_k*R_0_bar*K_ER_p;
k_ER_m = m_k*K_ER_m;

% Define initial mesh resolution.
nsteps = 2000; % Number of steps

% ---------------------- Calculated Parameters ----------------------
x_s = 5*sqrt(max([(1/a_L) (k_D_E/a_E) (k_D_EL/a_EL)])); % sink location
% x_s = 5*sqrt(max([(1/a_L_sink) (k_D_E_sink/a_E_sink) (k_D_EL_sink/a_EL_sink)])); % sink location
h = (x_s+x_m)/nsteps; % step size

% ---------------------- Initialization ----------------------------
% Define mesh in each region
xinit_region1 = -x_m:h:0;
xinit_region2 = 0:h:1;
xinit_region3 = 1:h:x_s;

% Since a region may not be divisble by h, ensure that the boundaries
% of the region are defined
if( xinit_region1(end) ~= 0 )
xinit_region1 = [ xinit_region1 0 ];
end
if( xinit_region2(end) ~= 1 )
xinit_region2 = [ xinit_region2 1 ];
end
if( xinit_region3(end) ~= x_s )
xinit_region3 = [ xinit_region3 x_s ];
end

% Initialize mesh.
xinit = [ xinit_region1 xinit_region2 xinit_region3 ];

% Constant initial guess for the solution
yinit = [0.1; 0.1; 0.1; 0.1; 0.1; 0.1];

% The initial profile
sol = bvpinit(xinit,yinit);

% --------------------- Simulation ------------------------------------
options_bvp = bvpset('RelTol', 1e-9, 'AbsTol', 1e-12);
sol = bvp4c(@f,@bc,sol,options_bvp);

% --------------------- Visualization ---------------------------------
% debug lines
% sol.x
% sol.y

tCalculateR = tic;
r_mem = zeros(1,length(sol.x));
for i = 1:1:length(sol.x) % loop through all mesh points
if(sol.x(i) <=1)
[r_mem(i),~,exitflag] = fzero(@(r_mem) receptor_equation(r_mem,sol.x(i),sol.y(:,i)), [-0.0000001 1000]); % find r numerically (can't find explicit solution)
if(exitflag ~= 1)
error('Abedsays:rNotFoundFinal', 'r was not found. Values: exitflag: %d, y(1): %f, y(3): %f', exitflag, sol.y(1,:), sol.y(3,:));
end
else
r_mem(i) = 0;
end
end
TimeToCalculateR = toc(tCalculateR)
b_er_mem = (k_ER_p*sol.y(3,:).*r_mem)/(k_ER_m + a_ER);
b_lr_mem = (k_LR_p*sol.y(1,:).*r_mem)/(k_LR_m + a_LR);

% debug lines
% 'r_final'
% r
% 'b_er'
% b_er
% 'b_lr'
% b_lr
figure
plot(sol.x,sol.y(1,:),'b',sol.x,sol.y(3,:),'r',sol.x,sol.y(5,:),'g',sol.x,r_mem,'--k',sol.x,b_er_mem,'--r',sol.x,b_lr_mem,'--b')
legend('l','e','el','r','er','lr');
TimeElapsed = toc(tStart)

% --------------------- Nested Functions ------------------------------
% Numerical solution for receptor concentration (can't find explicit
% solution)
function rec = receptor_equation(r,x,y)
% H = (p^2*heaviside(-x))+heaviside(x)-heaviside(x-1);
H = (p^2*heaviside(-x))+heaviside(x);
rec = ( (H*v_R_max)/(1+(k_R*((k_LR_p*y(1)*r)/(k_LR_m+a_LR))^n_R)) ) + r*( -(k_LR_p*y(1)) + ((k_LR_m*k_LR_p*y(1))/(k_LR_m+a_LR)) - (k_ER_p*y(3)) + ((k_ER_m*k_ER_p*y(3))/(k_ER_m+a_ER)) - a_R );
% sprintf('x: %f, H: %f, l: %f, e: %f, r: %f', x, H, y(1), y(3), r)
% % debug line
end

% Define differential equations
% eventually, this should be vectorized to improve speed - see
% http://www.mathworks.com/help/techdoc/ref/bvpset.html
% cannot be vectorized due to use of fzero function - see
function dydx = f(x,y,~)

if (x <= 1)
% algebraic equations for r, b_er, and b_lr
[r,~,exitflag] = fzero(@(r) receptor_equation(r,x,y), [-0.0000001 1000]); % find r numerically (can't find explicit solution)
if(exitflag ~= 1)
error('Abedsays:rNotFound', 'r was not found. Values: exitflag: %d, y(1): %f, y(3): %f', exitflag, y(1), y(3));
end
b_er = (k_ER_p*y(3)*r)/(k_ER_m + a_ER);
b_lr = (k_LR_p*y(1)*r)/(k_LR_m + a_LR);

% production of species
v_L = v_L_max * heaviside(-x); %ligand production
v_E = (v_E_max/(1+(k_E*(b_lr)^n_E)))*(heaviside(x)); %expander production
else
% sink region
r = 0;
b_er = 0;
b_lr = 0;
v_L = 0;
v_E = 0;
end

% differential equations
% y(1):=l, y(2):=l', y(3):=e, y(4):=e', y(5):=b_el, y(6):=b_el'
dydx = [ y(2)
-v_L + (k_EL_p*y(3)*y(1)) - (k_EL_m*y(5)) + (k_LR_p*y(1)*r) - (k_LR_m*b_lr) + (a_L*y(1))
y(4)
( -v_E + (k_ER_p*y(3)*r) - (k_ER_m*b_er) + (k_EL_p*y(3)*y(1)) - (k_EL_m*y(5)) + (a_E*y(3)) ) / k_D_E
y(6)
( -(k_EL_p*y(3)*y(1)) + (k_EL_m*y(5)) + (a_EL*y(5)) ) / k_D_EL ];
end

% Define boundary conditions
function res = bc(YL,YR)
res = [ YL(2,1) % dl(-xmin)/dx = 0
YR(1,end) % l(xs) = 0
YL(4,1) % de(-xmin)/dx = 0
YR(3,end) % e(xs) = 0
YL(6,1) % d(b_el(-xmin))/dx = 0
YR(5,end) % b_el = 0
YR(1,1)-YL(1,2) % continuity of l at x=0
YR(2,1)-YL(2,2) % continuity of l' at x=0
YR(3,1)-YL(3,2) % continuity of e at x=0
YR(4,1)-YL(4,2) % continuity at e' at x=0
YR(5,1)-YL(5,2) % continuity of b_el at x=0
YR(6,1)-YL(6,2) % continuity of b_el' at x=0
YR(1,2)-YL(1,3) % continuity of l at x=1
YR(2,2)-YL(2,3) % continuity of l' at x=1
YR(3,2)-YL(3,3) % continuity of e at x=1
YR(4,2)-YL(4,3) % continuity at e' at x=1
YR(5,2)-YL(5,3) % continuity of b_el at x=1
YR(6,2)-YL(6,3) ]; % continuity of b_el' at x=1
end
end