# BVP4C Unable to solve the collocation equations -- a singular Jacobian encountered

44 visualizaciones (últimos 30 días)
Luke Maran el 25 de Abr. de 2021
Comentada: MOSLI KARIM el 8 de Oct. de 2022
Hello, I am receiving this error attempting to solve a 2nd order ODE with BVP4c.
I am solving for the internal temperature distribution of a nuclear fuel rod using the ode:
with Boundary Conditions:
where R = 0.2m and is the outer diameter
Let me know if you need anything else.
q and k are constants
r is a radius variable
%These are my setup functions.
function dydx = rod_BV(r,y) % Details ODE to be solved
dydx = zeros(2,1);
dydx(1) = y(2)
dydx(2) = -(q/k)-((1/r)*y(1)) % This equation is invalid at r = 0
end
function res = rod_BC(ya,yb) % Details boundary conditions
res = zeros(2,1);
res(1) = ya(2);
res(2) = yb(1) - T_s;
end
% here is my implementation of BVP4C
solinit = bvpinit(linspace(0,R,4),[0,0]); % I have tried setting these 0 values to realmin
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@rod_BV,@rod_BC,solinit,options)
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Respuestas (2)

Divija Aleti el 24 de Jun. de 2021
Hi Luke,
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.
A detailed explanation and possible solutions for this error are given in the link below:
Also, based on the picture of the equation you are solving, I believe the rod_BV function should be written as follows:
function dydx = rod_BV(r,y) % Details ODE to be solved
dydx = zeros(2,1);
dydx(1) = y(2)
dydx(2) = -(q/r)-((1/r)*y(2)) % This equation is invalid at r = 0
end
This is because T = y(1), T' = y(2), which implies, T'' = -(q/r)-(1/r)*y(2).
Hope this helps!
Regards,
Divija
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

MOSLI KARIM el 31 de Jul. de 2022
Editada: Torsten el 31 de Jul. de 2022
Hi, here is the solution that I propose to you
%These are my setup functions.
R=0.2;
q=1;
k=1;
eps=1.0000e-06;
r=linspace(0,R,4);
solinit = bvpinit(linspace(0,R,4),[0,0]); % I have tried setting these 0 values to realmin
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)rod_BV(r,y,q,k),@rod_BC,solinit,options);
y=deval(sol,r);
plot(r,y(1,:))
function dydx = rod_BV(r,y,q,k) % Details ODE to be solved
dydx = zeros(2,1);
dydx(1) = y(2);
dydx(2) = -(q/k)-((1/(r+eps))*y(1)) ; % This equation is invalid at r = 0
end
function res = rod_BC(ya,yb) % Details boundary conditions
%res = zeros(2,1);
% res(1) = ya(2);
% res(2) = yb(1) - 200;
res=[ya(2);yb(1)-200];
end
% here is my implementation of BVP4C
% solinit = bvpinit(linspace(0,R,4),[0,0]); % I have tried setting these 0 values to realmin
% options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
% sol = bvp4c(@rod_BV,@rod_BC,solinit,options)
##### 2 comentariosMostrar NingunoOcultar Ninguno
Torsten el 31 de Jul. de 2022
I think it must be
dydx(2) = -(q/k)-((1/(r+eps))*y(2)) ; % This equation is invalid at r = 0
dydx(2) = -(q/k)-((1/(r+eps))*y(1)) ; % This equation is invalid at r = 0
MOSLI KARIM el 8 de Oct. de 2022
yes sir torsten I forget to pute y(2)

Iniciar sesión para comentar.

### Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

### Community Treasure Hunt

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

Start Hunting!

Translated by