1
$\begingroup$

I have the following code to calculate a numerical integral for any given a, however it takes a very long time, even with adaptivemontecarlo, which is not accurate enough:

 Needs["NumericalCalculus`"]

 a=15

 g[t_] := {-(a + 2*Cos[2*t])*Sin[3*t], (a + 2*Cos[2*t])*Cos[3*t], 2*Sin[2*t]}

 dg[t_] := If[t - 2*Pi <= 0, g'[t], g'[2*Pi]];

 tfn :=
    NDSolveValue[
      {t'[s] == 1/Norm[dg[t[s]]], t[0] == 0, 
        WhenEvent[t[s] == 2*Pi, "StopIntegration"]}, 
      t, {s, 0, 2*Pi + NIntegrate[Norm[g'[t]], {t, 0, 2*Pi}]}];

 l := NIntegrate[Norm[D[g[t], t]], {t, 0, 2*Pi}]

 c1[s_] := g[tfn[s]]

 j[s_] := Normalize[Cross[D[c1[s], s], D[D[c1[s], s], s]]]
 v[s_] := Normalize[Cross[j[s], Normalize[D[c1[s], s]]]]

 T2[s_, y_] := c1[s] + 0.4*(j[s]*Sin[y] - v[s]*Cos [y])

 s := p
 y := 500*Pi*p/l

 der[x_?NumericQ] := (Norm[ND[T2[s, y], {p, 2}, x]])^2

 NIntegrate[der[x], {x, 0, l}, Method -> "AdaptiveMonteCarlo"]/l

How could I improve the efficiency of this code?

Edit

I forgot to mention that the integral will be evaluated for many values of "a", and that's why I used a lot of := signs.

Disclaimer: I have opened a new thread since my question is different this time.

Improve this question
$\endgroup$
6
  • $\begingroup$ you can actually remove all of the delayed definitions except for der. Change the s in the NDSolveValue expression to something else, so it doesnt get confused by your later unrleated assignment of s for some other purpose. $\endgroup$
    – george2079
    Commented May 22, 2015 at 16:11
  • $\begingroup$ OK, I am trying your suggestions, but it takes incredibly long with = signs. The integration was taking 56 seconds with maxrecursion 10 and 23 with maxrecursion 1, where both gave the same result to a sufficient accuracy. But I will calculate with 50 different values of a. $\endgroup$
    – Baran Cimen
    Commented May 22, 2015 at 16:20
  • $\begingroup$ ND[T2[s, y] is the numerical derivative of T2. tfn is an interpolation function, so it is a function of s. $\endgroup$
    – Baran Cimen
    Commented May 22, 2015 at 16:41
  • $\begingroup$ george2079, the s is the variable I used for arclength. It is same for all expressions. $\endgroup$
    – Baran Cimen
    Commented May 22, 2015 at 16:50
  • $\begingroup$ Wrap NDSolveValue in Module : Module[{s}, NDSolveValue[] ] if you really must re-use the same symbol as the ode independent variable as you use elsewhere. $\endgroup$
    – george2079
    Commented May 22, 2015 at 18:00

1 Answer 1

Reset to default
2
$\begingroup$

This is a stab at cleaning it up. I put in table form so you can see how to loop over a.

 Needs["NumericalCalculus`"]
 Table[
   g[t_] = {-(a + 2*Cos[2*t])*Sin[3*t], (a + 2*Cos[2*t])*Cos[3*t], 
       2*Sin[2*t]};
   dg[t_] = If[t - 2*Pi <= 0, g'[t], g'[2*Pi]];
   tfn = Module[{s}, NDSolveValue[{t'[s] == 1/Norm[dg[t[s]]], t[0] == 0,
       WhenEvent[t[s] == 2*Pi, "StopIntegration"]}, t,
          {s, 0, 2*Pi + NIntegrate[Norm[g'[t]], {t, 0, 2*Pi}]}]];
   el = NIntegrate[Norm[D[g[t], t]], {t, 0, 2*Pi}];
   c1[s_] = g[tfn[s]];
   j[s_] = Normalize[Cross[D[c1[s], s], D[D[c1[s], s], s]]];
   v[s_] = Normalize[Cross[j[s], Normalize[D[c1[s], s]]]];
   T2[s_, y_] = c1[s] + 0.4*(j[s]*Sin[y] - v[s]*Cos[y]);
   y = 500*Pi*s/el;
   der[x_?NumericQ] := (Norm[ND[T2[s, y], {s, 2}, x]])^2;
    {a, NIntegrate[der[x], {x, 0, el}, MaxRecursion -> 0]/el}  ,
      {a, 10, 10} ]

NOTE for testing purpose I put MaxRecursion -> 0 so that NIntegrate finishes in some reasonable time (the result is highly inaccurate ).

This is your integrand:

enter image description here

You can increase MaxRecursion, but i doubt you will get any sort of accuracy.

Improve this answer
$\endgroup$
0

Your Answer

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

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