%% Copyright (c) 2025 Duong Ngoc Son <son@duongnotes.com>

%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
%% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
%% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
%% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
%% HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
%% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
%% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%' @title Simpson 1/3 rule.
%'
%' @description
%' Find approximately the integral of a function by Simpson 1/3 rule

function integral = Simpson(f, a, b, m)
	% Tham số: f - hàm lấy tích phân; a, b - hai cận lấy tích phân;
	% n = 2 * m + 1 - số các khoảng chia.
	
	if m == 0
		disp('Số khoảng chia không phù hợp.');
		return;
	end
	
	n = 2 * m + 1;
	h = (b - a) / (2 * m);
	
	% Chia đoạn [a, b] thành n đoạn nhỏ bằng nhau với độ dài h
	
	x = linspace(a, b, n);

	% Tính giá trị của hàm số dưới dấu tích phân tại các điểm chia
	
	y = zeros(1, n);
	for i = 1:n
		y(i) = feval(f, x(i));
	end
    
	% Bắt đầu thực hiện tính tích phân theo công thức Simpson 1/3
	
	temp = y(1) + y(n) + 4 * y(n - 1);
	for i = 2:2:(n - 3)
		temp = temp + 4 * y(i) + 2 * y(i + 1);
	end
	integral = temp * (h / 3);
end

% Ví dụ áp dụng: Hàm dưới dấu tích phân, f(x) = 1 / (1 + x), giá trị 
% đúng của tích phân I = ln(2) = 0,6931471806 ... 

f = @(x) 1 / (1 + x);
w = Simpson(f, 0, 1, 10);
fprintf('%d\n', w)

