Numerically compute second integral (not two dimensional integral)

[richard.physics]

Dear Experts,

I would like to obtain numerically a second/double integral with ROOT but not sure how to go about this. To elaborate, suppose I had a TF1 or double-valued function
f(x) = x
and wanted the second integral over the domain [0,1]. The result would be
Int dx’ [Int dx f(x)] = 1/6

I do not want an analytical result, only a number in the end. Since what I want implicitly involves indefinite integration, I suspect that it might not be possible, but maybe I am just not seeing something obvious that involves ROOT::Math::Functor1D and ROOT:::Math::Integrator

Any thoughts or suggestions would be deeply appreciated!

Best,
Richard

---

Please read tips for efficient and successful posting and posting code

_ROOT Version:6.22/06
_Platform: Ubuntu 20.04 (LTS)

---

[Wile_E_Coyote]

ROOT cannot calculate antiderivatives.

[richard.physics]

Hi @Wile_E_Coyote,

Of course ROOT cannot (at least to my knowledge) symbolic algebra/calculus. However, I am not asking for indefinite integrals. As I wrote in the post, I am asking how to perform a definite integral numerically.

One (tedious) way of doing what I want is computing TF1->Integral(0,x) vs x for x = [0,1], and then getting the area under this curve. However, I am hoping that this is a smarter way of doing this.

cheers

[Wile_E_Coyote]

A straight way:

{
  TF1 *f = new TF1("f", [=](double *x, double *p){ return x[0]; }, 0., 1., 0);
  // new TCanvas("c_f", "c_f"); f->Draw();
  TF1 *F = new TF1("F", [=](double *x, double *p){ return f->Integral(p[0], x[0]); }, 0., 1., 1);
  // F->SetParName(0, "lower limit of integration");
  F->SetParameter(0, 0.); // lower limit of integration
  std::cout << F->Integral(F->GetParameter(0), 1.) << std::endl;
  // new TCanvas("c_F", "c_F"); F->Draw();
}

Actually, it seems to me that one can use the Cauchy formula for repeated integration:

{
  TF1 *f = new TF1("f", [=](double *x, double *p){ return x[0]; }, 0., 1., 0);
  // new TCanvas("c_f", "c_f"); f->Draw();
  TF1 *F = new TF1("F", [=](double *x, double *p){ return (p[0] - x[0]) * f->Eval(x[0]); }, 0., 1., 1);
  // F->SetParName(0, "upper limit of integration");
  F->SetParameter(0, 1.); // upper limit of integration
  std::cout << F->Integral(0., F->GetParameter(0)) << std::endl;
  // new TCanvas("c_F", "c_F"); F->Draw();
}

or simply:

{
  TF1 *F = new TF1("F", [=](double *x, double *p){ return (p[0] - x[0]) * x[0]; }, 0., 1., 1);
  // F->SetParName(0, "upper limit of integration");
  F->SetParameter(0, 1.); // upper limit of integration
  std::cout << F->Integral(0., F->GetParameter(0)) << std::endl;
  // new TCanvas("c_F", "c_F"); F->Draw();
}

[richard.physics]

Hi @Wile_E_Coyote,

Thank you for the suggestions. These are great and gives me what I want. Though I should have realized, I did not know that one could define with an inline function, i.e., [=](double *x, double *p){ return ...; },

Thank you for the help!