1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Forward difference method for heat equation

  1. Apr 11, 2009 #1
    I don't be able to convert the following code(HEAT EQUATION BACKWARD-DIFFERENCE ALGORITHM in the Burden-faires numerical analysis book).I need heat EQUATION FORWARD-DIFFERENCE ALGORITHM C like following code.I dont be able to convert FORWARD-DIFFERENCE the following code .Please help me.

    if I write VV = - [ALPHA * ALPHA * K / ( H * H )]; instead of VV = ALPHA * ALPHA * K / ( H * H );in the heat EQUATION BACKWARD-DIFFERENCE ALGORITHM then have I obtained HEAT EQUATION FORWARD-DIFFERENCE ALGORITHM?

    Code (Text):

    /*
    *   HEAT EQUATION BACKWARD-DIFFERENCE ALGORITHM 12.2
    *
    *   To approximate the solution to the parabolic partial-differential
    *   equation subject to the boundary conditions
    *                  u(0,t) = u(l,t) = 0, 0 < t < T = max t,
    *   and the initial conditions
    *                  u(x,0) = F(x), 0 <= x <= l:
    *
    *   INPUT:   endpoint l; maximum time T; constant ALPHA; integers m, N.
    *
    *   OUTPUT:  approximations W(I,J) to u(x(I),t(J)) for each
    *            I = 1, ..., m-1 and J = 1, ..., N.
    */

    #include<stdio.h>
    #include<math.h>
    #define pi 4*atan(1)
    #define true 1
    #define false 0

    double F(double X);
    void INPUT(int *, double *, double *, double *, int *, int *);
    void OUTPUT(double, double, int, double *, double);

    main()
    {
       double W[25], L[25], U[25], Z[25];
       double FT,FX,ALPHA,H,K,VV,T,X;
       int N,M,M1,M2,N1,FLAG,I1,I,J,OK;

       INPUT(&OK, &FX, &FT, &ALPHA, &N, &M);
       if (OK) {
          M1 = M - 1;
          M2 = M - 2;
          N1 = N - 1;
          /* STEP 1 */
          H = FX / M;
          K = FT / N;
          VV = ALPHA * ALPHA * K / ( H * H );
          /* STEP 2 */
          for (I=1; I<=M1; I++) W[I-1] = F( I * H );
          /* STEP 3 */
          /* STEPS 3 through 11 solve a tridiagonal linear system
             using Algorithm 6.7 */
          L[0] = 1.0 + 2.0 * VV;
          U[0] = -VV / L[0];
          /* STEP 4 */
          for (I=2; I<=M2; I++) {
             L[I-1] = 1.0 + 2.0 * VV + VV * U[I-2];
             U[I-1] = -VV / L[I-1];
          }
          /* STEP 5 */
          L[M1-1] = 1.0 + 2.0 * VV + VV * U[M2-1];
          /* STEP 6 */
          for (J=1; J<=N; J++) {
             /* STEP 7 */
             /* current t(j) */
             T = J * K;
             Z[0] = W[0] / L[0];
             /* STEP 8 */
             for (I=2; I<=M1; I++)
                Z[I-1] = ( W[I-1] + VV * Z[I-2] ) / L[I-1];
             /* STEP 9 */
             W[M1-1] = Z[M1-1];
             /* STEP 10 */
             for (I1=1; I1<=M2; I1++) {
                I = M2 - I1 + 1;
                W[I-1] = Z[I-1] - U[I-1] * W[i];
             }
          }
          /* STEP 11 */
          OUTPUT(FT, X, M1, W, H);
       }
       /* STEP 12 */
       return 0;
    }

    /* Change F for a new problem */
    double F(double X)
    {
       double f;

       f =  sin(pi * X);
       return f;
    }

    void INPUT(int *OK, double *FX, double *FT, double *ALPHA, int *N, int *M)
    {
       int FLAG;
       char AA;

       printf("This is the Backward-Difference Method for Heat Equation.\n");
       printf("Has the function F been created immediately\n");
       printf("preceding the INPUT procedure? Answer Y or N.\n");
       scanf("\n%c", &AA);
       if ((AA == 'Y') || (AA == 'y')) {
          printf("The lefthand endpoint on the X-axis is 0.\n");
          *OK =false;
          while (!(*OK)) {
             printf("Input the righthand endpoint on the X-axis.\n");
             scanf("%lf", FX);
             if (*FX <= 0.0)
                printf("Must be positive number.\n");
             else *OK = true;
          }
          *OK = false;
          while (!(*OK)) {
             printf("Input the maximum value of the time variable T.\n");
             scanf("%lf", FT);
             if (*FT <= 0.0)
                printf("Must be positive number.\n");
             else *OK = true;
          }
          printf("Input the constant alpha.\n");
          scanf("%lf", ALPHA);
          *OK = false;
          while (!(*OK)) {
             printf("Input integer m = number of intervals on X-axis\n");
             printf("and N = number of time intervals - separated by a blank.\n");
             printf("Note that m must be 3 or larger.\n");
             scanf("%d %d", M, N);
             if ((*M <= 2) || (*N <= 0))
                printf("Numbers are not within correct range.\n");
             else *OK = true;
          }
       }  
       else {
          printf("The program will end so that the function F can be created.\n");
          *OK = false;
       }  
    }

    void OUTPUT(double FT, double X, int M1, double *W, double H)
    {
       int I, J, FLAG;
       char NAME[30];
       FILE *OUP;

       printf("Choice of output method:\n");
       printf("1. Output to screen\n");
       printf("2. Output to text file\n");
       printf("Please enter 1 or 2.\n");
       scanf("%d", &FLAG);
       if (FLAG == 2) {
          printf("Input the file name in the form - drive:name.ext\n");
          printf("for example:   A:OUTPUT.DTA\n");
          scanf("%s", NAME);
          OUP = fopen(NAME, "w");
       }
       else OUP = stdout;
       fprintf(OUP, "THIS IS THE BACKWARD-DIFFERENCE METHOD\n\n");
       fprintf(OUP, "  I        X(I)    W(X(I),%12.6e)\n", FT);
       for (I=1; I<=M1; I++) {
          X = I * H;
          fprintf(OUP, "%3d %11.8f    %14.8f\n", I, X, W[I-1]);
       }
       fclose(OUP);
    }
     
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

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