Error using bvp4c: "Warning: Unable to meet the tolerance without using more than 5000 mesh points. The last mesh of 4921 points and the solution are available in the output argument."

Question

Jennifer Yang 2018-11-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/430112-error-using-bvp4c-warning-unable-to-meet-the-tolerance-without-using-more-than-5000-mesh-points

回答： Torsten 2018-11-16

Hello,

I'm using bvp4c to solve for a reaction-diffusion equation. When I set the initial concentration to 100, I get this error saying

"Warning: Unable to meet the tolerance without using more than 5000 mesh points.
 The last mesh of 4921 points and the solution are available in the output argument.
 The maximum residual is 0.512336, while requested accuracy is 0.001. 
> In bvp4c (line 323)
  In steady_state_bvp4c_FINAL (line 38) 
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In steady_state_bvp4c_FINAL (line 44) 
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In steady_state_bvp4c_FINAL (line 50) "

Am I setting my boundary conditions incorrectly? I want the left hand side (at x = 0, C_L = C_L0) and on the right hand side I want the flux to equal to 0 (at x = x_f, dC_L/dx = 0).

function steady_state_bvp4c_FINAL
clear all; clc;
%Diffusion coefficient
D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s
%Initial concentration at x=0
L0 = 100;
%Total length
x_f = 3500;
%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;
%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;
%Total Receptors
N_T =1.7;
L = L0./2;
solinit = bvpinit(linspace(0,x_f,11),[L 0]);
sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);
figure(1)
plot(sol.x,(sol.y(1,:)),'LineWidth',1)
xlabel('Distance (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])
figure(2)
plot(sol.x,(sol.y(1,:))./L0,'LineWidth',1)
xlabel('Distance (\mum)')
ylabel('Normalized Concentration (C_L/C_L_0)')
axis([0 x_f 0 1])
function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
C_L    = y(1);
dC_Ldx = y(2);
N_r = (-1 - (C_L./K_1) + sqrt((-1 - (C_L./K_1)).^2 -4.*(-(2./K_4) - ((2.*C_L)./(K_2.*K_4)) - ((2.*C_L)./(K_2.*K_4.*K_i)) - ((2.*C_L.^2)/(K_2.*K_3.*K_4.*K_i))).*N_T))/(4.*((1./K_4) + (C_L./(K_2.*K_4)) + (C_L./(K_2.*K_4.*K_i)) + (C_L.^2./(K_2.*K_3.*K_4.*K_i))));
N_R = (C_L./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;
r = N_T-(2.*N_R2)-(4.*N_pR)-N_R-(4.*N_rR)-(2.*N_r2);
R = N_T-(2.*N_R2)-(4.*N_pR)-N_r-(4.*N_rR)-(2.*N_r2);
r_2 = (N_T./2)-(N_R./2)-N_R2-(2.*N_pR)-(N_r./2)-(2.*N_rR);
rR = (N_T./4)-(N_R./4)-(N_R2./2)-N_pR-(N_r./4)-(N_r2./2);
pR = (N_T./4)-(N_R./4)-(N_R2./2)-(N_r./4)-N_rR-(N_r2./2);
R_2 = (N_T./2)-(N_R./2)-(2.*N_pR)-(N_r./2)-(2.*N_rR)-N_r2;
R_L = (2.*d_1.*(R))-(2.*k_1.*C_L.*r)+...
    ((2.*d_2.*(rR)))-((k_2.*C_L.*(r_2)))+...
    (d_3.*(R_2))-(2.*k_3.*C_L.*(pR));
dy = zeros(2,1);
dy(1) = dC_Ldx;
dy(2) = (-R_L/D_ij);
end
function res = twobc(ya,yb)
res = [ ya(1)-(L0); yb(2) ];
end
end