' program "demo_chebyshev.bas" April 13, 2017

' This program demonstrates the procedures for calling
' the Chebyshev subroutines. These subroutines can be
' used to approximate the integral, derivative, and
' function value of a user-defined analytic function.
' This same information can be determined for tabular
' data of the form y = f(x).

' This program demonstrates the use of the Chebyshev
' subroutines for evaluating information about

' y = x^2 * (x^2 - 2.0) * sin(x)

'''''''''''''''''''''''''''''''''''

option default float

option base 1

const pidiv2 = 0.5 * pi

' begin main program

print " "
print "chebyshev approximations"
print "------------------------"

do

print " "
print "please input the maximum degree of the chebyshev approximation"
input nval%

loop until (nval% > 0)

dim c(nval%), cder(nval%), cinq(nval%)

print " "
print "please input the lower limit of the evaluation interval"
input xl

print " "
print "please input the upper limit of the evaluation interval"
input xu

' compute chebyshev coefficients

cheby_fit(xl, xu, c(), nval%)

do

print " "
print "please input the number of terms in the chebyshev approximation"
input n%

loop until (n% > 0 and n% <= nval%)

do

print " "
print "please input the x argument for evaluation"
input x

loop until (x >= xl and x <= xu)

' evaluate chebyshev approximation of user-defined function

cheby_eval(xl, xu, c(), n%, x, fval1)

' compute analytic value of function

user_func(x, fval2)

' compute derivative chebyshev coefficients

cheby_der(xl, xu, c(), n%, cder())

' evaluate chebyshev derivative approximation

cheby_eval(xl, xu, cder(), n%, x, dfdx1)

' compute analytic value of derivative

dfdx(x, dfdx2)

' compute quadrature chebyshev coefficients

cheby_quad(xl, xu, c(), n%, cinq())

' evaluate chebyshev quadrature approximation

cheby_eval(xl, xu, cinq(), n%, x, fint1)

' compute analytic value of quadrature

intx(x, fx)

intx(xl, fl)

fint2 = fx - fl

' print results

print " "
print "chebyshev approximations"
print "------------------------"
print " "
print "lower limit of interval = "; xl
print "upper limit of interval = "; xu
print " "
print "function argument = "; x
print " "
print "chebyshev function value = "; str$(fval1, 0, 8)
print "analytic function value = "; str$(fval2, 0, 8)
print "error = "; abs(fval1 - fval2)
print " "
print "chebyshev derivative = "; str$(dfdx1, 0, 8)
print "analytic derivative = "; str$(dfdx2, 0, 8)
print "error = "; abs(dfdx1 - dfdx2)
print " "
print "chebyshev integral = "; str$(fint1, 0, 8)
print "analytic integral = "; str$(fint2, 0, 8)
print "error = "; abs(fint1 - fint2)
print " "
print "maximum degree of approximation = "; nval%
print " "
print "degree of chebyshev fit = "; n%
print " "

end

'''''''''''''''''''''''''''
'''''''''''''''''''''''''''

sub cheby_fit (a, b, c(), n%)

' compute chebshev coefficients subroutine

' input

' a = lower limit of evaluation interval
' b = upper limit of evaluation interval
' n% = maximum degree of chebyshev approximation

' output

' c() = array of chebyshev coefficients

' note: user-defined function coded as
' sub user_func(x, fx)

''''''''''''''''''''''''''''

local f(n%), bma, bpa, x, fx, fac, sum

bma = 0.5 * (b - a)

bpa = 0.5 * (b + a)

for k% = 1 to n%

x = cos(pi * (k% - 0.5) / n%) * bma + bpa

user_func(x, fx)

f(k%) = fx

next k%

fac = 2.0 / n%

for j% = 1 to n%

sum = 0.0

for k% = 1 to n%

sum = sum + f(k%) * cos((pi * (j% - 1)) * ((k% - 0.5) / n%))

next k%

c(j%) = fac * sum

next j%

end sub

'''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''

sub cheby_der (a, b, c(), n%, cder())

' chebyshev coefficients of the derivative subroutine

' input

' a = lower limit of evaluation interval
' b = upper limit of evaluation interval
' c() = array of chebyshev coefficients
' n% = size of c() array

' output

' cder() = array of chebyshev derivative coefficients

''''''''''''''''''''''''''''''''''''''''''''''''''''''

local con

cder(n%) = 0.0

cder(n% - 1) = 2.0 * (n% - 1) * c(n%)

if (n% >= 3) then

for j% = n% - 2 to 1 step -1

cder(j%) = cder(j% + 2) + 2.0 * j% * c(j% + 1)

next j%

end if

con = 2.0 / (b - a)

for j% = 1 to n%

cder(j%) = cder(j%) * con

next j%

end sub

''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''

sub cheby_eval (a, b, c(), m%, x, vcheb)

' evaluate chebshev coefficients subroutine

' input

' a = lower limit of evaluation interval
' b = upper limit of evaluation interval
' c() = array of chebshev coefficients
' m% = size of c() array
' x = argument

' output

' vcheb = function value at x

''''''''''''''''''''''''''''''

local d, dd, y, y2, sv

d = 0.0

dd = 0.0

y = (2.0 * x - a - b) / (b - a)

y2 = 2.0 * y

for j% = m% to 2 step -1

sv = d

d = y2 * d - dd + c(j%)

dd = sv

next j%

vcheb = y * d - dd + 0.5 * c(1)

end sub

''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''

sub cheby_quad (a, b, c(), n%, cinq())

' chebyshev quadrature coefficients subroutine

' input

' a = lower limit of evaluation interval
' b = upper limit of evaluation interval
' c() = array of chebyshev coefficients
' n% = size of c() array

' output

' cinq() = array of chebyshev quadrature coefficients

''''''''''''''''''''''''''''''''''''''''''''''''''''''

local con, sum, fac

con = 0.25 * (b - a)

sum = 0.0

fac = 1.0

for j% = 2 to n% - 1

cinq(j%) = con * (c(j% - 1) - c(j% + 1)) / (j% - 1)

sum = sum + fac * cinq(j%)

fac = -fac

next j%

cinq(n%) = con * c(n% - 1) / (n% - 1)

sum = sum + fac * cinq(n%)

cinq(1) = 2.0 * sum

end sub

''''''''''''''''''
''''''''''''''''''

sub user_func(x, fx)

' user-defined function subroutine

''''''''''''''''''''''''''''''''''

fx = (x * x) * (x * x - 2.0) * sin(x)

end sub

'''''''''''''
'''''''''''''

sub intx(x, fx)

' test function integral subroutine

'''''''''''''''''''''''''''''''''''

fx = 4.0 * x * (x * x - 7.0) * sin(x) - (x * x * x * x - 14.0 * x * x + 28.0) * cos(x)

end sub

'''''''''''''
'''''''''''''

sub dfdx(x, fx)

' test function derivative subroutine

'''''''''''''''''''''''''''''''''''''

fx = 4.0 * x * (x * x - 1.0) * sin(x) + x * x * (x * x - 2.0) * cos(x)

end sub