Integrasi: gauss quadrature, simpson, dan jumlah atas dan bawah

numeric_integartion.c

/*
	coder		: duken marga
	blog		: http://duken.info
	description	: integrasi numerik menggunakan motoda trapezoidal, 
			  gauss quadrature, simpson, dan jumlah atas dan bawah
*/
 
#include <stdio.h>
#include <math.h>
 
float jumlah_atas(float, float, float);
float jumlah_bawah(float, float, float);
float trapezoidal(float, float, float);
float simpson(float, float, float);
float gauss_quadrature(float, float);
 
float f(float);
 
int main(){
	float batas_atas;
	float batas_bawah;
	int jumlah_segmen = 20;
 
	printf("Masukkan batas atas : ");
	scanf("%f", &batas_atas);
	printf("Masukkan batas bawah : ");
	scanf("%f", &batas_bawah);
	printf("Masukkan jumlah segmen : ");
	scanf("%d", &jumlah_segmen);
 
	printf("Jumlah atas  : %20.18f \n",
		jumlah_atas(batas_bawah, batas_atas, jumlah_segmen));
	printf("Jumlah bawah : %20.18f \n",
		jumlah_bawah(batas_bawah, batas_atas, jumlah_segmen));
	printf("Rata-rata    : %20.18f \n\n",
		( jumlah_atas(batas_bawah, batas_atas, jumlah_segmen)
		+ jumlah_bawah(batas_bawah, batas_atas, jumlah_segmen) 
		) /2
		);
	printf("Jumlah trapezoidal  : %20.18f \n",
		trapezoidal(batas_bawah, batas_atas, jumlah_segmen));
	printf("Jumlah simpson  : %20.18f \n",
		simpson(batas_bawah, batas_atas, jumlah_segmen));
	printf("Jumlah gauss_quadrature  : %20.18f \n",
		gauss_quadrature(batas_bawah, batas_atas));
}
 
/*definisi fungsi yang akan dicari nilainya*/
float f(float x){
	return exp(x)*sin(x);
}
 
/* mencari nilai luas dengan metode jumlah atas */
float jumlah_atas(float bwh, float ats, float jlh){
	float i, luas = 0;
	float node1, node2;
	float y_node1, y_node2, y;
	float lebar_segmen = (ats - bwh) / jlh;
	for(i = bwh; i < ats; i += lebar_segmen){
		node1 = i;
		node2 = i + lebar_segmen;
		y_node1 = f(node1);
		y_node2 = f(node2);
		if(y_node1 < y_node2)
			y = y_node1;
		else
			y = y_node2;
 
		luas += lebar_segmen * y;
	}
	return luas;
}
 
/* mencari nilai luas dengan metode jumlah bawah */
float jumlah_bawah(float bwh, float ats, float jlh){
	float i, luas = 0;
	float node1, node2;
	float y_node1, y_node2, y;
	float lebar_segmen = (ats - bwh) / jlh;
	for(i = bwh; i < ats; i += lebar_segmen){
		node1 = i;
		node2 = i + lebar_segmen;
		y_node1 = f(node1);
		y_node2 = f(node2);
		if(y_node1 > y_node2)
			y = y_node1;
		else
			y = y_node2;
 
		luas += lebar_segmen * y;
	}
	return luas;
}
 
/* mencari nilai luas dengan metode trapezoidal */
float trapezoidal(float bwh, float ats, float jlh){
	int i;
	float luas = 0;
	float lebar_segmen = (ats - bwh) / jlh;
	float node, y;
 
	luas = 0.5 * lebar_segmen * (f(ats) + f(bwh));
	for(i = 1; i < jlh; i++){
		node = bwh + i * lebar_segmen;
		y = f(node);
		luas += lebar_segmen * y;
	}
	return luas;
}
 
/* mencari nilai luas dengan metode simpson 1/3 */
float simpson(float bwh, float ats, float jlh){
	int i;
	float luas = 0;
	float lebar_segmen = (ats - bwh) / jlh;
	float node, y;
 
	luas = lebar_segmen * (f(ats) + f(bwh)) / 3;
	for(i = 1; i < jlh; i++){
		node = bwh + i * lebar_segmen;
		y = f(node);
		if(i%2 == 1)
			y *= 4;
		else
			y *= 2;
		luas += lebar_segmen * y / 3;
	}
	return luas;
}
 
/* mencari nilai luas dengan metode gauss quadrature */
/* menggunakan orde 3, bisa diganti dengan orde lain */
/* yang perlu diganti : orde, weight[], x[] */
float gauss_quadrature(float bwh, float ats){
	int i;
	int orde = 3;
	float luas = 0;
	float t;
	float weight[] = { (float) 8/9, (float)5/9, (float)5/9 };
	float x[] = { 0, sqrt((float) 3/5), -sqrt((float) 3/5) };
 
	for(i = 0; i < orde; i++){
		t = ((ats - bwh) * x[i] + (ats + bwh)) / 2 ;
		luas += weight[i] * f(t);
	}
	luas *= (ats - bwh) / 2;
	return luas;
}