Approximate an Integral Using Simpson's Rule in C

Introduction

In this lab, you will learn how to approximate an integral using Simpson's Rule in C. The lab covers the step-by-step process of defining the function to integrate, the integration interval, and then implementing the Simpson's Rule algorithm to compute the approximate integral. By the end of this lab, you will have a working C program that can numerically integrate a given function over a specified interval using the Simpson's Rule method.

The lab consists of two main steps: defining the function and interval, and then using Simpson's Rule with even subintervals to compute the approximate integral. This hands-on exercise will help you gain practical experience in applying numerical integration techniques in C programming.

Skills Graph

Define f(x) and Interval [a,b]

In this step, we will set up the foundational components for implementing Simpson's Rule for numerical integration in C. We'll define the function to integrate and specify the integration interval.

First, let's create a new C file in the ~/project directory:

cd ~/project
nano simpson_integration.c

Now, let's write the initial code to define our function and interval:

#include <stdio.h>
#include <math.h>

// Function to integrate: f(x) = x^2
double f(double x) {
    return x * x;
}

int main() {
    // Define integration interval [a, b]
    double a = 0.0;  // Lower bound
    double b = 1.0;  // Upper bound

    printf("Integrating f(x) = x^2 from %f to %f\n", a, b);

    return 0;
}

Let's break down the code:

1. We include necessary headers: stdio.h for input/output and math.h for mathematical functions.
2. f(x) is defined as a simple quadratic function x^2, which we'll integrate.
3. In main(), we set the integration interval from 0 to 1.
4. We print out the interval for clarity.

Compile and run the code to verify:

gcc simpson_integration.c -o simpson_integration -lm
./simpson_integration

Example output:

Integrating f(x) = x^2 from 0.000000 to 1.000000

Use Simpson's Rule with Even Subintervals

In this step, we will implement the Simpson's Rule algorithm to approximate the definite integral of our function. We'll add the integration function to our existing C program.

Open the previous file:

cd ~/project
nano simpson_integration.c

Update the code to include the Simpson's Rule implementation:

#include <stdio.h>
#include <math.h>

// Function to integrate: f(x) = x^2
double f(double x) {
    return x * x;
}

// Simpson's Rule implementation
double simpsons_rule(double a, double b, int n) {
    double h = (b - a) / n;  // Width of each subinterval
    double sum = f(a) + f(b);  // First and last points

    // Compute sum of even and odd points
    for (int i = 1; i < n; i++) {
        double x = a + i * h;
        sum += (i % 2 == 0 ? 2 : 4) * f(x);
    }

    return (h / 3) * sum;
}

int main() {
    // Define integration interval [a, b]
    double a = 0.0;  // Lower bound
    double b = 1.0;  // Upper bound
    int n = 100;     // Number of subintervals (must be even)

    double integral = simpsons_rule(a, b, n);

    printf("Integral of f(x) = x^2 from %f to %f\n", a, b);
    printf("Approximation using Simpson's Rule with %d subintervals: %f\n", n, integral);

    return 0;
}

Compile and run the code:

gcc simpson_integration.c -o simpson_integration -lm
./simpson_integration

Example output:

Integral of f(x) = x^2 from 0.000000 to 1.000000
Approximation using Simpson's Rule with 100 subintervals: 0.333333

Let's break down the Simpson's Rule implementation:

1. simpsons_rule() takes the interval bounds and number of subintervals.
2. h calculates the width of each subinterval.
3. We start with the first and last points of the interval.
4. The loop adds weighted contributions from intermediate points.
5. Even points are multiplied by 2, odd points by 4.
6. The final result is scaled by h/3.

Print the Approximate Integral

In this final step, we'll enhance our Simpson's Rule implementation by adding more detailed output and comparing the numerical approximation with the exact integral value.

Open the previous file:

cd ~/project
nano simpson_integration.c

Update the code to include more comprehensive output:

#include <stdio.h>
#include <math.h>

// Function to integrate: f(x) = x^2
double f(double x) {
    return x * x;
}

// Exact integral calculation for f(x) = x^2 from a to b
double exact_integral(double a, double b) {
    return (pow(b, 3) - pow(a, 3)) / 3.0;
}

// Simpson's Rule implementation
double simpsons_rule(double a, double b, int n) {
    double h = (b - a) / n;  // Width of each subinterval
    double sum = f(a) + f(b);  // First and last points

    // Compute sum of even and odd points
    for (int i = 1; i < n; i++) {
        double x = a + i * h;
        sum += (i % 2 == 0 ? 2 : 4) * f(x);
    }

    return (h / 3) * sum;
}

int main() {
    // Define integration interval [a, b]
    double a = 0.0;  // Lower bound
    double b = 1.0;  // Upper bound
    int n = 100;     // Number of subintervals (must be even)

    // Compute approximation and exact integral
    double approx_integral = simpsons_rule(a, b, n);
    double exact_value = exact_integral(a, b);
    double error = fabs(exact_value - approx_integral);

    // Print detailed results
    printf("Integral Approximation Results:\n");
    printf("------------------------------\n");
    printf("Function: f(x) = x^2\n");
    printf("Interval: [%.2f, %.2f]\n", a, b);
    printf("Number of Subintervals: %d\n\n", n);

    printf("Approximation (Simpson's Rule): %.6f\n", approx_integral);
    printf("Exact Value: %.6f\n", exact_value);
    printf("Absolute Error: %.6f\n", error);
    printf("Relative Error: %.4f%%\n", (error / exact_value) * 100);

    return 0;
}

Compile and run the code:

gcc simpson_integration.c -o simpson_integration -lm
./simpson_integration

Example output:

Integral Approximation Results:
------------------------------
Function: f(x) = x^2
Interval: [0.00, 1.00]
Number of Subintervals: 100

Approximation (Simpson's Rule): 0.333333
Exact Value: 0.333333
Absolute Error: 0.000000
Relative Error: 0.0000%

Key additions in this step:

1. Added exact_integral() function to calculate the true integral value
2. Computed absolute and relative error
3. Created a formatted output showing detailed integration results

Summary

In this lab, we first defined the function to integrate, f(x) = x^2, and the integration interval [a, b] = [0, 1]. We then implemented the Simpson's Rule algorithm to approximate the definite integral of the function. The Simpson's Rule involves dividing the interval into an even number of subintervals, and computing the weighted sum of the function values at the endpoints and midpoints of these subintervals. Finally, we will print the approximate integral value obtained using the Simpson's Rule.