I already asked this on stack overflow but was redirected here.

I am trying to solve the (highly non-linear) following system of equations :


and $g_0$, $\Lambda$ are positive constants which i can attribute a numerical value to. $\sinh^{-1}$ is the inverse hyperbolic sine function.

Mathematica gives me no result when I try to use the function Solve. So I thought it would be a good idea to first try to numerically evaluate the integrals $I_n$ in terms of the parameters $\Delta_1$, $\Delta_2$ and $\Delta_3$ before putting it back into the equations.

My question : is there any function designed to give a numerical approximation of an integral which depends on some real parameters ? It could be for instance a polynomial on these parameters, or something more complicated...

Thanks in advance for your time.

EDIT : Here is the (corrected according to comments) code I've tried to use :

f0[d0_?NumericQ, d1_?NumericQ, d2_?NumericQ] := NIntegrate[
ArcSinh[(10*d0)/(d0 + d1*Cos[2*t]^2 + d2*Cos[2*t])], {t, 0, 2*Pi}]
f1[d0_?NumericQ, d1_?NumericQ, d2_?NumericQ] := 
Cos[2*t]*ArcSinh[(10*d0)/(d0 + d1*Cos[2*t]^2 + d2*Cos[2*t])], {t, 0,
f2[d0_?NumericQ, d1_?NumericQ, d2_?NumericQ] := 
Cos[2*t]^2*ArcSinh[(10*d0)/(d0 + d1*Cos[2*t]^2 + d2*Cos[2*t])], {t, 
0, 2*Pi}]
f3[d0_?NumericQ, d1_?NumericQ, d2_?NumericQ] := 
Cos[2*t]^3 ArcSinh[(10*d0)/(d0 + d1*Cos[2*t]^2 + d2*Cos[2*t])], {t, 
0, 2*Pi}]
f4[d0_?NumericQ, d1_?NumericQ, d2_?NumericQ] := 
Cos[2*t]^4*ArcSinh[(10*d0)/(d0 + d1*Cos[2*t]^2 + d2*Cos[2*t])], {t, 
0, 2*Pi}]

FindRoot[{d0 == 
d0*f2[d0, d1, d2] + d1*f4[d0, d1, d2] + d2*f3[d0, d1, d2], 
d1 == d0*f0[d0, d1, d2] + d1*f2[d0, d1, d2] + d2*f1[d0, d1, d2], 
d2 == 2*(d0*f1[d0, d1, d2] + d1*f3[d0, d1, d2] + 
d2*f2[d0, d1, d2])}, {{d0, 1.}, {d1, 1.}, {d2, 1.}}]

I've tried many different guesses for the initial values but it looks like it is always converging to zero. Any suggestion of a more efficient way to approach the problem would be very appreciated. Thanks !

  • $\begingroup$ You will need to use FindRoot[] and NIntegrate[] for this task; see if you can come up with good initial estimates for your roots. $\endgroup$ Jun 8, 2016 at 13:09
  • $\begingroup$ Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – user9660
    Jun 8, 2016 at 13:11
  • $\begingroup$ Please show us the code of what you have tried. Showing the nicely formatted equation is helpful but it is beneficial for people who try to answer your question if you also write the code in Mathematica format. You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. $\endgroup$ Jun 8, 2016 at 13:47
  • $\begingroup$ One correction is the spelling of Nintegrate. Should use NIntegrate. When corrected FindRoot yields an answer but I doubt that is is very good. All the numbers are very close to zero. $\endgroup$ Jun 8, 2016 at 22:15
  • $\begingroup$ I knew it might be something like that.. Thanks ! I'll try to find a correct initial guess for the parameters. $\endgroup$
    – Dimitri
    Jun 9, 2016 at 8:33

1 Answer 1



Realize that there is a difference between your original nicely printed equations and the code that you have typed.

The difference is in the numerator of the function to be integrated.

In the nicely typed set of equations you have Δ in the numerator where as in the code you have typed you have Δ0.


Mathematica graphics


Mathematica graphics

Which is correct?

Original Answer

As I interpret the function inside the integral I see

fun[n_, Δ_, Δ0_, Δ1_,  Δ2_, θ_] := 
 ArcSinh[Δ/( Δ0 + Δ1 Cos[2 θ]^2 + Δ2 Cos[2 θ] )] Cos[2 θ]^n

This function plots nicely when the denominator stays positive

Plot[fun[2, 10, 3, 2, 1, θ], {θ, 0, 2 π}]

Mathematica graphics

but has some nasty singularities when we allow the denominator to go negative

Plot[fun[2, 10, 1, 2, 3, θ], {θ, 0, 2 π}]

Mathematica graphics

Numerical Integration

Numerical integration works find for the nice parameters

NIntegrate[fun[2, 1, 3, 2, 1, θ], {θ, 0, 2 π}]

(* 0.725971 *)

but we get a complaint for the parameters that produce singularities

NIntegrate[fun[2, 1, 1, 2, 3, \[Theta]], {\[Theta], 0, 2 \[Pi]}]

Mathematica graphics

(* -5.6803 *)
  • $\begingroup$ Thanks for pointing that nasty behavior. For now I'm mostly concerned with figuring out the right way to implement the system of equations, I've just edited my question to include the code. $\endgroup$
    – Dimitri
    Jun 8, 2016 at 15:15
  • $\begingroup$ $\Lambda$ is a high-energy cutoff in my equations (which are gap equations for a certain type of unconventional superconductor). Which is why i set $\Lambda \sim 10 \Delta_0$. But now I've realized that it effectively makes $\Lambda$ a variable, and thus makes the problem more complicated than what it should be... I guess that I should first try to find a guess for what order of magnitude I expect for the $\Delta_i$'s then assign an appropriate numerical value to $\Lambda$. $\endgroup$
    – Dimitri
    Jun 10, 2016 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.