Chebyshev coefficients
Chebyshev coefficients are the basis of polynomial approximations of functions. Write a program to generate Chebyshev coefficients.
Calculate coefficients: cosine function, 10 coefficients, interval 0 1
C
C99. <lang C>#include <stdio.h>
- include <string.h>
- include <math.h>
- ifndef M_PI
- define M_PI 3.14159265358979323846
- endif
double test_func(double x) { return sin(cos(x)) * exp(-(x - 5)*(x - 5)/10); //return cos(x); }
// map x from range [min, max] to [min_to, max_to] double map(double x, double min_x, double max_x, double min_to, double max_to) { return (x - min_x)/(max_x - min_x)*(max_to - min_to) + min_to; }
void cheb_coef(double (*func)(double), int n, double min, double max, double *coef) { memset(coef, 0, sizeof(double) * n); for (int i = 0; i < n; i++) { double f = func(map(cos(M_PI*(i + .5f)/n), -1, 1, min, max))*2/n; for (int j = 0; j < n; j++) coef[j] += f*cos(M_PI*j*(i + .5f)/n); } }
// f(x) = sum_{k=0}^{n - 1} c_k T_k(x) - c_0/2 double cheb_approx(double x, int n, double min, double max, double *coef) { double a = 1, b = map(x, min, max, -1, 1), c; double res = coef[0]/2 + coef[1]*b;
x = 2*b; for (int i = 2; i < n; a = b, b = c, i++) // T_{n+1} = 2x T_n - T_{n-1} res += coef[i]*(c = x*b - a);
return res; }
int main(void) {
- define N 20
double c[N], min = 0, max = M_PI*3; cheb_coef(test_func, N, min, max, c);
printf("Coefficients:"); for (int i = 0; i < N; i++) printf(" %lg", c[i]);
puts("\n\nApproximation:\n x func(x) approx diff"); for (int i = 0; i <= 20; i++) { double x = map(i, 0, 20, min, max); double f = test_func(x); double approx = cheb_approx(x, N, min, max, c);
printf("% 10.8lf % 10.8lf % 10.8lf % 4.1le\n", x, f, approx, approx - f); }
return 0; }</lang>
J
From 'J for C Programmers: Calculating Chebyshev Coefficients [[1]] <lang J> chebft =: adverb define
f =. u 0.5 * (+/y) - (-/y) * 2 o. o. (0.5 + i. x) % x
(2 % x) * +/ f * 2 o. o. (0.5 + i. x) *"0 1 (i. x) % x
) </lang> Calculate coefficients: <lang J>
10 (2&o.) chebft 0 1
1.64717 _0.232299 _0.0537151 0.00245824 0.000282119 _7.72223e_6 _5.89856e_7 1.15214e_8 6.59629e_10 _1.00227e_11 </lang>
Perl
<lang perl>use constant PI => 3.141592653589793;
sub chebft {
my($func, $a, $b, $n) = @_; my($bma, $bpa) = ( 0.5*($b-$a), 0.5*($b+$a) );
my @pin = map { ($_ + 0.5) * (PI/$n) } 0..$n-1;
my @f = map { $func->( cos($_) * $bma + $bpa ) } @pin;
my @c = (0) x $n;
for my $j (0 .. $n-1) {
$c[$j] += $f[$_] * cos($j * $pin[$_]) for 0..$n-1;
$c[$j] *= (2.0/$n);
}
@c;
}
print "$_\n" for chebft(sub{cos($_[0])}, 0, 1, 10);</lang>
- Output:
1.64716947539031 -0.232299371615172 -0.0537151146220477 0.00245823526698163 0.000282119057433938 -7.72222915566001e-06 -5.89855645105608e-07 1.15214274787334e-08 6.59629917354465e-10 -1.00219943455215e-11
Perl 6
<lang perl6>sub chebft ( Code $func, Real $a, Real $b, Int $n ) {
my $bma = 0.5 * ( $b - $a ); my $bpa = 0.5 * ( $b + $a );
my @pi_n = ( (^$n).list »+» 0.5 ) »*» ( pi / $n );
my @f = ( @pi_n».cos »*» $bma »+» $bpa )».$func;
my @sums = map { [+] @f »*« ( @pi_n »*» $_ )».cos }, ^$n;
return @sums »*» ( 2 / $n );
}
say .fmt('%+13.7e') for chebft &cos, 0, 1, 10;</lang>
- Output:
+1.6471695e+00 -2.3229937e-01 -5.3715115e-02 +2.4582353e-03 +2.8211906e-04 -7.7222292e-06 -5.8985565e-07 +1.1521427e-08 +6.5962992e-10 -1.0021994e-11
Racket
<lang racket>#lang typed/racket (: chebft (Real Real Nonnegative-Integer (Real -> Real) -> (Vectorof Real))) (define (chebft a b n func)
(define b-a/2 (/ (- b a) 2)) (define b+a/2 (/ (+ b a) 2)) (define pi/n (/ pi n)) (define fac (/ 2 n))
(define f (for/vector : (Vectorof Real)
((k : Nonnegative-Integer (in-range n)))
(define y (cos (* pi/n (+ k 1/2))))
(func (+ (* y b-a/2) b+a/2))))
(for/vector : (Vectorof Real)
((j : Nonnegative-Integer (in-range n)))
(define s (for/sum : Real
((k : Nonnegative-Integer (in-range n)))
(* (vector-ref f k)
(cos (* pi/n j (+ k 1/2))))))
(* fac s)))
(module+ test
(chebft 0 1 10 cos))
- Tim Brown 2015</lang>
- Output:
'#(1.6471694753903137 -0.2322993716151719 -0.05371511462204768 0.0024582352669816343 0.0002821190574339161 -7.722229155637806e-006 -5.898556451056081e-007 1.1521427500937876e-008 6.596299173544651e-010 -1.0022016549982027e-011)
REXX
The following REXX program is a translation of the C program plus added optimizations. <lang rexx>/*REXX program calculates ten Chebyshev coefficients for the range 0 ──► 1.*/ /*Pafnuty Lvovich Chebysheff: Chebysheff [English transliteration] */ /*─────────────────────────── Chebyshov [ " " ] */ /*─────────────────────────── Tchebychev [French " ] */ /*─────────────────────────── Tchebysheff [ " " ] */ /*─────────────────────────── Tschebyschow [German " ] */ /*─────────────────────────── Tschebyschev [ " " ] */ /*─────────────────────────── Tschebyschef [ " " ] */ /*─────────────────────────── Tschebyscheff [ " " ] */ numeric digits length(pi()) /*DIGITS default is nine, but use 63. */ parse arg a b n . /*obtain optional arguments from the CL*/ if a== | a==',' then a=0 /*A not specified? Then use default.*/ if b== | b==',' then b=1 /*B " " " " " */ if n== | n==',' then n=10 /*N " " " " " */ bma =(b-a)/2 /*calculate one─half of the difference.*/ bpa =(b+a)/2 /* " " " " sum. */ pi@n=pi/n /*calculate a handy─dandy value. */
do k=0 for n
f.k=cos(cos(pi@n * (k+.5)) * bma + bpa)
end /*k*/
fac=2/n /*calculate another handy─dandy value. */
do j=0 for n; $=0
z=pi@n * j
do m=0 for n
$=$ + f.m * cos(z * (m+.5))
end /*m*/
cheby.j=fac*$ /*uses 63 digs, but only shows 30 ►───┐*/
say right(j,length(n)+3) ' Chebyshev coefficient is:', /*│*/
left(,cheby.j>=0) format(cheby.j,,30) /* ◄───────────┘*/
end /*j*/
exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ cos: procedure; parse arg x; numeric digits digits()+10; numeric fuzz 5
z=1; _=1; x=r2r(x); a=abs(x); q=x*x
if a=pi then return -1; if a=pi*.5 | a=pi*2 then return 0; pi3=pi/3
if a=pi3 then return .5; if a=2*pi3 then return -.5
do ?=2 by 2 until p=z; p=z; _=-_*q/(?*(?-1)); z=z+_; end /*?*/
return z
/*────────────────────────────────────────────────────────────────────────────*/ pi: pi=3.1415926535897932384626433832795028841971693993751058209749446;return pi /*────────────────────────────────────────────────────────────────────────────*/ r2r: return arg(1) // (pi()*2) /*normalize radians ───► a unit circle.*/</lang> output when using the default inputs:
0 Chebyshev coefficient is: 1.647169475390313686961473816798
1 Chebyshev coefficient is: -0.232299371615171942121038341178
2 Chebyshev coefficient is: -0.053715114622047555071596203933
3 Chebyshev coefficient is: 0.002458235266981479866768882753
4 Chebyshev coefficient is: 0.000282119057434005702410217295
5 Chebyshev coefficient is: -0.000007722229155810577892832847
6 Chebyshev coefficient is: -5.898556452177103343296676960522E-7
7 Chebyshev coefficient is: 1.152142733310315857327524390711E-8
8 Chebyshev coefficient is: 6.596300035120132380676859918562E-10
9 Chebyshev coefficient is: -1.002259170944625675156620531665E-11
Sidef
<lang ruby>func chebft (callback, a, b, n) {
var bma = (0.5 * b-a); var bpa = (0.5 * b+a);
var pi_n = ((0..(n-1) »+» 0.5) »*» (Math::PI / n));
var f = (pi_n »cos»() »*» bma »+» bpa «call« callback);
var sums = (0..(n-1) «run« {|i| f »*« ((pi_n »*» i) »cos»()) «+» });
sums »*» (2/n);
}
chebft(func(v){v.cos}, 0, 1, 10).each { |v|
say ("%+.10e" % v);
};</lang>
- Output:
+1.6471694754e+00 -2.3229937162e-01 -5.3715114622e-02 +2.4582352670e-03 +2.8211905743e-04 -7.7222291558e-06 -5.8985564522e-07 +1.1521427333e-08 +6.5963000351e-10 -1.0022591709e-11
zkl
<lang zkl>var [const] PI=(1.0).pi; fcn chebft(a,b,n,func){
bma,bpa,fac := 0.5*(b - a), 0.5*(b + a), 2.0/n;
f:=n.pump(List,'wrap(k){ (PI*(0.5 + k)/n).cos():func(_*bma + bpa) });
n.pump(List,'wrap(j){
fac*n.reduce('wrap(sum,k){ sum + f[k]*(PI*j*(0.5 + k)/n).cos() },0.0);
})
} chebft(0.0,1.0,10,fcn(x){ x.cos() }).enumerate().concat("\n").println();</lang>
- Output:
L(0,1.64717) L(1,-0.232299) L(2,-0.0537151) L(3,0.00245824) L(4,0.000282119) L(5,-7.72223e-06) L(6,-5.89856e-07) L(7,1.15214e-08) L(8,6.5963e-10) L(9,-1.00219e-11)