A javascript code for Romberg integration

As I am writing backend JavaScript for simulations in teaching Numerical Methods, I have also started developing some functions for some numerical techniques. Here is a function for Romberg integration.

function auto_integrator_trap_romb_hnm(func,a,b,nmax,tol_ae,tol_rae)
// INPUTS
// func=integrand
// a= lower limit of integration
// b= upper limit of integration
// nmax = number of partitions, n=2^nmax
// tol_ae= maximum absolute approximate error acceptable (should be >=0)
// tol_rae=maximum absolute relative approximate error acceptable (should be >=0)
// OUTPUTS
// integ_value= estimated value of integral

{
//Checking for input errors
	if (typeof a !== 'number') 
		{
		  throw new TypeError('<a> must be a number');
		}
    if (typeof b !== 'number') 
		{
		  throw new TypeError('<b> must be a number');
		}
    if ((!Number.isInteger(nmax)) || (nmax<1))
		{
		  throw new TypeError('<nmax> must be an integer greater than or equal to one.');
		}
	if ((typeof tol_ae !== 'number') || (tol_ae<0)) 
		{
		  throw new TypeError('<tole_ae> must be a number greater than or equal to zero');
		}
	if ((typeof tol_rae !== 'number') || (tol_rae<=0)) 
		{
		  throw new TypeError('<tole_ae> must be a number greater than or equal to zero');
		}
    
	var h=b-a
	// initialize matrix where the values of integral are stored
	
	var Romb = []; // rows
	for (var i = 0; i < nmax+1; i++) 
	{
		Romb.push([]);
		for (var j = 0; j < nmax+1; j++) 
		{
			Romb[i].push(math.bignumber(0)); 
		}
	}
	
	//calculating the value with 1-segment trapezoidal rule
	Romb[0][0]=0.5*h*(func(a)+func(b))
	var integ_val=Romb[0][0]
	
	for (var i=1; i<=nmax; i++)
	// updating the value with double the number of segments
	// by only using the values where they need to be calculated
	// See https://autarkaw.org/2009/02/28/an-efficient-formula-for-an-automatic-integrator-based-on-trapezoidal-rule/
	{
		h=0.5*h
		var integ=0
		for (var j=1; j<=2**i-1; j+=2)
		{
			var integ=integ+func(a+j*h)
		}
	
		Romb[i][0]=0.5*Romb[i-1][0]+integ*h
		// Using Romberg method to calculate next extrapolatable value
		// See https://young.physics.ucsc.edu/115/romberg.pdf
		for (k=1; k<=i; k++)
		{   
			var addterm=Romb[i][k-1]-Romb[i-1][k-1]
			addterm=addterm/(4**k-1.0)
			Romb[i][k]=Romb[i][k-1]+addterm

			//Calculating absolute approximate error
			var Ea=math.abs(Romb[i][k]-Romb[i][k-1])
			
			//Calculating absolute relative approximate error
			var epsa=math.abs(Ea/Romb[i][k])*100.0
			
			//Assigning most recent value to the return variable
			integ_val=Romb[i][k]
			
			// returning the value if either tolerance is met
			if ((epsa<tol_rae) || (Ea<tol_ae))
			{
				return(integ_val)
			}
		}
	}
	// returning the last calculated value of integral whether tolerance is met or not
	return(integ_val)
}

Here we are testing it for a typical integrand of f(x)=1/x. Take it for a spin and see how well it works. Make it even better.

<!DOCTYPE html>
<meta content="text/html;charset=utf-8" http-equiv="Content-Type">
<meta content="utf-8" http-equiv="encoding">
<html>
<head>
	<title>A test for the automatic integrator based on Romberg integration and trapezoidal rule</title>
	https://cdnjs.cloudflare.com/ajax/libs/mathjs/5.1.2/math.min.js
	http://trap_romberg_2021.js
</head>
<body>
<script>
// This program is written to test the romberg integration scheme that is used
// as an automatic integrator
// INPUTS
// a= lower limit of integration
// b= upper limit of integraton
// nmax= number of partitions, segment is then 2^nmax 
// tol_ae= tolerance on absolute approximate error 
// tol_rae=tolerance on percentage absolute relative approximate error 
var a=0.001
var b=10
var nmax=20
var tol_ea=0.0
var tol_rae=0.0000000005

var abc=auto_integrator_trap_romb_hnm(func,a,b,nmax,tol_ea,tol_rae)
console.log("romberg "+abc)

var exact=math.log(b)-math.log(a)
console.log("exact "+exact)
function func(x)
{
   //val=math.exp(-x)
   //var pi=4*math.atan(1.0)
   //var val=2/math.sqrt(pi)*math.exp(-x*x)
   var val=1/x
   return(val)
}
</script>
</body>
</html>

____________________________

This post is brought to you by

MATLAB code for the efficient automatic integrator

In the previous post, we discussed why doubling the number of segments in the automatic integrator based on multiple-segment trapezoidal rule is more efficient than increasing the number of segments one at a time. But this advantage involves having to store the individual function values from previous calculations and then having to retrieve them properly. This drawback was circumvented very efficiently by using the formula derived in another previous post where there is no need to store individual function values.

