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."
15 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
Torsten
on 16 Nov 2018
If the code works for values of L0 smaller than 100, try approaching 100 by subsequently calling bvp4c with values L0start, L01, L02,...,L0=100.
Take a look at
https://de.mathworks.com/matlabcentral/answers/428271-how-to-write-matlab-bvp4c-code-for-moving-boundary-surface-problem
to see how this can be implemented.
Best wishes
Torsten.
0 Comments
See Also
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!