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

5 views (last 30 days)
Juan Camilo Santana Contreras
Moved: Torsten on 19 Jun 2024
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,
  10 Comments
Show 8 older comments
Torsten
Torsten on 16 Feb 2024
I would need to add some conditions more in the searching to get a smoother path: q(e^) = qinf and q(e_)=tau, but i don't know how incorporate it.
You have a second-order differential equation. So there is no way to incorporate more than two boundary conditions.
Juan Camilo Santana Contreras
Juan Camilo Santana Contreras on 19 Jun 2024
Moved: Torsten on 19 Jun 2024
Hi @Torsten, i was testing my code to get some solutions for a second order DE of a price q(x), but i'm getting some errors or unexpected jumps on the smooth solution (see attached file). At the smallest values of the domain, the smooth is not good, as you can see in the first and second order derivatives. The first 1000 points of the domain are ilustrated. Do you know how could i check the problem to get a better solution, a smoother solution. Thanks for any advice.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Matrix Indexing in Help Center and File Exchange

Tags

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!