The functional forms that we saw in Gauss Quadrature are restrictive. Like, for example the Gauss-Legendre formula is [-1,1]. Ofte in times, we want more flexibility in the intervals. The other orms integrate from and/or till infinity, which makes it not very useful to apply to simpler quadrature rules like the trapezoidal or Simpson’s. Moreover, many of the integrals are difficult to calculate as they are very complicated. There are also many functions taht cannot be evaluated at certain points due to singularity. Integration by substitution writes teh integral that needs to be evaluated, in a simpler form so that it is easier to compute such as by avoiding the singularity. It can also transform the integral to a desired range and also to avoid infinity. Moreover, it speeds up the computation by reducing the number of iterations needed for convergence. The substitution rule ays that we find a change of variable, 

 x=\phi(t)

such that,

 dx=\phi'(t)dt  

We can then replace the integral in  t with an integral in  x .

 \displaystyle \int_{a}^{b}f(\phi(t))\phi'(t)dx=\int_{\phi(a)}^{\phi(b)}f(x)du

Consider for example, suppose we want to calculate,

 \displaystyle I=\int_{0}^{2}t\,cos(t^2+1)dt

If we set ,  x=t^2+1 in the above equation, we have

 dx=2tdt

Or, more specifically

 \displaystyle \frac{1}{2}x=tdt

Now, substituting these in the original integral and eliminating  t , we have

 \displaystyle I=\int_{0}^{2}t\,cos(t^2+1)dx=\int_{x=1}^{x=5}\frac{1}{2}cosxdx

The new integral in  x is easier to compute than the original one.

 \displaystyle I=\frac{1}{2}\int_{1}^{5}cosxdx=\frac{1}{2}(sin5-sin1)

Standard Interval

This particular transformation maps integral interval from  [a,b] to [-1,1]. The substitution rule for the same is 

 \displaystyle x(t)=\frac{(b-a)t+(a+b)}{2}

 \displaystyle x'(t)=\frac{(b-a)}{2}

 \displaystyle \int_{a}^{b}f(x)dx=\int_{-1}^{1}\frac{(b-a)}{2}f\Bigg(\frac{(b-a)t+(a+b)}{2}\Bigg)dt

 

				
					val f: AbstractRealUnivariateFunction = new AbstractUnivariateFunction();
val I: Integrator = new Integrator();

double a = 0, b = 10; // the limits
Integrator integrator1 = new NewtonCotes(3, NewtonCotes.Type.OPEN, 1e-8, 10);
Integrator integrator2 = new ChangeofVariables(new StandardInterval(a, b), integrator1);

double I = integrator2.integrate(
       new AbstractRealUnivariateFunctionction() {
     override public double evaluate (double t) {
      return t; // the original integrand
     }
    },
       a, b // the original limits
    );
    System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[0,10] f(x) dx = 50.000000
				
			

Inverting Variable

This particular transformation is useful when there are three particular conditions, viz.

  1.  b\longrightarrow\infty,a>0
  2.  a\longrightarrow-\infty,b<0
  3. Any function that decreases towards infinity ster than  \displaystyle \frac{1}{x^2}
The integrator for this type of substitution should use an open formula to avoid computing for the end point where  t=0 . The substitution is
 \displaystyle x(t)=\frac{1}{t}
 \displaystyle x'(t)=-\frac{1}{t^2}
 \displaystyle \int_{a}^{b}f(x)x=\int_{1/b}^{1/a}\frac{1}{t^2}f(t)dt,\, ab>0
				
					val f: AbstractRealUnivariateFunction = new AbstractUnivariateFunction();
val I: Integrator = new Integrator();

double a = 1, b = Double.POSITIVE_INFINITY; // the limits
NewtonCotes integrator1 = new NewtonCotes(3, NewtonCotes.Type.OPEN, 1e-15, 10);
ChangeOfVariable integrator2 = new ChangeofVariable(new InvertingVariable(a, b), integrator1);

double I = integrator2.integrate( // I = 1 
       new AbstractRealUnivariateFunctionction() {
     override public double evaluate (double x) {
      return 1 / x / x; // the original integrand
     }
    },
       a, b // the original limits
    );
    System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[1,Infinity] f(x) dx = 1.000000
				
			

Exponential

This transformation is good for when the lower limit is finite, the upper limit is infinite and the integral falls off exponentially. The integrator for this substiution should use an open formula to avoid computing for the end point where  t=0 . The substitution for this method is
 x(t)=-log\,t
 \displaystyle x'(t)=-\frac{1}{t}
 \displaystyle \int_{a}^{\infty}f(x)dx=\int_{0}^{e^a}\frac{f(-log\,t)}{t}dt
				
					val f: AbstractRealUnivariateFunction = new AbstractUnivariateFunction();
val I: Integrator = new Integrator();

double a = 1, b = Double.POSITIVE_INFINITY; // the limits
NewtonCotes integrator1 = new NewtonCotes(3, NewtonCotes.Type.OPEN, 1e-15, 10);
ChangeOfVariable integrator2 = new ChangeofVariable(new Exponential(a), integrator1);

double I = integrator2.integrate( // I = sqrt(PI)/2  
       new AbstractRealUnivariateFunctionction() {
     override public double evaluate (double x) {
      return sqrt(x) * exp(-x); // the original integrand
     }
    },
       a, b // the original limits
    );
System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[0,Infinity] f(x) dx = 0.886227
				
			

Mixed Rule