The matlab file for finding a definite integral by directly using the multiple segment trapezoidal rule from this post is given here (matlab file, html file), while the matlab file that uses the more efficient formula from this post is given here (matlab file, html file).  Here are the inputs to the programs.

% a = Lower limit of integration
% b = Upper limit of integration
%  nmax = Maximum number of segments
% tolerance = pre-specified tolerance in percentage
% f = inline function as integrand

a=5.3;
b=10.7;
nmax=200000;
tolerance=0.000005;
f=inline(‘exp(x)*sin(2*x)’)

We ran both the program on a PC and found that the more efficient algorithm (51 seconds) ran in half the time as the other one (82 seconds).  This is expected, as only n function evaluations are made for 2n-segments rule with the efficient formula, while 2n+1 functions evaluations are made for the original formula.

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://www.youtube.com/numericalmethodsguy.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

An efficient formula for an automatic integrator based on trapezoidal rule

In the previous post, we discussed why doubling the number of segments in the automatic integrator based on multiple-segment trapezoidal rule is more efficient than increasing the number of segments one at a time. But this advantage involves having to store the individual function values from previous calculations and then having to retrieve them properly. This drawback can be circumvented very efficiently as explained below. What you will see is that there is no need to store individual function values.

Automatic Integrator Formula
Automatic Integrator Formula

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://www.youtube.com/numericalmethodsguy.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Why keep doubling the segments for an automatic integrator based on Trapezoidal rule?

Automatic Integrator
Automatic Integrator

Automatic Integrator
Automatic Integrator

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://www.youtube.com/numericalmethodsguy.  

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

An automatic integrator using Trapezoidal rule

How would you know how many segments to use in a Trapezoidal rule of integration to get an accurate value of the integral?  This can be done by applying the Trapezoidal rule for 1 segment rule, then 2 segment rule, followed by 4 segment rule and so on.  As soon as the absolute relative approximate error (page 5-6) between the consecutive answers becomes less than the pre-specified tolerance chosen by the user, you would have your integral within the accuracy you desired.

Here is a MATLAB program that does that for you.  The MATLAB program that can be downloaded at http://numericalmethods.eng.usf.edu/blog/trapezoidal_rule_automatic.m (better to download it as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is here: http://numericalmethods.eng.usf.edu/blog/html/trapezoidal_rule_automatic.html

% Simulation : Using Trapezoidal rule as an automatic integrator

% Language : Matlab 2007a

% Authors : Autar Kaw, http://numericalmethods.eng.usf.edu

% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/trapezoidal_rule_automatic.m

% Last Revised : October 12, 2008

% Abstract: This program uses multiple-segment Trapezoidal
% rule to integrate f(x) from x=a to x=b within a pre-specified tolerance

clc
clear all

disp(‘This program uses multiple-segment Trapezoidal rule as an automatic integrator’)
disp(‘to integrate f(x) from x=a to x=b within a pre-specified tolerance’)
disp(‘ ‘)
disp(‘Author: Autar K Kaw.’)
disp(‘https://autarkaw.wordpress.com&#8217;)
disp(‘http://numericalmethods.eng.usf.edu&#8217;)
disp(‘ ‘)

%INPUTS.  If you want to experiment, these are the only variables
% you should and can change.
% a = Lower limit of integration
% b = Upper limit of integration
% nmax = Maximum number of segments
% tolerance = pre-specified tolerance in percentage
% f = inline function as integrand
a=5.3;
b=10.7;
nmax=20000;
tolerance=0.005;
f=inline(‘exp(x)*sin(2*x)’);

% SIMULATION
disp(‘INPUTS’)
func=[‘     The integrand is =’ char(f)];
disp(func)
fprintf(‘     Lower limit of integration, a= %g’,a)
fprintf(‘\n     Upper limit of integration, b= %g’,b)
fprintf(‘\n     Maximum number of segments, nmax = %g’,nmax)
fprintf(‘\n     Pre-specified percentage tolerance, eps = %g’,tolerance)
disp(‘  ‘)
disp(‘  ‘)

% Doing the automatic integration
% Calculating the integral using 1-segment rule
previous_integral=(b-a)/2*(f(a)+f(b));
% Initializing ea as greater than pre-specified tolerance for loop to work
ea=2*tolerance;
% Starting with 2-segments inside the while loop
n=2;
while (ea>tolerance) & (n<=nmax)
h=(b-a)/n;
% Keeping track of used number of segments
nused=n;
current_integral=0;
for i=1:1:n-1
current_integral=current_integral+f(a+i*h);
end
current_integral=2*current_integral+f(a)+f(b);
current_integral=(b-a)/(2*n)*current_integral;
% Calculating the absolute relative approximate error
ea = abs((current_integral-previous_integral)/current_integral)*100;
previous_integral=current_integral;
% Doubling the number of segments for next estimate of the integral
n=n*2;
end

disp(‘OUTPUTS’)
fprintf(‘      Number of segments used  =%g’, nused)
fprintf(‘\n      Approximate value of integral is =%g’,current_integral)
fprintf(‘\n      Absolute percentage relative approximate error =%g’, ea)
if (ea>tolerance)
disp(‘  ‘)
disp(‘  ‘)
disp(‘     NOTE: The value of integral is not within the pre-specified tolerance’)
end
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

%d bloggers like this: