bvp4c - Unable to solve the collocation equations -- a singular Jacobian encountered.

7 次查看（过去 30 天）
Wes Fermanich
Wes Fermanich 2020-6-26
编辑： Wes Fermanich 2020-7-1
Hello,
I am working on solving a system of 3 ODEs using the bvp4c solver. I am receiving the common error "Unable to solve the collocation equations -- a singular Jacobian encountered." From what I have read, this is an error that can be due to a number of problems with the the equation function, boundary value function, or the initial guess. I have tried altering a few things but to no avail.
The .pdf attached details the system I am working with as well as the boundary conditions. I am working to recreate a figure from a research paper and have attached that figure which should give an idea what the solution hopefully will look like.
I believe that the error is likely in my initial guess as I am not quite sure on how to best create that initial guess. I have tried to use functions for the guess that roughly match the solutions plotted in the figures I mentioned above.
Any help or insight is greatly appreciated!
**update July 1st: found a sign error in the first boundary condition and flipped the subraction to addition
global d alpha phi D_0 D_I D_I3 k_edr k_eer k_rg k_ext I2_0 I_0 Keq b gamma theta beta m1 n_0 I3_0 e0
e0 = 1.602e-19; % elementary charge (Coulombs)
d = 12.5e-4; % cm
alpha = 0.2e-4; % cm^-1
phi = 6.1e16; % cm^-2 s^-1
D_0 = 0.4; % cm^2 s^-1
D_I = 4.4e-6; % cm^2 s^-1
D_I3 = 3.6e-6; % cm^2 s^-1
k_edr = 3e-11; % cm^3 s^-1
k_eer = 2e-16; % cm^3 s^-1
k_rg = 1.3e-15; % cm^3 s^-1
k_ext = 10^12; % cm s^-1
n_0 = 1e14; % cm^-3
I2_0 = 0.05; % M
I_0 = 0.6; % M
Keq = 10^8.3;
b = 1;
gamma = 1.4; % gamma/omega^2
theta = 0.45;
beta = 0.35;
m1 = 1.55;
% solving for the equilibrium concentrations of I3_0 and I_0
syms I3
eqn = solve(I3/((I_0 - I3)*(I2_0 - I3)) == Keq, I3);
I3_0 = double(eqn(1));
I_0 = I_0 - I3_0;
x = linspace(0,d,51);
solinit = bvpinit(x,@guess);
sol = bvp4c(@bvfcn,@bc,solinit);
%---Functions----------------------------------------
% u(1) = n, u(2) = dn/dx, u(3) = [I-], u(4) = d[I-]/dx, u(5) = [I3-],
% u(6)= d[I3-]/dx
%
% dudx(1) = dn/dx, dudx(2) = d2n/dx2, dudx(3) = d[I-]/dx,
% dudx(4) = d2[I-]/dx2, dudx(5) = d[I3-]/dx, dudx(6) = d2[I3-]/dx2
function dudx = bvfcn(x,u)
global phi alpha k_eer k_edr k_rg gamma theta D_0 D_I D_I3
S = phi.*alpha.*exp(-alpha.*x)./(k_rg.*u(3)+k_edr.*u(1));
dudx = [u(2)
(-phi.*alpha.*exp(-alpha.*x)+k_eer.*u(5).*u(1)+k_edr.*S.*u(1))./(gamma.*(1-theta).*D_0)
u(4)
(-1.5.*k_eer.*u(5).*u(1) + 1.5.*k_rg.*S.*u(3))./(gamma.*theta.*D_I./(1.5))
u(6)
(0.5.*k_eer.*u(5).*u(1) + 0.5.*k_rg.*S.*u(3))./(gamma.*theta.*D_I3./(1.5))];
end
function res = bc(ul,ur)
global e0 theta D_0 k_ext gamma
j = e0.*k_ext.*ul(1);
res = [ul(2)+(j./(e0.*(1-theta).*D_0.*gamma))
ur(2)
ul(4)
ur(4)+(1.5.*j./(e0.*(1-theta).*D_0.*gamma))
ul(6)
ur(6)-(0.5.*j./(e0.*(1-theta).*D_0.*gamma))];
end
function g = guess(x)
global n_0 I_0 I3_0
g = [n_0.*sqrt(x.*10^3)
n_0./sqrt(x.*10^3)
I_0.*(2+cos(x.*10^3+pi))
sin(x.*10^3)
I3_0.*cos(x.*10^3)
-sin(x.*10^3)];
end
  2 个评论
Alan Stevens
Alan Stevens 2020-6-26
There seems to be an inconsistency between the right side boundary conditions and the graph. The slope of I- at x = d is <= 0 (according to the constants specified in your program), yet it is clearly positive on the graph. The opposite is true of I3-. Am I missing something here?
Wes Fermanich
Wes Fermanich 2020-6-28
Ah yes sorry, I should have clarified. In plot a, d corresponds to the dashed line where x = 12.5 um. So for this instance I am looking just at the part of the graph that is left of the dashed line. It doesn't look readily true that the slope of I- and I3- are equal to zero at those points but I believe that there is supposed to be type of inflection-like point near the dashed line for both lines.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by