# Help me with a program i wrote that finds double integrals over general regions

1. Dec 27, 2009

### okkvlt

edit: i found the problem. a stupid but easy 2 miss typo. heres the fixed code anyway, in case anybody wants 2 use it. note that the edge smoothing on the border of the upper function is not implemented yet. ull have to do that urself.

#include <stdio.h>
#include <math.h>
/*define f[x]*/
long double z(long double x, long double y)
{
long double z=x*y;
return z;
}
long double y1(long double x)
{
long double y1=2*x;
return y1;
}
long double y0(long double x)
{
long double y0=x*x;
return y0;
}

/*program main*/

int main()
{
double x1,x0,dy,dx, test1,test2,test3;
long double totalx=0, totaly=0,total1top=0,total2top=0,total1bot=0,total2bot=0,areatop,areabot, centerx1top=0,centery1top=0, centerx2top=0,centery2top=0, centerx1bot=0,centery1bot=0, centerx2bot=0,centery2bot=0,centerxtop,centerytop,centerxbot,centerybot,voledge=0;
long int nx,ny,i=0,j=1, k=0;
printf("x1 y1\nS S f[x,y]dydx\nx0 y0\n\n");
printf("x1=");
scanf("%lf",&x1);
printf("x0=");
scanf("%lf",&x0);
printf("number of x increments");
scanf("%d",&nx);
dx=(x1-x0)/nx;
printf("size of dx=%lf\n",dx);
printf("size of dy=");
scanf("%lf",&dy);
long double xval,yval, ytop[nx];
/*calculations*/

while(i<=nx-1)
{
xval=x0+i*dx;
ny=(y1(xval)-y0(xval))/dy;
while(j<=ny-2)
{

totaly=totaly+dy*((z(xval,y0(xval)+j*dy)+z(xval+dx,y0(xval)+j*dy)+z(xval,y0(xval)+j*dy+dy)+z(xval+dx,y0(xval)+j*dy+dy))/4+z(xval+.5*dx,y0(xval)+j*dy+.5*dy))/2;

j=j+1;
}
ytop=y0(x0+i*dx)+(ny-1)*dy;
totalx=totalx+dx*totaly;
totaly=0;
j=1;
printf("%d/%d\r",i,nx);
i=i+1;
}
double integral=totalx;
printf("preliminary integral=%.13lf\n smoothing edges\n", integral);
long double h=dx/1000;
long int n=1000;
while(k<=nx-1)
{
i=1;
xval=x0+k*dx;

while(i<=n-1)
{
total1bot=h*y0(xval+i*h)+total1bot;
centerx1bot=h*(xval+i*h)*y0(xval+i*h)+centerx1bot;
centery1bot=h*pow(y0(xval+i*h),2)+centery1bot;
i=i+1;
}
i=0;
while(i<=n-1)
{
total2bot=h*y0(xval+(i+.5)*h)+total2bot;
centerx2bot=h*(xval+(i+.5)*h)*y0(xval+(i+.5)*h)+centerx2bot;
centery2bot=h*pow(y0(xval+(i+.5)*h),2)+centery2bot;
i=i+1;
}

areabot=total1bot/3+2*total2bot/3+h*(y0(xval)+y0(xval+dx))/6;
areabot=(y0(xval)+dy)*(x0+dx)-(y0(xval)+dy)*(x0)-areabot;
centerxbot=centerx1bot/3+2*centerx2bot/3+h*(y0(xval)*xval+y0(xval+dx)*(xval+dx))/6;
centerxbot=.5*(y0(xval)+dy)*pow(xval+dx,2)-.5*(y0(xval)+dy)*pow(xval,2)-centerxbot;
centerxbot=centerxbot/areabot;
centerybot=centery1bot/3+2*centery2bot/3+h*(pow(y0(xval),2)+pow(y0(xval+dx),2))/6;
centerybot=pow((y0(xval)+dy),2)*(xval+dx)-pow((y0(xval)+dy),2)*(xval)-centerybot;
centerybot=centerybot/(areabot*2);
voledge=voledge+fabs(areabot)*z(centerxbot,centerybot);
k=k+1;
centery1bot=0;
centery2bot=0;
centerx2bot=0;
centerx1bot=0;
total1bot=0;
total2bot=0;

}
integral=integral+voledge;
printf("%.13lf\n", integral);

getchar();
getchar();
}

Last edited: Dec 27, 2009