I face an error when apply bvp4c: Unable to solve the collocation equations -- a singular Jacobian encountered.

The function that I want to solve has been put into a first-order ODE as follows:

function wantq= bvp_Q(p,q)
t=0;
T=30;
n=35;
m=0.4;
m_c=0.7;
kappa=0.05;
hatrho=0.5;
k_1=0.4;
k_2=0.8;
r=0.02;  
mu=0.07;   
sigma=0.2;   
tildesigma=0.4;  
delta=0.03;
lambda_1=0.03;
lambda_2=0.03;

q1=exp(r*(T-t))+m*lambda(n,t)/m_c;
q2=delta*p*g2(t,T,n,m,m_c,kappa);
q3=(delta^2*p^2*(1-p)^2)/(2*tildesigma^2);
q4=lambda_2-(lambda_1+lambda_2)*p-(m+k_2)*delta*p*(1-p)*g2(t,T,n,m,m_c,kappa)...
    -(m*hatrho*delta*p*(1-p)*(mu-r-hatrho*sigma*tildesigma*m*g2(t,T,n,m,m_c,kappa)))/(sigma*tildesigma*(m+k_1));
q5=(delta^2*p^2*(1-p)^2)/(2*tildesigma^2)*(m+k_2-(m^2*hatrho^2)/(m+k_1));

wantq=[q(2); (q1*q(1)-q2-q4*q(2)+q5*(q(2))^2)/q3];
end

The boundary condition is constructed as the following code

function bc= bvp_B(Q1,Q2)
r=0.02;    
delta=0.03;
lambda_1=0.03;
lambda_2=0.03;
t=0;
T=30;
n=35;
m=0.4;
m_c=0.7;
kappa=0.05;

bc=[(exp(r*(T-t))+m*lambda(n,t)/m_c)*Q1(1)-lambda_2*Q1(2);
    (exp(r*(T-t))+m*lambda(n,t)/m_c)*Q2(1)-delta*g2(t,T,n,m,m_c,kappa)+lambda_1*Q2(2)];
end

Then I use the above two functions as follows:

init=bvpinit(linspace(0,1,10),[0 0]);
sol=bvp4c(@bvp_Q,@bvp_B,init);
p=linspace(0,1,100);
BS=deval(sol,p);
plot(p,BS(1,:))

Then it was saying "Unable to solve the collocation equations -- a singular Jacobian encountered."

[Stephan]

What is lambda - it doesnt appear here

[Ning Wang]

Thanks for your answer. lambda and g_2 are two functions not depending on p.

[Stephan]

We need it also to try some stuff... else we get errors...

[Ning Wang]

This is function lambda

function
ction lambda = lambda(n,t)
lambda=1/9.5*exp((n+t-86.3)/9.5);
end

The function g_2 is given as follows

function g2 = g2(t,T,n,m,m_c,kappa)
r=0.02;  % risk-free interest rate
syms s;
syms u;
b1=exp(r*(T-u))+m*lambda(n,u)/m_c+kappa;
b2=exp(r*(T-s))-int(b1,u,t,s);
g2=double(int(b2,s,t,T));
end

Answer 1

[Stephan]

It appears to work by making 2 small changes in the code:

init=bvpinit(linspace(0.001,0.9999,10),[0 0]);
sol=bvp5c(@bvp_Q,@bvp_B,init);
p=linspace(0.001,0.9999,100);
BS=deval(sol,p);
plot(p,BS(1,:))

1. I used bvp5c instead of bvp4c
2. I changed the interval from [0 1] to [0.001,0.9999]

The code is running on my machine since a while without erros - it appears to be a compute-intensive problem. Please let me know if you get a result or an error by trying this code. In this context you should also have a read here:

https://de.mathworks.com/matlabcentral/answers/120468-why-do-i-receive-an-error-about-a-singular-jacobian-when-i-use-the-bvp4c-function-within-matlab

[Ning Wang]

Hi, Stephan

Thanks for your answer and sharing the reading material with me. I used your code and it was running for 40 minutes but I did not get the results yet. Do you have some ideas to speed up the computation? Because my next step is going to have a look at how the solution changes with respect to some interesting parameters, for example, the parameter m. So I am worried I would not get the results easily for changing m if it goes so slow for a fixed m.

[Stephan]

I think the usage of symbolic calculations should be avoided. Symbolic calculations take very much time. You should try solve this by using numeric methods.

[Ning Wang]

yes, that is true. But bvp4c and bvp5c are numerical methods, right?

[Stephan]

Yes, the part of your Code involving syms and with it symbolic variables should get your attention.

[Ning Wang]

OK, I see. I use that to calculate the integrations in function g2. But g2 is indenpent of the function that I want to solve through bvp4c.