Solving system of integral equations - Error FSOLVE cannot continue

Dear all,

I am trying to solve a system of complicated integral equations in MATLAB.

The equations and initial values are attached in the figure.

Indeed, I want to solve Eqs. (35) and (37) with considering Eq. (36) vanishes in vacuum. Of course the next step will be solving all three equations with together.

I am using Fsolve routine in Matlab R2015a to obtaine "m" and "Delta_MF", But some thing is wrong and I get error massage in the command window. I would really appreciate that if someone could help me to solve my equations.

T=50.;
mu=0.;
mustar = 0.;
Lambda = 602.3;
% Lambda is the upper limit of integral.
GsLambda2 = 2.319;
Gs = 0.00000639;
m0 = 5.5;
Gv = 0.5 * Gs;
Gd = 0.75 * Gs;
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);
Deltagap = @(b) integral(@(p) 4.0/pi.^2*Gd*b(2)*((1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))), 0, 602.3) - b(2);
B = fsolve(@(b) [mgap(b), Deltagap(b)], rand(2,1))


Error using  / 
Matrix dimensions must agree.

Error in
@(p)4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)-mustar)/(sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))+(1-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)+mustar)/(sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))+((exp(-(sqrt(p.^2+b(1).^2)+mustar)/T)+1).^-1)-((exp((sqrt(p.^2+b(1).^2)-mustar)/T)+1).^-1))


Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in
@(b)integral(@(p)4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)-mustar)/(sqrt((sqrt(p.^2+b(1).^2)-mustar).^2+(abs(b(2))).^2))+(1-2.*((exp((sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))/T)+1).^-1))*(sqrt(p.^2+b(1).^2)+mustar)/(sqrt((sqrt(p.^2+b(1).^2)+mustar).^2+(abs(b(2))).^2))+((exp(-(sqrt(p.^2+b(1).^2)+mustar)/T)+1).^-1)-((exp((sqrt(p.^2+b(1).^2)-mustar)/T)+1).^-1)),0,602.3)-(b(1)-m0)


Error in @(b)[mgap(b),Deltagap(b)]


Error in fsolve (line 219)
            fuser = feval(funfcn{3},x,varargin{:});

Error in gapequations (line 41)
B = fsolve(@(b) [mgap(b), Deltagap(b)], rand(2,1))

Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Answer 1

[Matt J]

Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Before you use FSOLVE, you must make sure that your objective function code works.

T=50.;
mu=0.;
mustar = 0.;
Lambda = 602.3;
% Lambda is the upper limit of integral.
GsLambda2 = 2.319;
Gs = 0.00000639;
m0 = 5.5;
Gv = 0.5 * Gs;
Gd = 0.75 * Gs;
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);
Deltagap = @(b) integral(@(p) 4.0/pi.^2*Gd*b(2)*((1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))), 0, 602.3) - b(2);

fun=@(b) [mgap(b), Deltagap(b)];

fun(rand(2,1))

Error using /
Matrix dimensions must agree.

Error in solution (line 11)
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in solution (line 11)
mgap = @(b) integral(@(p) 4.0/pi.^2*Gs*b(1)*1.0/(sqrt(p.^2+b(1).^2))*((1.-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) - mustar)/(sqrt ((sqrt(p.^2+b(1).^2) - mustar).^2 + (abs (b(2))).^2))+(1-2.*((exp ((sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))/T) + 1).^-1))*(sqrt(p.^2+b(1).^2) + mustar)/(sqrt ((sqrt(p.^2+b(1).^2) + mustar).^2 + (abs (b(2))).^2))+((exp (-(sqrt(p.^2+b(1).^2) + mustar)/T) + 1).^-1)-((exp ((sqrt(p.^2+b(1).^2) - mustar)/T) + 1).^-1)), 0, 602.3) - (b(1)-m0);

Error in solution (line 14)
fun=@(b) [mgap(b), Deltagap(b)];

[Mahboubeh Shahrbaf]

Uh thank you Matt,

but the objective functions are introduced in this very popular paper (Eqs. 35, 36, 37)

and it seems lots of people in our field are using these functions.

Now I need to reproduce their results!

I am confused what is wrong!

[Matt J]

I'm happy to assume that the paper is popular for good reason. My point above though, is that you clearly haven't implemented what is in the paper accurately, because your objective function mgap() does not run. It just spits out a flurry of error messages. If I cannot run your objective function, then you cannot expect that fsolve() will have better luck. You need to test mgap() first and make sure it works before giving it to fsolve().

[Mahboubeh Shahrbaf]

Thanks Matt.

I will double check with the implementation of the function in code.

Hopefully I will find the mistake.

But in any case, if the function is true, should the same routine which I had used work? Even if the the upper limit of the integral or the result of it include one or two of unknown variables?

[Matt J]

For fsolve to be successful, the function needs to be differentiable (I can't tell if yours is because of the sqrt's). Also, fsolve could fail if your initial guess is not sufficiently close to a solution. The more accurately you can guess the solution, the better the chances of success.