Metoda prostokątów - całkowanie numeryczne.

Witam. Mialem do napisania program calkujacy. Zalozenia programu:

#include <stdio.h>
#include <stdlib.h>

float funkcja(float a, float x);
float pros(float xp, float xk, int n);
float calkowanie_pros(int n);


float a,p,xp,xk,dx,wartosc_calki;
int i,n;
int main()
{
 
  printf("Podaj poczatek przedzialu\n");
  scanf("%f", &xp);
  printf("Podaj koniec przedzialu\n");
  scanf("%f", &xk);

  printf("Podaj dokladnosc (ilosc prostokatow)\n");
  scanf("%d",&n);
  printf("\n\n");
  calkowanie_pros(n);
  system("PAUSE");   
  return 0;
}

float funkcja(float a, float x) //tutaj jest zawarta funkcja, ktora bedzie calkowana
{
  return x*5;
}
float pros(float xp, float xk, int n) //algorytm calkujacy moetoda prostokatow
{

  for(wartosc_calki=0,i=1;i<=n;i++)
  {
    dx=(xk-xp)/(float)n;                               
    wartosc_calki += funkcja(a,xp + i*dx);
  }
  wartosc_calki *= dx;
  return wartosc_calki;
 
}

   
float calkowanie_pros(int n)
// funkcja obliczajaca na ile prostokatow musi byc podzielona funkcja, aby uzyskac prog dokladnosci=0.1%
{
  do
  {
    p=(1-(pros(xp, xk, n) / pros(xp, xk, n-1)))*100;
    if (p >= 0.1) { ++n; continue;}
  } while (p>=0.1);
  printf("dla dokladnosci n=%d, wartosc calki(m. pros) jest rowna= %f \n",n,wartosc_calki);
}

Generalnie mam problem z progiem dokladnosci, nie jest wyliczany chyba zbyt dokladnie. Moglby ktos sprawdzic i wyrazic swoje opinie? Pozdrawiam

Błąd p który liczysz to tylko sprawdzenie o ile zmienia się przybliżana wartość całki przy zmianie ilości prostokątów o 1. Może się okazać, że dodając kolejny prostokąt szacowana wartość całki zmienia się już niewiele (np. poniżej 0.1%), ale nadal jesteś daleeeko od faktycznej (analitycznej) wartości calki.

Taki mały przykład do zabawy (obliczanie całki metodą prostokątów na pierwszym boku, na drugim i w środku, oraz porównanie z analityczną wartością).

 
#include <stdio.h>
#include <math.h>

float fun(float);

float calk_prost1(int);
float calk_prost2(int);
float calk_prost3(int);

float calk_analit();

float xp, xk;

int main() {
	int i;
	int n;

	// tymczasowe wartosci calki
	float Ct1, Ct2;

	// wartosci calki dla kolejnych metod
	float C1, C2, C3;

	// liczba prostokatow dla kolejnych metod
	int n1, n2, n3;

	float p;

	printf("Podaj poczatek przedzialu: ");
	scanf("%f", &xp);
	printf("Podaj koniec przedzialu: ");
	scanf("%f", &xk);
	printf("Podaj ilosc prostokatow: ");
	scanf("%d", &n);

	n1 = n;
	do {
		Ct1 = calk_prost1(n1);
		Ct2 = calk_prost1(n1-1);
		//p = ( 1 - Ct1/Ct2 ) * 100;
		p = fabsf( Ct2 - Ct1 ) / Ct2 * 100;
		if ( p > 0.1 ) {
			n1++;
			continue;
		}
	} while ( p > 0.1 );
	C1 = Ct1;

	n2 = n;
	do {
		Ct1 = calk_prost2(n2);
		Ct2 = calk_prost2(n2-1);
		//p = ( 1 - Ct1/Ct2 ) * 100;
		p = fabsf( Ct2 - Ct1 ) / Ct2 * 100;
		if ( p > 0.1 ) {
			n2++;
			continue;
		}
	} while ( p > 0.1 );
	C2 = Ct1;

	n3 = n;
	do {
		Ct1 = calk_prost3(n3);
		Ct2 = calk_prost3(n3-1);
		//p = ( 1 - Ct1/Ct2 ) * 100;
		p = fabsf( Ct2 - Ct1 ) / Ct2 * 100;
		if ( p > 0.1 ) {
			n3++;
			continue;
		}
	} while ( p > 0.1 );
	C3 = Ct1;

	printf("Po calkowaniu na pierwszym boku: C = %f dla n=%d\n", C1, n1 );
	printf("Po calkowaniu na drugim boku: C = %f dla n=%d\n", C2, n2 );
	printf("Po calkowaniu na w srodku: C = %f dla n=%d\n", C3, n3);
	printf("A tym czasem calka analityczna: C = %f\n", calk_analit() );

	return 0;
}

float fun(float x) {
	return x*x;
}

float calk_analit() {
	// calka analityczna: nalezy wpisac analityczny wzor na calke
	// z funkcji fun(x) w granicach xp , xk ...
	return (xk*xk*xk - xp*xp*xp)/3;
}

float calk_prost1(int n) {
	// calkowanie metoda prostokatow
	// z przyjeciem wartosci funkcji na pierwszym boku prostokata
	float dx = (xk-xp)/n;
	float x0 = 0;
	float calka = 0;
	float pole = 0;
	int i;

	//printf("dx = %f\n", dx);

	for (i=0; i<n; i++) {
		x0 = xp + i*dx; // x0 na pierwszym boku prostokata
		pole = dx * fun(x0);
		calka += pole;
		//printf("Pole %d prostokata = %f\n", i+1, pole);
	}

	return calka;
}

float calk_prost2(int n) {
	// calkowanie metoda prostokatow
	// z przyjeciem wartosci funkcji na drugim boku prostokata
	float dx = (xk-xp)/n;
	float x0 = 0;
	float calka = 0;
	float pole = 0;
	int i;

	//printf("dx = %f\n", dx);

	for (i=0; i<n; i++) {
		x0 = xp + dx + i*dx; // x0 na drugim boku prostokata
		pole = dx * fun(x0);
		calka += pole;
		//printf("Pole %d prostokata = %f\n", i+1, pole);
	}

	return calka;
}

float calk_prost3(int n) {
	// calkowanie metoda prostokatow
	// z przyjeciem wartosci funkcji w srodku prostokata
	float dx = (xk-xp)/n;
	float x0 = 0;
	float calka = 0;
	float pole = 0;
	int i;

	//printf("dx = %f\n", dx);

	for (i=0; i<n; i++) {
		x0 = xp + dx/2 + i*dx; // x0 w srodku prostokata
		pole = dx * fun(x0);
		calka += pole;
		//printf("Pole %d prostokata = %f\n", i+1, pole);
	}

	return calka;
}