Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.03573e-39

I'm trying to solve an boundary value problem with eigenvalues, following the article ( https://doi.org/10.1063/5.0155331 )

where α and At are certain numbers and Q(σ) is a known function, erf function in this case.

ϕ is the unknown function and Ψ the eigenvalue.

The boundary condition is , so I specified

For small alphas my code with bvp5c works well.

However for large alphas, my solution deviates from the result in the article, and I'm encountering messages like

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 2.03573e-39

I guess the boundary values may be too small and tried to scale my boundary values with C*exp(-3/alpha) but the problem remains.

How can I solve the problem? Fig.1 in https://doi.org/10.1063/5.0155331 can be referred as a correct solution to this problem. My solution looks similar to it but deviates in asymptotic values.

eigenvalues0999 = zeros(1, 30);

global alpha;
global At;
options = bvpset('RelTol',1e-5,'AbsTol',1e-5);

At = 0.999;
for i=1:30    
alpha=1/i;
Psi1 = 1;
solinit = bvpinit(linspace(-3,6,30000),@mat4init,Psi1);
sol = bvp5c(@mat4ode, @mat4bc, solinit,options);
 eigenvalues0999(1,i) = sol.parameters;
end

figure(Color="white");
plot(eigenvalues0999);
grid on;
legend('0','0.243','0.5','0.843','0.999','FontSize', 14)
xlabel('\alpha^{-1}','FontSize', 14);
ylabel('\Psi','FontSize', 14);
set(gca, 'FontSize', 14);

function dydx = mat4ode(x,y,Psi1) % equation being solved
global alpha;
global At;
dydx = [y(2)
        (-alpha^2*At*2/sqrt(pi)*exp(-x.^2)*y(2)+(1+At*erf(x)-alpha*Psi1*2/sqrt(pi)*exp(-x.^2))*y(1))/(alpha^2*(1+At*erf(x)))];
end

function res = mat4bc(ya,yb,Psi1) % boundary conditions
global alpha;
res = [ya(2)-1/alpha*exp(-3/alpha)
       abs(yb(1))+abs(yb(2))
       %yb(2)
       ya(1)-exp(-3/alpha)];
end

function yinit = mat4init(x) % initial guess function
global alpha;
yinit = [exp(-abs(x)/alpha)
        -sign(x)*1/alpha*exp(-abs(x)/alpha)];
end

Answer 1

[Satwik]

Hi,

I understand that you are attempting to solve the problem described in the article using the 'bvp5c' solver. You are encountering a warning and obtaining incorrect results for larger values of 'alpha'. I assume that you are achieving correct results for smaller values of 'alpha' and have already tried scaling the boundary values. Additionally, I am not certain what constitutes large values of 'alpha' in this context, so I assume values above 100 to be large.

In my knowledge, the warning message you are encountering suggests that the Jacobian matrix in your boundary value problem solver is becoming ill-conditioned, which can lead to inaccurate results. This issue often arises in stiff problems or when the parameters of the problem cause the solution to vary rapidly.

Here are a few suggestions to help mitigate this issue:

1. Refine the Mesh: Increase the number of points in the mesh to provide a better initial guess and more resolution for the solver.
2. Adjust Tolerances: Tighten the tolerances in bvpset to improve accuracy.
3. Initial Guess: Provide a better initial guess for the solution to help the solver converge more reliably.

Here's an example of how you might modify your code to refine the mesh and adjust tolerances :

eigenvalues0999 = zeros(1, 30);

global alpha;
global At;
options = bvpset('RelTol',1e-7,'AbsTol',1e-7,'NMax',50000);

At = 0.999;
for i=1:3000    
    alpha = 1/i;
    Psi1 = 1;
    solinit = bvpinit(linspace(-3, 6, 50000), @mat4init, Psi1);
    sol = bvp5c(@mat4ode, @mat4bc, solinit, options);
    eigenvalues0999(1, i) = sol.parameters;
end

The modifications made to the code are :

1. Refined Mesh: Increased the number of points in the mesh from 30,000 to 50,000.
2. Tighter Tolerances: Changed the relative and absolute tolerances in bvpset from 1e-5 to 1e-7.
3. Increased NMax: Increased the maximum number of mesh points (NMax) to allow the solver more flexibility.

I hope this gives you a direction for taking the next steps to resolve the issue.