Simpson's rule in C++, FORTRAN, and Python [algorithm] by 바죠


https://en.wikipedia.org/wiki/Simpson%27s_rule
http://incredible.egloos.com/3021478

#include <iostream>
#include <cmath>
using namespace std;

double func(double x)
{
        return exp(x);
}

double simpson(int n, double a, double b, double (* f) (double ) )
{
        int m, i;
        double x, h;
        double sss;


        m=2*n ;
        h=(b-a)/(double ) m ;


        sss=f(a);
        for (i=1 ; i < m ; i=i+2)
        {
                x=a+h * (double ) i ;
                sss=sss+4.e0* f(x);
        }
        for (i=2 ; i < m-1 ; i=i+2)
        {
                x=a+h * (double ) i ;
                sss=sss+2.e0* f(x);
        }
        sss=sss+f(b);
        sss=h*sss/3.e0 ;


        return sss;
}

int main()
{
        int n ;
        double a,b;
        double sss;

        n=1000;
        a=0.e0;
        b=1.e0;
        sss=simpson(n,a,b, func);


        sss=func(1.)-1.e0;

        cout << sss << endl;

        return 0 ;
}





 g++ simpson.cpp
[ihlee@helix CPP_train]$ a.out
1.71828



 

---------------------------------



#!/usr/bin/python
from math import exp

def func(x):
  return exp(x)

def simpson(a,b,n,f):
  m=2*n
  h=(b-a)/m
  sss=f(a)

  for i in range(1,m-0,2):
     x=a+h* i
     sss=sss+ 4* f(x)
  for i in range(2,m-1,2):
     x=a+h* i
     sss=sss+ 2* f(x)
  sss=sss+f(b)
  sss=h*sss/3.
  return sss


a=0.
b=1.
n=1000
print simpson(a,b,n,func)
print func(1.)-1.


simpson.py
1.71828182846
1.71828182846

---------------------------------


!234567890
       implicit none
       real*8, external :: func
       real*8 rslt,aa,bb
       integer n

       n=1000
       aa=0.0d0
       bb=1.0d0
       call simpson(func,n,aa,bb,rslt)
       write(6,*) rslt,' rslt'

       stop
       end

       subroutine simpson(f,n,aa,bb,rslt)
       implicit none
       real*8 f
       integer n
       real*8 rslt,aa,bb
       real*8 h,xx
       integer j
       logical lodd
       integer m

       m=2*n
       h=(bb-aa)/float(m)
       rslt=(f(aa)+f(bb))
       lodd=.true.
       do j=1,m-1
        xx=aa+h*float(j)
       if(lodd)then
       rslt=rslt+4.0d0*f(xx)
       else
         rslt=rslt+2.0d0*f(xx)
       endif
       lodd=(.not. lodd)
       enddo

       rslt=rslt*h/3.0d0
       return
       end

       real*8 function func(x)
       implicit none
       real*8 x

       func=exp(x)
       return
       end




 a.out
    1.718281828459045       rslt
FORTRAN STOP



# C++에서 C언어의 표준함수 호출하기

- c를 더하고 .h를 빼라.
#include <stdio.h>  →  #include <cstdio>
#include <stdlib.h> →  #include <cstdlib>
#include <math.h>  →  #include <cmath>
#include <string.h> →  #include <cstring>


Computer Programming 2008:
Computer_programming_2008.pdf

핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax