Singularity Problem trying to Solve a higher order ODE with boundary conditions by bvp4c

Hi, i’m trying to replicate the numerical solution of a recent paper "A macroeconomic model with
Heterogeneous and financially constrained intermediaries. Wouters, R (2019)"
(https://www.econstor.eu/bitstream/10419/207747/1/1067667695.pdf) but after many runs and different
modifications i always get the following issue:

elbarnew= 0.25; %(JCSC) Un valor de salida, siguiendo el rango prob de L&W (0.65,1.66)/(0.57, 3.02)
qlbarnew = 0.05; %(JCSC) Azar corto plazo, escoger inferior frente a qinf
eubar = 3; %(JCSC) Máximo valor de busqueda de e
KK = 500;
solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction_rw);
solF = bvp4c(@Fvecwithproduction_rw,@Fbcwithproduction_rw,solinitF,options);

Error using bvp4c
Unable to solve the collocation equations — a singular Jacobian encountered.

The problem is to solve a nonlinear ODE of second order with boundary conditions (q2prime), that’s the equation (32) of the
appendix A.4 from the paper. That is an expression of exogeneous parameters, the prime derivative (qprime) and the function itself
(q). The aim goal is to find a density function q(e) (it’s a price, so q>0) with argument "e" (e>0). The boundary conditions
are qprime(e_)=0 and qprime(e^)=0, where e_ is a lower limit (e_>0) and e^ a higher limit. The functions are:

function yinit = Finitwithproduction_rw(x)

% —- This function is the initial guess of the ODE solution, in the form of y(x) —–

global qprime elbarnew qlbarnew eubar qinf

