0

Currently, I'm writing a program to calculate the Integral of an equation by using the Trapezoidal rule and a combination of these to archive higher precision. Right now as you can see below, I have hard-coded the function. Is it possible to read in a mathematical equation and evaluate it? I know I could read in the input character list and then evaluate the function (like if char[i] == '+' do ....) but is there an easier way? Thank you in advance!

void Integral_TN (double* TN_ptr,double a,double b,int n,int my_rank,int p){
     int i;
     double a_loc;
     double b_loc;        
     double hN;

    *TN_ptr = 0;
     hN = (b-a)/n;

     a_loc = a + my_rank*n/p*hN;     
     b_loc = a + (my_rank+1)*n/p*hN;
    *TN_ptr += (function(a_loc)+function(b_loc))/2;     /*Evaluate f at the borders*/

     for(i = 1; i < n/p; i++){
          *TN_ptr += function(a_loc + i*hN);            /*Evaluate f at the inner nodes*/
     }
    *TN_ptr = *TN_ptr*hN;
}

double function(double x){
      double y;
      y = 1/(1+x*x);
      return y; 
}
Swr7der
  • 849
  • 9
  • 28
Noveel
  • 65
  • 5
  • 2
    Perhaps have another go at formatting the code and make it more readable. PS What does _Thank you in advantage_ mean? – Ed Heal Nov 19 '16 at 13:42
  • Done. The important part is the function "function". Right now the function is hard coded I would like to read it in as a parameter, e.g. from the command line. – Noveel Nov 19 '16 at 13:54
  • Sounds like you need a parser. – EOF Nov 19 '16 at 13:57
  • This is a lot of work to be able to enter in arbitrary maths functions. Parser, reverse polish etc.. – Ed Heal Nov 19 '16 at 13:59
  • 1
    There is no easy way to do this since it would involve parsing. The simplest way might be to write a function which evaluates expressions given in space-delimited Reverse Polish Notation since it raises no precedence problems and is parentheses-free. – John Coleman Nov 19 '16 at 14:02
  • [Evaluate Mathematical Function from String](http://stackoverflow.com/q/9977045/995714) – phuclv Nov 20 '16 at 16:15

2 Answers2

1

There is no easier way for what you wish to achieve. If you want to apply a specific formula to some values, you have to define a function for it and input the values.

If you wish to type in the whole expression (with values and operators) as input and get the desired result as output, you would have to go beyond basic programming. You would need to create a parser.

For instance, if you provide 3+2*4 as input, you might expect 11 as output without reading in separate values 3, 2 and 4. This can be made possible by implementing a custom parser in a parser-generator like YACC. Basically, you would create and define new rules on how the input should be interpreted.

skrtbhtngr
  • 2,223
  • 23
  • 29
1

skrtbhtngr already answered the stated question, but I'd like to address the underlying problem.

Apply the Unix philosophy. Use one tool to generate the data, and another to consume it -- here, to compute the integral using the trapezoid rule.

The easiest format you can use is the same supported by Gnuplot:

  • Empty lines are ignored
  • Lines beginning with # are ignored, and can be used for comments
  • Each line defines one sample

Essentially, you could describe a sine curve, very roughly, using

#x  sin(x)
0.000  0.000000000
0.100  0.099833417
0.200  0.198669331
0.300  0.295520207
0.400  0.389418342
0.500  0.479425539
0.600  0.564642473
0.700  0.644217687
0.800  0.717356091
0.900  0.783326910
1.000  0.841470985
1.100  0.891207360
1.200  0.932039086
1.300  0.963558185
1.400  0.985449730
1.500  0.997494987
1.600  0.999573603
1.700  0.991664810
1.800  0.973847631
1.900  0.946300088
2.000  0.909297427
2.100  0.863209367
2.200  0.808496404
2.300  0.745705212
2.400  0.675463181
2.500  0.598472144
2.600  0.515501372
2.700  0.427379880
2.800  0.334988150
2.900  0.239249329
3.000  0.141120008
3.100  0.041580662
3.200 -0.058374143
3.300 -0.157745694
3.400 -0.255541102
3.500 -0.350783228
3.600 -0.442520443
3.700 -0.529836141
3.800 -0.611857891
3.900 -0.687766159
4.000 -0.756802495
4.100 -0.818277111
4.200 -0.871575772
4.300 -0.916165937
4.400 -0.951602074
4.500 -0.977530118
4.600 -0.993691004
4.700 -0.999923258
4.800 -0.996164609
4.900 -0.982452613
5.000 -0.958924275
5.100 -0.925814682
5.200 -0.883454656
5.300 -0.832267442
5.400 -0.772764488
5.500 -0.705540326
5.600 -0.631266638
5.700 -0.550685543
5.800 -0.464602179
5.900 -0.373876665
6.000 -0.279415498
6.100 -0.182162504
6.200 -0.083089403

You could use eg. awk to generate that, like I did:

awk 'BEGIN { printf "#x sin(x)\n" ; for (x=0.0; x<6.3; x+=0.1) printf "%.3f %11.9f\n", x, sin(x) }'

If you save that to a file (appending > data.txt to the above command), you can plot it in Gnuplot using

plot "data.txt" using 1:2 notitle with lines

Such data is easy to read in a C program. Because I only use POSIX.1 systems (Linux, BSDs, macOS), and POSIX.1 provides the very useful getline() function -- it lets you read in lines of any length, dynamically allocating a large enough buffer --, this particular implementation also requires POSIX.1 support. In other words, it works basically everywhere except in Windows.

#define _POSIX_C_SOURCE 200809L
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>

/* Read x y(x) values, one pair per line, from a stream.
   The arrays are dynamically allocated, and pointers stored
   to *xptr and *yptr. The size of the arrays is stored at *nptr.
   They are initially cleared to NULL/zero.
   The function returns 0 if success, an errno error code otherwise:
      EINVAL:  Invalid function parameters
      EIO:     Read error
      ENOMEM:  Out of memory
      EBADMSG: Malformed line
*/
int read_xy(double **xptr, double **yptr, size_t *nptr, FILE *in)
{
    /* Line input buffer variables. */
    char   *line = NULL;
    size_t  size = 0;
    ssize_t len;

    /* Data array variables. */
    double *x = NULL;
    double *y = NULL;
    size_t  n = 0;    /* Entries in x[] and y[] */
    size_t  nmax = 0; /* Entries allocated */

    /* Temporary variables. */
    double  xval, yval, *newx, *newy;

    /* We clear the output parameters to NULL or zero,
       in case the caller is careless and does not check
       the return value. Clearing them ensures they do
       not contain garbage in such a case. */
    if (xptr)
        *xptr = NULL;
    if (yptr)
        *yptr = NULL;
    if (nptr)
        *nptr = 0;

    /* We need in and nptr, and at least one of xptr and yptr. */
    if (!in || !nptr || (!xptr && !yptr))
        return errno = EINVAL;

    /* If an error has already occurred in 'in',
       we do not even try to read from it. */
    if (ferror(in))
        return EIO;

    while (1) {

        /* Read next input line. */
        len = getline(&line, &size, in);

        /* End of input or error? */
        if (len < 0)
            break;

        /* Skip empty and comment lines. */
        if (len == 0 ||
            line[0] == '\n' || (line[0] == '\r' && line[1] == '\n') ||
            line[0] == '#')
            continue;

        /* Parse the line. */
        if (sscanf(line, " %lf %lf", &xval, &yval) != 2)
            break;

        /* Need to grow the dynamically allocated arrays? */
        if (n >= nmax) {

            /* Allocation policy.
               We allocate room for at least 16 doubles,
               then double the size up to 1048576 (=2^20),
               then adjust to the next full multiple of 1048576.
               This is not 'the best', but it is robust,
               and not too wasteful.
            */
            if (n < 16)
                nmax = 16;
            else
            if (n < 1048576)
                nmax = n * 2;
            else
                nmax = (n | 1048575) + 1048576;

            /* Note: realloc(NULL, size) is equivalent to malloc(size).
               If the realloc() call fails, it returns NULL,
               but the original array is still valid.
               Also note that free(NULL) is safe, and does nothing.
            */
            newx = realloc(x, nmax * sizeof x[0]);
            newy = realloc(y, nmax * sizeof y[0]);
            if (newx)
                x = newx;
            if (newy)
                y = newy;
            if (!newx || !newy) {
                /* One or both of the allocations failed. */
                free(line);
                free(x);
                free(y);
                return ENOMEM;
            }
        }

        /* Save the parsed values to the arrays. */
        x[n] = xval;
        y[n] = yval;
        n++;
    }

    /* We no longer need the line buffer. */
    free(line);

    /* Did a read error occur? */
    if (ferror(in)) {
        free(x);
        free(y);
        return EIO;
    }

    /* Was there no data to read? */
    if (n < 1) {
        free(x);
        free(y);
        return 0;
    }

    /* Reallocate the arrays to their exact sizes
       (actually, allow for one extra double at the end,
        because it is often useful to copy the initial
        ones there if the data is considered cyclic).
    */
    nmax = n + 1; /* One extra just because it is so often useful. */
    newx = realloc(x, nmax * sizeof x[0]);
    newy = realloc(y, nmax * sizeof y[0]);
    if (newx)
        x = newx;
    if (newy)
        y = newy;
    if (!newx || !newy) {
        free(x);
        free(y);
        return ENOMEM;
    }

    /* Save the array pointers. */
    if (xptr)
        *xptr = x;
    else
        free(x);

    if (yptr)
        *yptr = y;
    else
        free(y);

    /* Save the number of samples read. */
    *nptr = n;

    /* If feof(in) is true, then we read everything
       up to end of input. Otherwise, we stopped at
       a line we could not parse.
    */   
    if (!feof(in))
        return EBADMSG;

    return 0;
}

That function, or something like it, should be in the course materials for every numerical computation course. They are just so darn useful. This particular one has no inherent limits to the size of data it can read, other than possible memory allocation limits set by the system administrator for each process. I know for a fact that it happily reads billions of lines of data, successfully, if you just have enough RAM.

Using the function is very simple. Here is an example main() that just reads such data from standard input -- remember that you can make it read from a file by appending < file to the command when you run it --, and prints the data out.

int main(void)
{
    double *x, *y;
    size_t  i, n;
    int     result;

    result = read_xy(&x, &y, &n, stdin);

    switch (result) {
    case 0: /* No errors */
        break;

    case EBADMSG:
        if (n > 1)
            fprintf(stderr, "Invalid line after %zu data samples.\n", n);
        else
            fprintf(stderr, "Cannot parse first input line.\n");
        return EXIT_FAILURE;

    case ENOMEM:
        fprintf(stderr, "Out of memory.\n");
        return EXIT_FAILURE;

    case EIO:
        fprintf(stderr, "Read error.\n");
        return EXIT_FAILURE;

    case EINVAL:
        fprintf(stderr, "Invalid parameters to the read_xy() function!\n");
        return EXIT_FAILURE;

    default:
        fprintf(stderr, "%s.\n", strerror(result));
        return EXIT_FAILURE;
    }

    printf("Read %zu samples:\n", n);
    for (i = 0; i < n; i++)
        printf("%.9f %.9f\n", x[i], y[i]);

    return EXIT_SUCCESS;
}

Note that most of it is the switch (result) { .. } error-reporting code. This example is very careful to tell you, if an error occurs; the reason is, you, as the user, need to know when the program knows they might spew garbage, and will in real life prefer for the programs to abort than spew that garbage silently -- perhaps making you believe it is valid.

While it is possible to modify the above code to work even on Windows (you replace the getline() with fgets(), for example, and hope that the buffer size you use for it will suffice; and also may have to change some of the errno error codes). However, there is a reason there are no Windows machines in the Top 500 supercomputer list: POSIXy systems (Unix and Linux) are just better suited for scientific computing. So, if you do intend to work on some scientific computing, you might as well just set up a Linux or BSD virtual machine, and do your development there.

Community
  • 1
  • 1
Nominal Animal
  • 38,216
  • 5
  • 59
  • 86
  • Thank you this rely helps me as a starting point. Since I use Kubuntu anyway I don't have to change the code. – Noveel Nov 20 '16 at 14:34