Singular Jacobian in bvp4c while Mathematica has no problem

Hi,

I am trying to solve the following boundary value problem:

,

with boundary conditions at and at . The parameters are all positive real numbers and .

Sadly, Matlab gives me the classic error: Unable to solve the collocation equations -- a singular Jacobian encountered. I paste my code below.

Note that there is no particular reason this should be ill behaved for , the system has an analytical solution which is regular over the mentioned interval. Mathematica's boundary value problem solver has no problem:

Could you please help me figure out what is going wrong in my Matlab code?

A=0.1;
B=7;
mu=1;
gammaR=2;
gammaTheta=5;

solinit = bvpinit( linspace(A,B,100),@(x)guess(x, gammaR, gammaTheta, B, mu));
sol = bvp4c(@(x,y)bvpfcn(x,y, gammaR, gammaTheta, mu), @(ya, yb)bcfcn(ya,yb, A), solinit);

--------------

function dydx = bvpfcn(x,y, gammaR, gammaTheta, mu)


dydx = zeros(2,1);
dydx = [gammaR*gammaTheta*x/y(1)
       2*mu*gammaR/(x*gammaTheta)*(1-x^4*gammaTheta^4/y(1)^4 )];
end

--------------

function res = bcfcn(ya,yb, A)

res = [ya(1)-A
       yb(1)];
end

--------------

function g = guess(x, gammaR, gammaTheta, B, mu)
g = [mu*(gammaR/gammaTheta - gammaTheta / gammaR)*log(x/B)
    x*sqrt(gammaR*gammaTheta)];
end

[Star Strider]

The problem occurs in ‘bvpfcn’ when x=6.965151515151 .... so the first time this occurs, dydx=[NaN;NaN] and in the next iteration, ‘dydx(1)’ is 0 and the following iteration throws the error.

None of the approaches I attempted to fix it worked (one of which was defining ‘y(1)’ to be the maximum of ‘y(1)’ and 1E-12 ) so the fix apparently depends on knowing more about how bvp4c works.

I tweaked ‘bvpfcn’ to be —

function dydx = bvpfcn(x,y, gammaR, gammaTheta, mu)

x

dydx = zeros(2,1);

dydx = [gammaR*gammaTheta*x/y(1)

2*mu*gammaR/(x*gammaTheta).*(1-x.^4*gammaTheta^4/y(1)^4 )]

fprintf('\n----------------------------\n')

end

in order to see if that would help define the problem.

It defines it, however offers no insight as to what is actually throwing the error.

I’m not posting this as an Answer because it didn’t solve the problem.

.

[Alexander Erlich]

Thank you! My friend tested this same bvp in Julia which ran successfully in less than 2s. It appears to be a Matlab issue. I would be curious to know how to solve this in Matlab.

[Star Strider]

My pleasure!

Since no other person has posted anything, Contact Support since this may well be a bug that hasn’t previoously been discovered. Include the URL of this thread in the message so that you don’t have to repeat everything.

Please post back here with what MathWorks suggests, since others may have the same problem and would be searching for it on Answers. (I would also be interested in knowing what the problem is and what the fix for it is.)

EDIT — (28 Dec 2021 at 2:08)

I submitted a bug report to MathWorks on this. The reply I got last week was not helpful, so I am not posting it.

.

[Torsten]

function res = bcfcn(ya,yb, A)
res = [ya(1)-A
       yb(2)];
         ...
end