% Note: Finitwithproduction returns a linear interpolation of functions q(e)
% (JCSC), returns [q(e), q'(e)]

qprime = (qinf – qlbarnew)/(eubar-elbarnew);
yinit = [qlbarnew + qprime*x
qprime];
end

function res=Fbcwithproduction_rw(ya, yb);

%(JCSC) Boundary Limits. The conditions
%required by the author (q'(e_) = 0, q'(e^-) = 0)
global qlbarnew

res=[%ya(1)-qlbarnew
ya(1)
yb(1)-0];
end

function dh=Fvecwithproduction_rw(x,y)

%(JCSC, 07-02-2024) Wouters Model
% Entrada x: valor de e; y=[q(e), q'(e)]
% Salida [q'(e), q”(e)]

global k delta kappa A iota m lambda zeta phi rho theta …
sigma alpha_ti alpha_sb alpha_h i w sigma_r sigma_e …
rr qprime q q2prime eta

%(JCSC) Here x is e in the paper. y is a vector y(x) = [q(x), q'(x)]

q=y(1,:);
qprime=y(2,:);
i = (delta + (q-1)./kappa); %(JCSC, 040224)
w = q.*k; %(JCSC ok)
%alpha_h = min(1, ((x.*k)./((1-lambda).*w))); %(JCSC, 040224)
alpha_h = min((1-lambda), ((x.*k)./w)); %(JCSC, 040224)
alpha_ti = ((((1-phi).*alpha_h.*(1-lambda)).^(-1)).*(zeta*(1-(qprime/q).*x).*sigma + …
(qprime./q).*x.*phi)-(phi./(1-phi)))./(zeta.*(1-(qprime./q).*x).*sigma + (1-m).*(qprime./q).*x.*phi); %(JCSC, 270124 ok)
sigma_e = (x.*(((1-phi).*m.*alpha_ti – 1).*sigma + (phi./zeta)))./(1-(1-phi).*m.*alpha_ti.*x.*(qprime./q)); %(JCSC, revisado 040224)
sigma_r = sigma + (qprime./q).*sigma_e;
alpha_sb = 1./(zeta.*sigma_r); %(JCSC,ok)
f = A – delta + ((1 – q.^2)./(2.*kappa)); %(JCSC, nuevo modelo ajustado con la inclusión del empleo)
X1 = (1-phi).*m.*alpha_ti + phi.*alpha_sb; %(JCSC, revisado ok)
FF = (q./x) – X1.*qprime; %(JCSC, revisado ok)
G = kappa.*f.*FF + q.*((phi + m.*(1-phi)) – X1).*theta.*q.*qprime; %(JCSC, revisado, 270124)
rr = rho + theta.*(i-delta)-((0.5.*theta.*qprime.*qprime./(kappa.*f)).*(sigma_e.^2)) – …
0.5.*theta.*(theta+1).*((sigma-(q.*qprime.*sigma_e./(kappa.*f))).^2);
A0 = 0.5.*(sigma_e.^2) + ((qprime.*kappa.*f)./G).*(0.5.*X1.*(sigma_e.^2) – …
0.5.*theta.*(q.*q./kappa.*f).*(sigma.^2).*(m-phi.*m+phi-X1)).*((q./G).*(0.5.*theta.*q.*qprime.*X1.*(sigma_e.^2) + …
0.5.*theta.*q.*FF.*(sigma_e.^2)));
A1 = -q.*(eta+i-delta) + X1.*(-delta.*q+A.*(1-iota)) + q.*(m-m.*phi+phi-X1).*rr;
A2 = kappa.*f.*FF.*rr – theta.*q.*q.*qprime.*(-eta-i+delta) – theta.*q.*qprime.*X1.*(A.*(1-iota) – delta.*q);
q2prime = (1./A0).*(alpha_ti.*(sigma_r.^2)./q) – A.*(1-iota) + delta.*q – qprime.*(kappa.*f./G).*A1 + (q./G).*A2;

dh=[y(2,:)
q2prime];
end

In the paper, the author identified values of "e" between 0.65 and 1.66 but running this range
on my code it doesn’t work, so i had problems either identifiying the initial values or some equation inside
the Fvecwithproduction function, but i checked each one to avoid errors or anything else. The Fvecwithproduction
function just has the different equations required to get the equation (32) in the paper append.

From a numerical point of view, I would like to get some help or recommendation to establish the origin of the error, or some
help on matlab to find or explore initial values in this kind of exercises.

Thanks a lot
Juan Camilo,Hi, i’m trying to replicate the numerical solution of a recent paper "A macroeconomic model with
Heterogeneous and financially constrained intermediaries. Wouters, R (2019)"
(https://www.econstor.eu/bitstream/10419/207747/1/1067667695.pdf) but after many runs and different
modifications i always get the following issue:

Error using bvp4c
Unable to solve the collocation equations — a singular Jacobian encountered.

function yinit = Finitwithproduction_rw(x)

% —- This function is the initial guess of the ODE solution, in the form of y(x) —–

global qprime elbarnew qlbarnew eubar qinf

% Note: Finitwithproduction returns a linear interpolation of functions q(e)
% (JCSC), returns [q(e), q'(e)]

qprime = (qinf – qlbarnew)/(eubar-elbarnew);
yinit = [qlbarnew + qprime*x
qprime];
end

function res=Fbcwithproduction_rw(ya, yb);

%(JCSC) Boundary Limits. The conditions
%required by the author (q'(e_) = 0, q'(e^-) = 0)
global qlbarnew

res=[%ya(1)-qlbarnew
ya(1)
yb(1)-0];
end

function dh=Fvecwithproduction_rw(x,y)

%(JCSC, 07-02-2024) Wouters Model
% Entrada x: valor de e; y=[q(e), q'(e)]
% Salida [q'(e), q”(e)]

global k delta kappa A iota m lambda zeta phi rho theta …
sigma alpha_ti alpha_sb alpha_h i w sigma_r sigma_e …
rr qprime q q2prime eta

%(JCSC) Here x is e in the paper. y is a vector y(x) = [q(x), q'(x)]

dh=[y(2,:)
q2prime];
end

Thanks a lot
Juan Camilo, Hi, i’m trying to replicate the numerical solution of a recent paper "A macroeconomic model with
Heterogeneous and financially constrained intermediaries. Wouters, R (2019)"
(https://www.econstor.eu/bitstream/10419/207747/1/1067667695.pdf) but after many runs and different
modifications i always get the following issue:

Error using bvp4c
Unable to solve the collocation equations — a singular Jacobian encountered.

function yinit = Finitwithproduction_rw(x)

% —- This function is the initial guess of the ODE solution, in the form of y(x) —–

global qprime elbarnew qlbarnew eubar qinf

% Note: Finitwithproduction returns a linear interpolation of functions q(e)
% (JCSC), returns [q(e), q'(e)]

qprime = (qinf – qlbarnew)/(eubar-elbarnew);
yinit = [qlbarnew + qprime*x
qprime];
end

function res=Fbcwithproduction_rw(ya, yb);

%(JCSC) Boundary Limits. The conditions
%required by the author (q'(e_) = 0, q'(e^-) = 0)
global qlbarnew

res=[%ya(1)-qlbarnew
ya(1)
yb(1)-0];
end

function dh=Fvecwithproduction_rw(x,y)

%(JCSC, 07-02-2024) Wouters Model
% Entrada x: valor de e; y=[q(e), q'(e)]
% Salida [q'(e), q”(e)]

global k delta kappa A iota m lambda zeta phi rho theta …
sigma alpha_ti alpha_sb alpha_h i w sigma_r sigma_e …
rr qprime q q2prime eta

%(JCSC) Here x is e in the paper. y is a vector y(x) = [q(x), q'(x)]

dh=[y(2,:)
q2prime];
end

Thanks a lot
Juan Camilo, doe, ordinary differential equations, singularity, initial values ode MATLAB Answers — New Questions

Cart

Cart