Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered.

Hello

I'm trying to solve an ODE, but getting the following error:

Error using bvp4c (line 251)

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

Here is my code:

Diff=7.0000e-11;
koffPart=5.1910e-06;
konPart=6.7702e-18;
CT=1.2044e+17;
LT=1.2044e+16;

xmesh = linspace(0,1000,100);
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@(t,C) odefcn(t,C,Diff,konPart,koffPart,CT,LT), @(Ca,Cb) bcfcn(Ca,Cb,CT), solinit);

Here is my ODE function:

function dCdt = odefcn(t,C,Diff,konPart,koffPart,CT,LT)
    dCdt=zeros(2,1);
    dCdt(1)=C(2);
    dCdt(2)=(C(2)+konPart*C(1)*(LT-CT+C(1))-koffPart*(CT-C(1)))/Diff;
end

Here are my boundary conditions:

function res = bcfcn(Ca,Cb,CT)
res = [Ca(1)-CT
       Cb(2)];
end

And my guess:

function g = guess(x)
g = [exp(-x)
     -exp(-x)];
end

I'd really appreiate any help.

Thanks!

Ben

[Torsten]

Your boundary condition for Ca(1) at the left endpoint does not match with your initial condition for C(1).

Try res(1) = Ca(1) - 1 to test.

[Ben Miller]

Hi Torsten

Thanks for your answer. I'm trying to have initial condition C(1)=CT. Not sure if I've coded that right?

I tried res(1)=Ca(1)-1, but got same error.

Ben

[Torsten]

Simplify your model (moderate model constants in size, smaller integration interval,...) until the error vanishes.

[Muhammad Safdar Nadeem]

function newpapercode

m=0;

zeta=0.5;

xi=1;

rho_s=997;

rho_f=4179;

M=30;

C_T=0;

P_rf=0.1;

E_c=0.5;

P_rf=1;

K_f=1;

K_fn=.1;

n=0.006;

eta=0.5;

lamda=1;

rho_cps=993;

rho_cpf=0.001;

gamma=1;

omega=1;

F_n=1;

lambda_u=2;

K_f=0;

K_fn=1;

N=1;

S=1;

sol = bvpinit(linspace(0,4,100),[0 0 0 0 0 ]);

sol1 = bvp4c(@bvpexam2, @bcexam2, sol);

x1 = sol1.x;

y1 = sol1.y;

figure (1)

plot(x1, y1(2,:));

hold on

% figure (2)

% plot(x1, y1(4,:));

% hold on

% figure (3)

% plot(x1, y1(6,:));

% hold on

% figure (4)

% plot(x1, y1(8,:));

% hold on

z = y1(3,1);

p = y1(7,1);

function res = bcexam2(y0, yinf)

res=[y0(1)+2*S/m+1;y0(2)-0;y0(4)-1;yinf(2)-1;yinf(4)-0];

end

function ysol = bvpexam2(~,y)

yy1 = ((2/m+1)*(1-zeta)^2.5*(xi)*(xi)/(1-zeta+zeta*(rho_s/rho_f)))*((M+lamda/(1-zeta)^2.5)*(y(2)-1)+(xi/1-m)*(1-zeta+zeta*(rho_cps/rho_cpf))*(gamma)*(sin(omega/2)))+(1-zeta+zeta*(rho_s/rho_f))*(1-zeta)^2.5*(((2)*(m-F_n)*(y(2)*y(2)-1))/(m+1)+y(1)*y(3)-lambda_u*(2-eta*y(3)-2*y(2)));

yy2 = (-4/3)*(K_f/K_fn)*(N)*((C_T+y(4))^3*y(5))-((P_rf)*(E_c)*y(3)*y(3))/(1-zeta)^2.5+(P_rf)*(1-zeta+zeta*(rho_cps/rho_cpf))*(K_f/K_fn)*(((2)*(n)*(y(2))*(y(4)))/(m+1)-(y(1)*y(5))+(lambda_u*eta*y(5)));

ysol = [y(2);y(3);yy1;y(5);yy2;];

end

end

[Muhammad Safdar Nadeem]

please help me in this error

Error using bvp4c (line 251)

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

Error in newpapercodeH (line 28)

sol1 = bvp4c(@bvpexam2, @bcexam2, sol);

[Torsten]

You use

res=[y0(1)+2*S/m+1;y0(2)-0;y0(4)-1;yinf(2)-1;yinf(4)-0]; 

in your code, but m = 0. So S/m means that you divide by 0.

Maybe you meant S/(m+1) instead of S/m+1.

Answer 1

[Harshendra Shah]

Hi Ben,

A singular Jacobian indicates that the initial guess causes the solution to diverge. The BVP4C function finds the solution by solving a system of nonlinear algebraic equations. Nonlinear solvers are only as effective as the initial guess they start with, so changing your starting guess may help. Also, BVP4C is responsible for providing a guess for the solution from one iteration to another.

