qsimp.c
#include<stdio.h>
#include<math.h>
float func(float x)
{
return pow((1-x*x),3.0/2.0) * cos(x);
}
float qsimp(float (*f)(float),float a,float b) {
/* obliczenie calki funkcji f() (zadanej przez uzytkownika) w
granicach od a do b przy uzyciu wzoru Simpsona. wymaga
<math.h>. */
int n,k;
float x,h,t,i,oi,ot,s;
/* maxymalna dopuszczalna liczba przedzialow (musi byc potega 2) */
const int m=1024;
/* wstepne obliczenie pola trapezu opartego o granice calkowania */
n=1;
h=b-a;
t=h*(f(a)+f(b))/2;
i=t;
/* wzor Simpsona na kolejnych podprzedzialach */
do {
n *= 2;
ot=t;
oi=i;
h /= 2;
s = 0;
for(k=0; k < n/2; k++)
s += f(a+(2*k+1)*h);
t /= 2;
t += h*s;
i = (4*t-ot)/3;
} while(n <= m || fabs((i-oi)/i) > 1e-4);
return i;
}
int main( void )
{
printf( "wynik == %f\n", qsimp( func, -1, 0 ) );
printf( "wynik == %f\n", qsimp( func, 0, 1 ) );
return 0;
}
qtrap.c
#include<stdio.h>
#include<math.h>
float func(float x)
{
return pow((1-x*x),3.0/2.0) * cos(x);
}
float qtrap(float (*f)(float),float a,float b) {
/* obliczenie calki funkcji f() (zadanej przez uzytkownika) w
granicach od a do b przy uzyciu wzoru trapezow. wymaga
<math.h>. */
int n,k;
float x,h,i,oi,s;
/* maxymalna dopuszczalna liczba przedzialow (musi byc potega 2) */
const int m=128;
/* obliczenie pola trapezu opartego o granice calkowania */
n=1;
h=b-a;
i=h*(f(a)+f(b))/2;
/* wzor trapezow na kolejnych podprzedzialach */
do {
n *= 2;
oi=i;
h /= 2;
s = 0;
for(k=0; k < n/2; k++)
s += f(a+(2*k+1)*h);
i /= 2;
i += h*s;
} while(n <= m || fabs((i-oi)/i) > 1e-4);
return i;
}
int main( void )
{
printf( "wynik == %f\n", qtrap( func, -1, 1 ) );
return 0;
}
Add a comment: