!  NUMINT   TRU   Demonstrates numerical integration.
!       type run
!       to run numerical demos:
!       call rdemo
!       call tdemo
!       call sdemo

def  f(x)=4/(1+x*x)
let a = 0          
let b = 1
let truevalue = pi
let maxiteration = 20

def riemann 
  let h = (b-a)/n
  let sum = 0
  let x = a
  for k = 0 to n-1
    let sum = sum + f(x)
    let x = x + h
  next k
  let riemann = h * sum
end def

def  trapazoid
  let h = (b-a)/n
  let sum = ( f(a) + f(b) ) / 2
  let x = a + h
  for k = 1 to n-1
    let sum = sum + f(x)
    let x = x + h
  next k
  let trapazoid = h * sum
end def

def  simpson 
  let h = (b-a)/n
  let evensum = 0
  let oddsum  = 0
  let x = a + h
  for k = 1 to n-1 step 2
    let oddsum  = oddsum  +  f(x) 
    let x = x + h
    let evensum = evensum + f(x)
    let x=  x + h
  next k
  let evensum = evensum - f(b)
  let simpson = (h/3) * (f(a) + 4*oddsum + 2*evensum + f(b))
end def

sub rdemo
clear
let format$="  ####  #.#######  +#.#^^^^  +#.#^^^^"
print       "              Riemann Sum           "
print
print       "   N     RIEMSUM    ERROR   ERROR/H"
print
let n = 1
for i = 1 to maxiteration
  let compvalue = riemann
  let error = compvalue-truevalue
  print using format$: n,compvalue,error,error/h
  let n = 2 * n
next i
line input a$
end sub

sub tdemo
clear
let format$="  ####  #.#######  +#.#^^^^  +#.#^^^^ "
print       "              Trapazoid Rule          "
print
print       "   N     TRAPSUM    ERROR   ERROR/H^2"
print
let n = 1
for i = 1 to maxiteration
  let compvalue = trapazoid
  let error = compvalue-truevalue
  print using format$: n,compvalue,error,error/(h*h)
  let n = 2 * n
next i
line input a$
end sub

sub sdemo
clear
let format$="####  #.###########  +#.#^^^^  +#.#^^^^"
print       "              Simpson Rule             "
print
print       "  N     SIMPSUM       ERROR   ERROR/H^4"
print
let n = 1
for i = 1 to maxiteration
  let compvalue = simpson
  let error = compvalue-truevalue
  print using format$: n,compvalue,error,error/(h*h*h*h)
  let n = 2 * n
next i
line input a$
end sub

sub macheps
  let k = 1
  let eps =1
  do   until (1+ eps)= 1
     print k,eps
     let k=k+1
     let eps = eps/2
  loop
end sub   

end