Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Dec 27, 2009 #1
    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
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted