1
$\begingroup$

I'm struggling with the following code. I have to solve some self consistent equations with FindRoot, but I always get some problem related to slow convergence and/or problems related with NIntegrate.

ClearAll["Global`*"];

U = 0.5;
S = 5;
w1 = 0;
a = 0.1;
R = 0;
a1 = a/S;
b = 2;
b1 = b*a/S^2;

ϕ[z_?NumericQ] := (1 + Erf[z/Sqrt[2]])/2;
x[q_?NumericQ,m_?NumericQ] := x[q,m] = Sqrt[(1 - a1)/b1]*(m/(Sqrt[q] + a1*R));
y[q_?NumericQ,m_?NumericQ] := y[q,m] = Sqrt[(1 - a1)/b1]*(U - w1 (S - 1)/(2*S) - b1*R/(2*Sqrt[q]))*(1/(Sqrt[q] + a1*R));


fq1[q_?NumericQ, m_?NumericQ, w_?NumericQ] := fq1[q,m,w] = NIntegrate[Exp[-z^2/2]/Sqrt[2 π]*(ϕ[z])^S, {z, y[q,m]*a1 + x[q,m] - I*Sqrt[a1]*w, ∞}];
fq2[q_?NumericQ, m_?NumericQ, w_?NumericQ] := fq2[q,m,w] = NIntegrate[Exp[-z^2/2]/Sqrt[2 π]*(ϕ[z + y[q,m]])^S, {z, -y[q,m]*(1 - a1) + x[q,m] - I*Sqrt[a1]*w, ∞}];
fq3[q_?NumericQ, m_?NumericQ, w_?NumericQ] := fq3[q,m,w] = NIntegrate[Exp[-z^2/2]/Sqrt[2 π]*ϕ[z - y[q,m]]*(ϕ[z])^(S - 1), {z, y[q,m]*a1 + x[q,m]-I*Sqrt[a1]*w, ∞}];
fq[q_?NumericQ,m_?NumericQ] := fq[q,m] = (1 - a)/a1*NIntegrate[Exp[-w^2/2]/Sqrt[2 π]*fq1[q,m,w], {w, -∞, ∞}] + NIntegrate[Exp[-w^2/2]/Sqrt[2 π]*fq2[q,m,w], {w, -∞, ∞}] + (S - 1)*NIntegrate[Exp[-w^2/2]/Sqrt[2 π]*fq3[q,m,w], {w, -∞, ∞}];

fm[q_?NumericQ, m_?NumericQ] := fm[q, m] = 1/(1 - a1)*NIntegrate[Exp[-w^2/2]/Sqrt[2 π]*fq2[q, m, w], {w, -∞, ∞}] - fq[q, m]*a1 /(1 - a1);

FindRoot[{q - fq[q, m] == 0, m - fm[q, m] == 0}, {{q, 1}, {m, 1}}]

In particular it gives the errors

Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

The integral failed to converge after 1012 integrand evaluations. NIntegrate

Any hints and suggestions are welcome.

$\endgroup$
  • 1
    $\begingroup$ You define x[] and y[] to be functions of q and m, but then you use x and y apparently as simple variables that aren't functions. If you replace each of your "bare" x and y with x[q,m] and y[q,m] then that might be a first start. It is difficult to guess whether that is the real problem or not. $\endgroup$ – Bill Mar 5 '17 at 22:27
  • 2
    $\begingroup$ Figure out which integral fails to converge and work on that one. If you can't make it work, update your question with the specific integral. Block[{res}, Check[res = NIntegrate[...], Print["It's me: <include function name, inputs>"]; res]] is one way to catch the troublemaker. $\endgroup$ – Michael E2 Mar 5 '17 at 22:45
  • $\begingroup$ @Bill, you are right. In my code actually x and y are x[q,m] and y[q,m], I didn't write them correctly here. $\endgroup$ – Michelangelo Mar 6 '17 at 8:42
  • $\begingroup$ @MichaelE2 thank you, I will try to see which integral fails to converge. $\endgroup$ – Michelangelo Mar 6 '17 at 8:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.