I am trying to calculate a integral in C and i am trying to use Simpson's 3/8 Rule (Also known as Simpson's second Rule) to calculate it. I wrote a program but it gives the wrong result. When i used the same rule on paper it works. I don't know what's wrong.
simpsons.cpp :

#include <iostream>

double f(int x) {
  return 3*(x*x) 2*x 5;
}

float integralf(int lower_end, int high_end) {
  int a = lower_end;
  int b = high_end;
  double area = 0.0;
  area = ((b - a)/8)*(f(a)   3*f((2*a b)/3)   3*f((a 2*b)/3)   f(b));
  return area;
}

int main() {
  int lower = -3;
  int high = 5;
  std::cout << integralf(lower, high) << std::endl;

  return 0;
}

It gives 194 but answer is 208. I tried the same integral with Reimann sum technique and it gave 166.24.
reimann.cpp :

#include <iostream>

double f(int x) {
  return 3*(x*x) 2*x 5;
}

double integral(double low_end, double high_end, double steps) {
  double step = (high_end - low_end)/steps;
  double area = 0.0;
  double y = 0;
  double x = 0;

  for(int i=0;i<=steps;  i) {
    x = low_end   i * step;
    y = f(x);
    area  = y * step;
  }

  return area;
}

int main() {
  int low_end = -3;
  int high_end = 5;
  int steps = 100;
  std::cout << integral(low_end, high_end, steps) << std::endl;

  return 0;
}

I don't know what is wrong. Is there a better way to calculate integrals in cpp? Every advice is welcome.

SOLUTION

Edit : I did it. Final version of code (simpsons.cpp) :

#include <iostream>

double f(double x) {
  return 3*(x*x) 2*x 5;
}

float integralf(double lower_end, double high_end) {
  double a = lower_end;
  double b = high_end;
  double h = (b - a)/3;
  double area = 0.0;
  area = ((3*h)/8)*(f(a)   3*f((2*a b)/3)   3*f((a 2*b)/3)   f(b));
  return area;
}

int main() {
  int lower = -3;
  int high = 5;
  std::cout << integralf(lower, high) << std::endl;

  return 0;
}

It gave the answer 208.

CodePudding user response：

Final version of code (simpsons.cpp) :

#include <iostream>

double f(double x) {
  return 3*(x*x) 2*x 5;
}

float integralf(double lower_end, double high_end) {
  double a = lower_end;
  double b = high_end;
  double h = (b - a)/3;
  double area = 0.0;
  area = ((3*h)/8)*(f(a)   3*f((2*a b)/3)   3*f((a 2*b)/3)   f(b));
  return area;
}

int main() {
  int lower = -3;
  int high = 5;
  std::cout << integralf(lower, high) << std::endl;

  return 0;
}

It gave the answer 208.