MATLAB: Different results for solving integral equation in each run

fsolveintegral equation

Dear all,

I am trying to solve the following integral equation:

T=100.;

mustar = 400.;

Lambda = 602.3;

GsLambda2 = 2.319;

Gs = 0.00000639;

m = 100.;

Gv = 0.5 * Gs;

Gd = 0.75 * Gs;

Deltagap = @(d) integral(@(p) 4.0/pi.^2*Gd*d.*p.^2.*(((1.-2.*((exp (( sqrt (( sqrt (p.^2.+m.^2.) – mustar ).^2 + …

d.^2))/T) +1).^(-1)))./(sqrt (( sqrt (p.^2.+m.^2.) – mustar).^2 + d.^2))) + …

((1. – 2.*((exp (( sqrt (( sqrt (p.^2 + m.^2) + mustar).^2 + d.^2))/T)…

+ 1).^(-1)))./(sqrt (( sqrt (p.^2+m.^2) + mustar).^2 + d.^2)))), 0, 602.3) – d;

B = fsolve(@(d) Deltagap(d), rand(1,1))

Each time when I run the code, I got a different answer.

I don't know which one is correct.

What should I do?

I have to run the code for the different value of the mustar as well as T thus I must manage to get correct result for every run.

Thank you in advance.

Best Answer

the fsolve funciton is a robust root-finding algorithm. If you plot ‘Deltagap’, you will find that it has no zero-crossings at all for positive values of ‘d’, and for ‘d’ itself crossing 0 (so allowing for negative and positive values), it has exactly 1 zero-crossing, that being exactly at 0.

Related Solutions

MATLAB: Solving system of integral equations – Error FSOLVE cannot continue

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)];

MATLAB: How to solve this equation which contains the complementary error function

As Walter has pointed out, it is surely not possible for 'solve' to find an explicit solution for t in terms of the other variables, since your expression involves the 'erfc' function. The best you can hope for is a numerical solution. To accomplish this I would suggest you make the following substitutions in G to simplify matters. Let x = sqrt(d*t)/h and k = 4*h^2/a^2 so that

 G = g(x) = ((1-2*x^2)*exp(x^2)*erfc(x)+2*x/sqrt(pi))/(1+k*x^2)

Then

 dG/dd = dg/dx*dx/dd = g'(x)*sqrt(t/d)/(2*h)

and

 d(dG/dd)/dt = g"(x)*dx/dt*sqrt(t/d)/(2*h)+g'(x)/sqrt(d*t)/(4*h)
   = g"(x)*(sqrt(d/t)/(2*h))*sqrt(t/d)/(2*h)+g'(x)/sqrt(d*t)/(4*h)
   = g"(x)/(4*h^2)+g'(x)/(4*h^2)/x 
   = (g"(x)*x+g'(x))/(4*h^2*x)

Hence the equation you want to solve is

 g"(x)*x+g'(x) = 0

For this purpose you can use matlab's 'fzero'. For any given value of k=4*h^2/a^2 you can find the x root or roots of this equation. Having found it or them, you can immediately use the solution(s) to solve for t as

 t = h^2*x^2/d

You can use your symbolic toolbox to evaluate and simplify g'(x) and g"(x) in setting up the expression g"(x)*x+g'(x) for 'fzero' to solve.

Best Answer

Related Solutions

MATLAB: Solving system of integral equations – Error FSOLVE cannot continue

MATLAB: How to solve this equation which contains the complementary error function

Related Question