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

Author: Autar Kaw

Autar Kaw (http://autarkaw.com) is a Professor of Mechanical Engineering at the University of South Florida. He has been at USF since 1987, the same year in which he received his Ph. D. in Engineering Mechanics from Clemson University. He is a recipient of the 2012 U.S. Professor of the Year Award. With major funding from NSF, he is the principal and managing contributor in developing the multiple award-winning online open courseware for an undergraduate course in Numerical Methods. The OpenCourseWare (nm.MathForCollege.com) annually receives 1,000,000+ page views, 1,000,000+ views of the YouTube audiovisual lectures, and 150,000+ page views at the NumericalMethodsGuy blog. His current research interests include engineering education research methods, adaptive learning, open courseware, massive open online courses, flipped classrooms, and learning strategies. He has written four textbooks and 80 refereed technical papers, and his opinion editorials have appeared in the St. Petersburg Times and Tampa Tribune.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s