博客

用Javascript实现了函数积分的算法。

从计算结果可以看出，算法精度受随机数发生器的影响较大。

//Monte carro求给定函数的定积分
//给定一个函数
function f1(x, xmin, xmax)
{
	if (x >= xmin && x <= xmax)
		return x * x * Math.sqrt(1 + x * x * x) - 1;	
}

function f2(x, xmin, xmax)
{
	if (x >= xmin && x <= xmax)
		return x * Math.log(1 + x);	
}

function f3(x, xmin, xmax)
{
	if (x >= xmin && x <= xmax)
		return Math.sqrt(1 - x * x);
}

function f4(x, xmin, xmax)
{
	if (x >= xmin && x <= xmax)
		return Math.cos(x) - 0.9;
}


//取得函数的上下极限——其精确度影响很大---由于用于积分，故若y>0，底限0；若y<0，则高限0；若y穿越x轴，则不修改默认极限
function getLimit(func, xmin, xmax)
{	
	//边界极限
	var t1 = func(xmin, xmin, xmax), t2 = func(xmax, xmin, xmax);
	var limit = {y0: Math.min(t1, t2), y1: Math.max(t1, t2)};			
	var delta = 1e-3, lastX, lastY, y;
	for(var x = xmin, i = 0; x <= xmax; x += delta, i++)
	{		
		y = func(x, xmin, xmax);			
		if (y < limit.y0)
			limit.y0 = y;		
		if (y > limit.y1)
			limit.y1 = y;
		if (i == 0)
		{
			lastX = x;
			lastY = y;
		}
		else
		{
			var tx = (lastX + x) / 2;
			var ty = func(tx, xmin, xmax);			
			if (Math.abs(lastY) > Math.abs(y))
			{
				if (Math.abs(lastY) < Math.abs(ty))
				{
					lastX = tx;
					lastY = ty;				
					if (ty < 0 && limit.y0 > ty)
						limit.y0 = ty;
					else if (ty >= 0 && limit.y1 <= ty)
						limit.y1 = ty;
				}				
			}			
			if (Math.abs(lastY) < Math.abs(y))
			{
				if (Math.abs(y) < Math.abs(ty))
				{
					lastX = tx;
					lastY = ty;				
					if (ty < 0 && limit.y0 > ty)
						limit.y0 = ty;
					else if (ty >= 0 && limit.y1 <= ty)
						limit.y1 = ty;
				}				
			}						
		}
	}
	if (limit.y0 > 0)
		limit.y0 = 0;
	if (limit.y1 < 0)
		limit.y1 = 0;
	return limit;
}

//积分
function cumulate(func, N, xmin, xmax)
{
	var limit = getLimit(func, xmin, xmax);
	//console.log('Limit=', limit.y0, limit.y1);
	var x, yr, y;
	var cntA = 0, cntB = 0;
	for(var i = 0; i < N; i++)
	{
		x = xmin + (xmax - xmin) * Math.random();
		yr = limit.y0 + (limit.y1 - limit.y0) * Math.random();
		y = func(x, xmin, xmax);
		if ( y >= 0 && yr >= 0 && yr <= y)
			cntA++;
		else
		if (y <= 0 && yr <= 0 && yr >= y)
			cntB++;
	}
	//console.log('counts=', cntA, cntB);	
	return (cntA - cntB) * (xmax - xmin) * (limit.y1 - limit.y0) / N;
}

//测试
var t1, t2;
var rst;
var N = 1e7;
var fs = [f1, f2, f3, f4];
var ans = [-0.593682861167513, 0.25, Math.PI/4, -0.058529015192103493347497678369701];
for(var i = 0; i < 4; i++)
{
	t1 = new Date();
	rst = cumulate(fs[i], N, 0, 1);
	t2 = new Date();
	console.log('[', i, ']=', ans[i].toFixed(15), rst.toFixed(15), t2 - t1, 'ms');
}