!************************************************************** !* Program to demonstrate integration of a real function F(x) * !* by Simpson's method * !* ---------------------------------------------------------- * !* REFERENCE: "Mathematiques en Turbo-Pascal (Part 1) * !* By Marc Ducamp and Alain Reverchon, Eyrolles, * !* Paris, 1987" [BIBLI 03]. * !* ---------------------------------------------------------- * !* SAMPLE RUN: * !* (Integrate sin(x) from x=0 to x=1) * !* * !* Integral of a function F(X) by Simpson's method * !* * !* Input begin and end values for x variable: * !* * !* X0 = 0 * !* X1 = 1 * !* * !* Number of integration steps: 100 * !* * !* * !* Value of integral: 0.459697694133457 * !* * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !************************************************************** PROGRAM Test_Simpson real*8 x0 !begin x value real*8 x1 !end x value integer nstep !number of integration steps real*8 result !result of integral print *,' ' print *,' Integral of a function F(X) by Simpson''s method' print *,' ' print *,' Input begin and end values for x variable:' print *,' ' write(*,10,advance='no'); read *, x0 write(*,20,advance='no'); read *, x1 print *,' ' write(*,30,advance='no'); read *, nstep call Integral_Simpson(x0,x1,nstep,result) print *,' ' print *,' ' print *,' Value of integral: ', result print *,' ' print *,' ' 10 format(' X0 = ') 20 format(' X1 = ') 30 format(' Number of integration steps: ') END ! Given function to integrate real*8 Function FUNC(x) real*8 x FUNC = dsin(x) end !******************************************************* !* Integral of a function FUNC(X) by Simpson's method * !* --------------------------------------------------- * !* INPUTS: * !* a begin value of x variable * !* b end value of x variable * !* n number of integration steps * !* * !* OUTPUT: * !* res the integral of FUNC(X) from a to b * !* * !******************************************************* Subroutine Integral_Simpson(a, b, n, res) real*8 a,b,res integer n, i real*8 step, r step=(b-a)/2/n r=FUNC(a) res=(r+FUNC(b))/2 do i=1, 2*n-1 r=FUNC(a+i*step); if (MOD(i,2).ne.0) then res = res + r+r else res = res + r end if end do res = res * step*2.d0/3.d0 return end ! End of file tsimpson.f90