This issue can occur due to multiple reasons. You can find the same in this MATLAB Answers page: https://www.mathworks.com/matlabcentral/answers/120468-why-do-i-receive-an-error-about-a-singular-jacobian-when-i-use-the-bvp4c-function-within-matlab

This MATLAB Answer also suggests some workarounds to solve this issue.

I hope this helps.

Answer 2

[Hifza]

hifza()

Error encountered while solving the BVP:
Unable to solve the collocation equations -- a singular Jacobian encountered.

function hifza
format long g

% ==== PARAMETERS ====
delta1 = 0.5; lamda = 0.01; sinr = 0.1; M = 0.5; epsilon = 0.01;
Rd = 0.1; Q = 0.01; Nb = 0.01; Nt = 0.01; Le = 0.3;
Bi = 0.1; B = 0.1; epsilon1 = 0.01; delta2 = 0.5; delta3 = 0.1;
delta4 = 0.1; n = 0.1; E = 1;

% Adding Prandtl number (Pr)
Pr = 0.7; % Typical value for air; adjust as needed

% ==== INITIAL MESH AND GUESS ====
xmesh = linspace(0, 10, 200); % Reduce mesh points for better stability
init_guess = @(x) [...
    0.5 * sin(x) + 0.1;  
    exp(-x) + 0.1;       
    0.1;                 
    cos(x) + 0.1;        
    -exp(-x) / 3 + 0.1;  
    sin(x) + 0.1;        
    -exp(-x) / 4 + 0.1   
    ];
solinit = bvpinit(xmesh, init_guess);

% ==== SOLVER OPTIONS ====
options = bvpset('RelTol', 1e-5, 'AbsTol', 1e-7, 'NMax', 200000, 'Stats', 'on');

% ==== SOLVE BVP ====
try
    sol = bvp4c(@bvpexam2, @bcexam2, solinit, options);
catch ME
    disp('Error encountered while solving the BVP:');
    disp(ME.message);
    return;
end

% ==== EXTRACT SOLUTION ====
x = sol.x;
y = sol.y;

% ==== PLOT ====
figure;
plot(x, y(2,:), 'r-', 'LineWidth', 2);
xlabel('x'); ylabel('f(x)'); title('Velocity Profile'); 
grid on;

% ==== VALUES AT x = 0 ====
value = deval(sol, 0);
disp('Solution at x = 0:');
disp(vpa(value, 6));

% ==== SHERWOOD AND NUSSELT NUMBERS ====
sherwood = real(y(7,1));
nusselt = real((1 + (4/3)*Rd) * y(5,1) * sqrt(Pr)); % Including Prandtl number effect
fprintf('Sherwood number: %.6f\n', sherwood);
fprintf('Nusselt number: %.6f\n', nusselt);
fprintf('Prandtl number: %.6f\n', Pr); % Displaying Pr value

% ==== ODE SYSTEM ====
    function dydx = bvpexam2(~, y)
        f = y(1); fp = y(2); fpp = y(3);
        theta = y(4); thetap = y(5);
        phi = y(6); phip = y(7);
        
        den1 = max(1e-2, 1 + lamda * fpp); % Avoid singularity
        den2 = max(1e-2, 1 + (4/3)*Rd + epsilon1 * theta); % Stability check
        %den1 = 1 + lamda * fpp;
        %den2 = 1 + 4/3*Rd + epsilon1 * theta;
        
        theta_safe = max(-0.99, theta);
        %theta_safe = theta;
        T = delta4 * (1 + theta_safe);
        exp_term = exp(-E / T);
        exp_term = min(exp_term, 1e3); % Cap exponentiation
        
        yy1 = (1 / den1 + (1 + 1/B) + epsilon * (1 + delta1 * fpp^2)) * ...
            (-f * fpp + fp * (M * sinr + fp));
        yy2 = -(Pr / den2) * (-Nt * thetap^2 / 5 - Nb * thetap * phip - thetap * f - Q * theta - epsilon1 * thetap^2);
        yy3 = (Nt / Nb) * yy2 - Le * f * phip + Le * delta2 * (1 + theta * delta3^n) * exp_term * phi;
        
        dydx = [fp; fpp; yy1; thetap; yy2; phip; yy3];
    end

% ==== BOUNDARY CONDITIONS ====
    function res = bcexam2(ya, yb)
        res = [ya(1); 
            ya(2) - 1;
            ya(5) + Bi * (ya(4) - 1);
            ya(6) - 1;
            ya(7) + 0.5;
            yb(4);
            yb(6)];
    end
end

this is my code when run the code the ''Error encountered while solving the BVP:

Unable to solve the collocation equations -- a singular Jacobian encountered'' . please guide me how can remove the error in this code

[Torsten]

It's hard to say what the problem is.

Maybe your system does not have a solution. Maybe the initial guesses for the solution variables are too bad. Maybe the discontinuities you introduce into the equations with your max and min operations hinder convergence.