The mixed rule is good for functions that fall off quckily at infinity, for example,  e^x or  e^{x^2} . The interval is  [0,\infty] . Figuring out a good range for  t is the tricky part of using this particular transformation. The substitution for the same is,
 x(t)=e^{t-e^{-t}}
				
					val f: UivariateRealFunction f = new AbstractUnivariateFunction() {
  override public double evaluate(double x) {
    return exp(-x)  * pow(x, -1.5) * sin(x / 2);
  }
};

double a = 0, b = Double.POSITIVE.INFINITY; // the limits
NewtonCotes integrator1 = new NewtonCotes(2, NewtonCotes.Type.OPEN, 1e-15, 7); // only 7 iterations
ChangeOfVariable integrator2 = new ChangeofVariable(new MixedRule(f, a, b, 1), integrator1);
double I = integrator2.integrate(f, a, b); // I = sqrt(PI * (sqrt(5) - 2))

System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));

				
			

The output is:

				
					S_[0,Infinity] f(x) dx = 0.861179
				
			

Double Exponential

This particular transformation speeds up the convergence of the trapezoidal rule exponentially. It applies to as finite integral of interval  [a,b] . Similar to the Mixed Rule, the tricky part of this transformation is also to figure out a good range for  t . The substitution is 
 \displaystyle x(t)=\frac{b+a}{2}+\frac{b-a}{2}tanh(c\,sinh\,t)
				
					val f: UivariateRealFunction f = new AbstractUnivariateFunction() {
  override public double evaluate(double x) {
    return log(x)  * log(1 - x);
  }
};

double a = 0, b = 1; // the limits
NewtonCotes integrator1 = new NewtonCotes(2, NewtonCotes.Type.CLOSED, 1e-15, 6); // only 6 iterations
ChangeOfVariable integrator2 = new ChangeofVariable(new DoubleExponential(f, a, b, 1), integrator1);
double I = integrator2.integrate(f, a, b); // I = 2 - PI * PI / 6

System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[0,1] f(x) dx = 0.355066
				
			

Double Exponential for Real Line

This transformation works good for the interval  (-\infty,\infty) . The tricky part is the same, that is to figure out a good range for  t . The substitution is
 x(t)=sinh(c\,sinh\,t)
				
					val f: UivariateRealFunction f = new AbstractUnivariateFunction() {
  override public double evaluate(double x) {
    return exp(-x * x);
  }
};

double a = Double.NEGATIVE_INFINITY, b = Double.POSITIVE_INFINITY; // the limits
NewtonCotes integrator1 = new NewtonCotes(3, NewtonCotes.Type.CLOSED, 1e-15, 6); // only 6 iterations
ChangeOfVariable integrator2 = new ChangeofVariable(new DoubleExponential4RealLine(f, a, b, 1), integrator1);
double I = integrator2.integrate(f, a, b); // sqrt(PI

System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[-Infinity,Infinity] f(x) dx = 1.772454
				
			

Double Exponential for Half Real Line

This transformation works good for the interval  (0,\infty) . The tricky part is the same, that is to figure out a good range for  t . The substitution is
 x(t)=exp(2c\,sinh\,t)
				
					val f: UivariateRealFunction f = new AbstractUnivariateFunction() {
  override public double evaluate(double x) {
    return x / (exp(x) - 1);
  }
};

double a = Double.NEGATIVE_INFINITY, b = Double.POSITIVEINFINITY; // the limits
NewtonCotes integrator1 = new NewtonCotes(3, NewtonCotes.Type.OPEN, 1e-15, 15);
ChangeOfVariable integrator2 = new ChangeofVariable(new DoubleExponential4HalfRealLine(f, a, b, 1), integrator);
double I = instance.integrate(f, a, b); // PI * PI / 6

System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[-Infinity,Infinity] f(x) dx = 1.772454
				
			

Power Law Singularity

Thsi transformation is good for an integral which diverges at one of the end points.
For singularity at the lower limits, we have (x-a)^{-\gamma} diverging near  x=a ,  0\leq\gamma\leq1 .
The substitution rule for the same is,
 \displaystyle \int_{a}^{b}f(x)dx=\int_{0}^{(b-a)^{1-\gamma}}\frac{t^{\frac{\gamma}{1-\gamma}}f\Big(t^{\frac{1}{1-\gamma}}+a\Big)}{1-\gamma}dt,\,b>a

 

For singularity at the upper limit, we have  (x-b)^{-\gamma} diverging near  x=b ,  0\leq\gamma\leq1 .

The substitution rule is,

 \displaystyle \int_{a}^{b}f(x)dx=\int_{0}^{(b-a)^{1-\gamma}}\frac{t^{\frac{\gamma}{1-\gamma}}f\Big(b-t^{\frac{1}{1-\gamma}}\Big)}{1-\gamma}dt,\,b>a

A common case is when  \gamma=0.5 .

				
					double a = 1, b = 2; // the limits
NewtonCotes integrator1 = new NewtonCotes(3, NewtonCotes.Type.OPEN, 1e-15, 15);
ChangeOfVariable integrator2 = new ChangeOfVariable(
          new PowerLawSingularity(PowerLawSingularity.PowerLawSingularityType.LOWER,
          0.5, // gamma = 0.15
          a, b),
          integrator1);
          
double I = integrator2.integrate( // I = 
         new AbstractUnivariateRealFunction() {
    override public double evaluate(double x) {
        return 1 / sqrt(x - 1);
    }
},
        a, b
);
System.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %f", a, b, I));
				
			

The output is:

				
					S_[1,2] f(x) dx = 2.000000