cdeagle Senior Member
Joined: 22/06/2014 Location: United StatesPosts: 261 |
Posted: 03:59am 13 Apr 2017 |
Copy link to clipboard |
Print this post |
|
This program demonstrates the procedures for calling several 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)
Here's a typical user interaction and the results produced by the software.
chebyshev approximations ------------------------ please input the maximum degree of the chebyshev approximation ? 20 please input the lower limit of the evaluation interval ? 1 please input the upper limit of the evaluation interval ? 2 please input the number of terms in the chebyshev approximation ? 15 please input the x argument for evaluation ? 1.5 chebyshev approximations ------------------------ lower limit of interval = 1 upper limit of interval = 2 function argument = 1.5 chebyshev function value = 0.56109093 analytic function value = 0.56109093 error = 1.11022e-15 chebyshev derivative = 7.52100208 analytic derivative = 7.52100208 error = 8.88178e-15 chebyshev integral = -0.23929577 analytic integral = -0.23929577 error = 2.77556e-16 maximum degree of approximation = 20 degree of chebyshev fit = 15 Here's the source code listing.
' 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
|