0
$\begingroup$

Problem 1
My code is given as follows

Clear["Global`*"]
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
int1[s_?NumericQ, a_] := 
 NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(
  k^3 (k^2 - 1)), {k, ((d + s^2)/(1 + d))^(1/2), s}]
int2[x_?NumericQ, a_] := NIntegrate[int1[s, a], {s, 1, x}]
Plot3D[int2[x, a] - (1/2) p*(x^2 - 1) + (1/4) x^2*
    Log[(d + x^2)/x^2] y^2, {x, 0.5, 9.5}, {y, -1.5, 1.5}, 
  PlotRange -> {-0.3, 0.15}] // AbsoluteTiming

It runs about 209s in my computer. How to improve the speed? Any suggestion is much appreciated!
Problem 2
When we vary a, the code is

Clear["Global`*"]
d = 0.8;
p = 0.26;
b = 0.03;

ContourPlot3D[
 NIntegrate[
    NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(
     k^3 (k^2 - 1)), {k, ((d + s^2)/(1 + d))^(1/2), s}], {s, 1, x}] - 
   1/2 p*(x^2 - 1) + 1/4 x^2*Log[(d + x^2)/x^2] y^2 == 0.029744, {x, 
  0.5, 13}, {y, -1.5, 1.5}, {a, 1.95, 2.08}]

For this problem, I can't get a result. How can we deal with this problem?

$\endgroup$
1
  • 1
    $\begingroup$ It takes so much time, because int1 does integration for overlapping k ranges multiple times for different s. If you want to do it with NIntegrate, form a double integral like doubleint[x_, a_] := NIntegrate[(2 (k^4 - 1) (1 + a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)), {s, 1, x}, {k, ((d + s^2)/(1 + d))^(1/2), s}] . Attention, the integration which depends on s must be more outside. This is 20 times faster. $\endgroup$
    – Akku14
    Oct 8 '20 at 13:04
2
$\begingroup$

If you want to do or you must do integration numericaly (e.g. a is not an integer), you can get antiderivatives as interpolation functions with NDSolve.

Clear["Global`*"];
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
f[k_] = (2 (k^4 - 1) (1 + 
     a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)
) //   Simplify;

Get the first antiderivative i1sol. The intial conditions are not decisive, only differences count as shown in the answer of @cvgmt.

(i1sol = 
  i1 /. First@
 NDSolve[{i1'[k] == f[k], i1[1/2] == 0}, i1, {k, 1/2, 10}];

i2sol = 
  i2 /. First@
 NDSolve[{i2'[s] == i1sol[s] - i1sol[((d + s^2)/(1 + d))^(1/2)], 
   i2[1] == 0}, i2, {s, 1/2, 10}];

Plot3D[((i2sol[x] - i2sol[1]) - (1/2) p*(x^2 - 1) + (1/4) x^2*
  Log[(d + x^2)/x^2] y^2), {x, 0.5, 9.5}, {y, -1.5, 1.5}, 
   PlotRange -> {-0.3, 0.15}, ImageSize -> 300]
) // Timing

This is even faster than the analytical approach.

Edit to Problem 2

In higher versions of MMA you would do this, when using the above method, with ParametricNDSolve with a as parameter. Since i use version 8.0, i have to do it this way:

d = 0.8;
p = 0.26;
b = 0.03;
f[k_, a_] = (2 (k^4 - 1) (1 + 
   a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)) // Simplify

i1sol = i1 /. 
  First@NDSolve[{Derivative[1, 0][i1][k, a] == f[k, a], 
 i1[1/2, a] == 0}, i1, {k, 1/3, 13}, {a, 1.95, 2.08}]

i2sol = i2 /. 
  First@NDSolve[{Derivative[1, 0][i2][s, a] == 
  i1sol[s, a] - i1sol[((d + s^2)/(1 + d))^(1/2), a], 
 i2[1, a] == 0}, i2, {s, 1/2, 13}, {a, 1.95, 2.08}]

ContourPlot3D[(i2sol[x, a] - i2sol[1, a]) - 1/2 p*(x^2 - 1) + 
   1/4 x^2*Log[(d + x^2)/x^2] y^2 == 0.029744, {x, 0.5, 13}, {y, -1.5,
   1.5}, {a, 1.95, 2.08}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thanks, It is a smart method! And do you have any ideas of a similar for Problem 2 I added in my question when a is not fixed? $\endgroup$
    – keanhy14
    Oct 9 '20 at 0:36
1
$\begingroup$

Use Newton-Leibniz Formula $\int_a^b f(t)\,dt=F(b)-F(a)$

We use F[x_]=Integrate[f[x],x] and subtract F[b] and F[a] instead of Integrate[f[x],{x,a,b} to calculate $\int_a^b f(t)\,\mathrm{d}t$ since Integrate[f[x],x] is easy to calculate.

Clear["Global`*"]
d = 0.8;
p = 0.26;
a = 2;
b = 0.03;
intA[k_] = 
 Integrate[(2 (k^4 - 1) (1 + 
       a b (k^2 + k^-2 - 2)^(a - 1)))/(k^3 (k^2 - 1)), k]
int1[s_, a_] = intA[s] - intA[((d + s^2)/(1 + d))^(1/2)];
intB[s_] = Integrate[int1[s, a], s];
int2[x_, a_] = intB[x] - intB[1];
Plot3D[int2[x, a] - (1/2) p*(x^2 - 1) + (1/4) x^2*
    Log[(d + x^2)/x^2] y^2, {x, 0.5, 9.5}, {y, -1.5, 1.5}, 
  PlotRange -> {-0.3, 0.15}] // AbsoluteTiming

enter image description here

$\endgroup$
2
  • $\begingroup$ It is a useful method. But I wonder if there are other ways to involve the 'NIntegrate', since that when I evaluate any points in the figure, for instance, 'int2[1.5, a] // AbsoluteTiming', it gives the result very fast. $\endgroup$
    – keanhy14
    Oct 8 '20 at 12:14
  • 1
    $\begingroup$ This method seems not to work if a is not fixed. The code is given in Problem2, do you have any suggestions? Thanks! $\endgroup$
    – keanhy14
    Oct 9 '20 at 0:33

Your Answer

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

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