Chebyshev coefficients: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Perl}}: formatted output) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
Line 13: | Line 13: | ||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang=11l>F test_func(Float x) |
||
R cos(x) |
R cos(x) |
||
Line 56: | Line 56: | ||
V f = test_func(x) |
V f = test_func(x) |
||
V approx = cheb_approx(x, n, minv, maxv, c) |
V approx = cheb_approx(x, n, minv, maxv, c) |
||
print(‘#.3 #.10 #.10 #.’.format(x, f, approx, format_float_exp(approx - f, 2, 9)))</ |
print(‘#.3 #.10 #.10 #.’.format(x, f, approx, format_float_exp(approx - f, 2, 9)))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 102: | Line 102: | ||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
Given the limitations of the language, only 8 coefficients are calculated |
Given the limitations of the language, only 8 coefficients are calculated |
||
< |
<syntaxhighlight lang=BASIC256>a = 0: b = 1: n = 8 |
||
dim cheby(n) |
dim cheby(n) |
||
dim coef(n) |
dim coef(n) |
||
Line 118: | Line 118: | ||
print i; " : "; cheby[i] |
print i; " : "; cheby[i] |
||
next i |
next i |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>0 : 1.64716947539 |
<pre>0 : 1.64716947539 |
||
Line 133: | Line 133: | ||
{{works with|QuickBasic|4.5}} |
{{works with|QuickBasic|4.5}} |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
< |
<syntaxhighlight lang=qbasic>pi = 4 * ATN(1) |
||
a = 0: b = 1: n = 10 |
a = 0: b = 1: n = 10 |
||
DIM cheby!(n) |
DIM cheby!(n) |
||
Line 150: | Line 150: | ||
PRINT USING " # : ##.#####################"; i; cheby(i) |
PRINT USING " # : ##.#####################"; i; cheby(i) |
||
NEXT i |
NEXT i |
||
END</ |
END</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 166: | Line 166: | ||
==={{header|FreeBASIC}}=== |
==={{header|FreeBASIC}}=== |
||
< |
<syntaxhighlight lang=freebasic>Const pi As Double = 4 * Atn(1) |
||
Dim As Double i, w, j |
Dim As Double i, w, j |
||
Dim As Double a = 0, b = 1, n = 10 |
Dim As Double a = 0, b = 1, n = 10 |
||
Line 183: | Line 183: | ||
Print i; " : "; cheby(i) |
Print i; " : "; cheby(i) |
||
Next i |
Next i |
||
Sleep</ |
Sleep</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> 0 : 1.647169475390314 |
<pre> 0 : 1.647169475390314 |
||
Line 198: | Line 198: | ||
==={{header|Yabasic}}=== |
==={{header|Yabasic}}=== |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
< |
<syntaxhighlight lang=yabasic>a = 0: b = 1: n = 10 |
||
dim cheby(n) |
dim cheby(n) |
||
dim coef(n) |
dim coef(n) |
||
Line 214: | Line 214: | ||
print i, " : ", cheby(i) |
print i, " : ", cheby(i) |
||
next i |
next i |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>0 : 1.64717 |
<pre>0 : 1.64717 |
||
Line 230: | Line 230: | ||
=={{header|C}}== |
=={{header|C}}== |
||
C99. |
C99. |
||
< |
<syntaxhighlight lang=C>#include <stdio.h> |
||
#include <string.h> |
#include <string.h> |
||
#include <math.h> |
#include <math.h> |
||
Line 296: | Line 296: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang=csharp>using System; |
||
using System.Collections.Generic; |
using System.Collections.Generic; |
||
using System.Linq; |
using System.Linq; |
||
Line 408: | Line 408: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 452: | Line 452: | ||
The wrapper class ChebyshevApprox_ supports very terse user code. |
The wrapper class ChebyshevApprox_ supports very terse user code. |
||
< |
<syntaxhighlight lang=CPP> |
||
#include <iostream> |
#include <iostream> |
||
#include <iomanip> |
#include <iomanip> |
||
Line 554: | Line 554: | ||
} |
} |
||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|D}}== |
=={{header|D}}== |
||
This imperative code retains some of the style of the original C version. |
This imperative code retains some of the style of the original C version. |
||
< |
<syntaxhighlight lang=d>import std.math: PI, cos; |
||
/// Map x from range [min, max] to [min_to, max_to]. |
/// Map x from range [min, max] to [min_to, max_to]. |
||
Line 620: | Line 620: | ||
writefln("%1.3f % 10.10f % 10.10f % 4.2e", x, f, approx, approx - f); |
writefln("%1.3f % 10.10f % 10.10f % 4.2e", x, f, approx, approx - f); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 700: | Line 700: | ||
=={{header|EasyLang}}== |
=={{header|EasyLang}}== |
||
<lang>numfmt 0 5 |
<syntaxhighlight lang=text>numfmt 0 5 |
||
a = 0 |
a = 0 |
||
b = 1 |
b = 1 |
||
Line 716: | Line 716: | ||
cheby[i] = w * 2 / n |
cheby[i] = w * 2 / n |
||
print cheby[i] |
print cheby[i] |
||
.</ |
.</syntaxhighlight> |
||
=={{header|Go}}== |
=={{header|Go}}== |
||
Line 724: | Line 724: | ||
Two variances here from the WP presentation and most mathematical presentations follow other examples on this page and so keep output directly comparable. One variance is that the Kronecker delta factor is dropped, which has the effect of doubling the first coefficient. This simplifies both coefficient generation and polynomial evaluation. A further variance is that there is no scaling for the range of function values. The result is that coefficients are not necessarily bounded by 1 (2 for the first coefficient) but by the maximum function value over the argument range from min to max (or twice that for the first coefficient.) |
Two variances here from the WP presentation and most mathematical presentations follow other examples on this page and so keep output directly comparable. One variance is that the Kronecker delta factor is dropped, which has the effect of doubling the first coefficient. This simplifies both coefficient generation and polynomial evaluation. A further variance is that there is no scaling for the range of function values. The result is that coefficients are not necessarily bounded by 1 (2 for the first coefficient) but by the maximum function value over the argument range from min to max (or twice that for the first coefficient.) |
||
< |
<syntaxhighlight lang=go>package main |
||
import ( |
import ( |
||
Line 787: | Line 787: | ||
} |
} |
||
return x1*t - s + .5*c.c[0] |
return x1*t - s + .5*c.c[0] |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 818: | Line 818: | ||
=={{header|Groovy}}== |
=={{header|Groovy}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang=groovy>class ChebyshevCoefficients { |
||
static double map(double x, double min_x, double max_x, double min_to, double max_to) { |
static 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 |
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to |
||
Line 846: | Line 846: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 862: | Line 862: | ||
=={{header|J}}== |
=={{header|J}}== |
||
From 'J for C Programmers: Calculating Chebyshev Coefficients [[http://www.jsoftware.com/learning/a_first_look_at_j_programs.htm#_Toc191734318]] |
From 'J for C Programmers: Calculating Chebyshev Coefficients [[http://www.jsoftware.com/learning/a_first_look_at_j_programs.htm#_Toc191734318]] |
||
<syntaxhighlight lang=J> |
|||
<lang J> |
|||
chebft =: adverb define |
chebft =: adverb define |
||
: |
: |
||
Line 868: | Line 868: | ||
(2 % x) * +/ f * 2 o. o. (0.5 + i. x) *"0 1 (i. x) % x |
(2 % x) * +/ f * 2 o. o. (0.5 + i. x) *"0 1 (i. x) % x |
||
) |
) |
||
</syntaxhighlight> |
|||
</lang> |
|||
Calculate coefficients: |
Calculate coefficients: |
||
<syntaxhighlight lang=J> |
|||
<lang J> |
|||
10 (2&o.) chebft 0 1 |
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 |
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 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Java}}== |
=={{header|Java}}== |
||
Partial translation of [[Chebyshev_coefficients#C|C]] via [[Chebyshev_coefficients#D|D]] |
Partial translation of [[Chebyshev_coefficients#C|C]] via [[Chebyshev_coefficients#D|D]] |
||
{{works with|Java|8}} |
{{works with|Java|8}} |
||
< |
<syntaxhighlight lang=java>import static java.lang.Math.*; |
||
import java.util.function.Function; |
import java.util.function.Function; |
||
Line 914: | Line 914: | ||
System.out.println(d); |
System.out.println(d); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 934: | Line 934: | ||
'''Preliminaries''' |
'''Preliminaries''' |
||
< |
<syntaxhighlight lang=jq>def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .; |
||
def rpad($len; $fill): tostring | ($len - length) as $l | . + ($fill * $l)[:$l]; |
def rpad($len; $fill): tostring | ($len - length) as $l | . + ($fill * $l)[:$l]; |
||
Line 950: | Line 950: | ||
| ((if $ix then $s[0:$ix] else $s end) | lpad) + "." + |
| ((if $ix then $s[0:$ix] else $s end) | lpad) + "." + |
||
(if $ix then $s[$ix+1:] | .[0:right] else "" end) |
(if $ix then $s[$ix+1:] | .[0:right] else "" end) |
||
end;</ |
end;</syntaxhighlight> |
||
'''Chebyshev Coefficients''' |
'''Chebyshev Coefficients''' |
||
< |
<syntaxhighlight lang=jq>def mapRange($x; $min; $max; $minTo; $maxTo): |
||
(($x - $min)/($max - $min))*($maxTo - $minTo) + $minTo; |
(($x - $min)/($max - $min))*($maxTo - $minTo) + $minTo; |
||
Line 993: | Line 993: | ||
| join(" ") ); |
| join(" ") ); |
||
task</ |
task</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,037: | Line 1,037: | ||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang=julia>mutable struct Cheb |
||
c::Vector{Float64} |
c::Vector{Float64} |
||
min::Float64 |
min::Float64 |
||
Line 1,086: | Line 1,086: | ||
approx = evaluate(c, x) |
approx = evaluate(c, x) |
||
@printf("%.1f %12.8f %12.8f % .3e\n", x, computed, approx, computed - approx) |
@printf("%.1f %12.8f %12.8f % .3e\n", x, computed, approx, computed - approx) |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,116: | Line 1,116: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang=scala>// version 1.1.2 |
||
typealias DFunc = (Double) -> Double |
typealias DFunc = (Double) -> Double |
||
Line 1,164: | Line 1,164: | ||
System.out.printf("%1.3f %1.8f %1.8f % 4.1e\n", x, f, approx, approx - f) |
System.out.printf("%1.3f %1.8f %1.8f % 4.1e\n", x, f, approx, approx - f) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,207: | Line 1,207: | ||
=={{header|Lua}}== |
=={{header|Lua}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang=lua>function map(x, min_x, max_x, min_to, max_to) |
||
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to |
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to |
||
end |
end |
||
Line 1,244: | Line 1,244: | ||
main() |
main() |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 1,260: | Line 1,260: | ||
=={{header|Microsoft Small Basic}}== |
=={{header|Microsoft Small Basic}}== |
||
{{trans|Perl}} |
{{trans|Perl}} |
||
< |
<syntaxhighlight lang=smallbasic>' N Chebyshev coefficients for the range 0 to 1 - 18/07/2018 |
||
pi=Math.pi |
pi=Math.pi |
||
a=0 |
a=0 |
||
Line 1,279: | Line 1,279: | ||
EndIf |
EndIf |
||
TextWindow.WriteLine(i+" : "+t+cheby[i]) |
TextWindow.WriteLine(i+" : "+t+cheby[i]) |
||
EndFor</ |
EndFor</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,296: | Line 1,296: | ||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang=Nim>import lenientops, math, strformat, sugar |
||
type Cheb = object |
type Cheb = object |
||
Line 1,344: | Line 1,344: | ||
let computed = fn(x) |
let computed = fn(x) |
||
let approx = cheb.eval(x) |
let approx = cheb.eval(x) |
||
echo &"{x:.1f} {computed:12.8f} {approx:12.8f} {computed-approx: .3e}"</ |
echo &"{x:.1f} {computed:12.8f} {approx:12.8f} {computed-approx: .3e}"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,375: | Line 1,375: | ||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang=perl>use constant PI => 2 * atan2(1, 0); |
||
sub chebft { |
sub chebft { |
||
Line 1,391: | Line 1,391: | ||
} |
} |
||
printf "%+13.7e\n", $_ for chebft(sub{cos($_[0])}, 0, 1, 10);</ |
printf "%+13.7e\n", $_ for chebft(sub{cos($_[0])}, 0, 1, 10);</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>+1.6471695e+00 |
<pre>+1.6471695e+00 |
||
Line 1,406: | Line 1,406: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<!--< |
<!--<syntaxhighlight lang=Phix>(phixonline)--> |
||
<span style="color: #008080;">function</span> <span style="color: #000000;">Cheb</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">cmin</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cmax</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">ncoeff</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nnodes</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">function</span> <span style="color: #000000;">Cheb</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">cmin</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cmax</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">ncoeff</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nnodes</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ncoeff</span><span style="color: #0000FF;">),</span> |
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ncoeff</span><span style="color: #0000FF;">),</span> |
||
Line 1,450: | Line 1,450: | ||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%.1f %12.8f %12.8f %10.3e\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">calc</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">est</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">calc</span><span style="color: #0000FF;">-</span><span style="color: #000000;">est</span><span style="color: #0000FF;">})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%.1f %12.8f %12.8f %10.3e\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">calc</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">est</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">calc</span><span style="color: #0000FF;">-</span><span style="color: #000000;">est</span><span style="color: #0000FF;">})</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,481: | Line 1,481: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang=python>import math |
||
def test_func(x): |
def test_func(x): |
||
Line 1,532: | Line 1,532: | ||
return None |
return None |
||
main()</ |
main()</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 1,574: | Line 1,574: | ||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang=racket>#lang typed/racket |
||
(: chebft (Real Real Nonnegative-Integer (Real -> Real) -> (Vectorof Real))) |
(: chebft (Real Real Nonnegative-Integer (Real -> Real) -> (Vectorof Real))) |
||
(define (chebft a b n func) |
(define (chebft a b n func) |
||
Line 1,597: | Line 1,597: | ||
(module+ test |
(module+ test |
||
(chebft 0 1 10 cos)) |
(chebft 0 1 10 cos)) |
||
;; Tim Brown 2015</ |
;; Tim Brown 2015</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,615: | Line 1,615: | ||
{{trans|C}} |
{{trans|C}} |
||
<lang |
<syntaxhighlight lang=raku line>sub chebft ( Code $func, Real \a, Real \b, Int \n ) { |
||
my \bma = ½ × (b - a); |
my \bma = ½ × (b - a); |
||
Line 1,627: | Line 1,627: | ||
} |
} |
||
say chebft(&cos, 0, 1, 10)».fmt: '%+13.7e';</ |
say chebft(&cos, 0, 1, 10)».fmt: '%+13.7e';</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,656: | Line 1,656: | ||
The numeric precision is dependent on the number of decimal digits specified in the value of '''pi'''. |
The numeric precision is dependent on the number of decimal digits specified in the value of '''pi'''. |
||
< |
<syntaxhighlight lang=rexx>/*REXX program calculates N Chebyshev coefficients for the range 0 ──► 1 (inclusive)*/ |
||
numeric digits length( pi() ) - length(.) /*DIGITS default is nine, but use 71. */ |
numeric digits length( pi() ) - length(.) /*DIGITS default is nine, but use 71. */ |
||
parse arg a b N . /*obtain optional arguments from the CL*/ |
parse arg a b N . /*obtain optional arguments from the CL*/ |
||
Line 1,685: | Line 1,685: | ||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164;return pi |
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164;return pi |
||
r2r: return arg(1) // (pi() * 2) /*normalize radians ───► a unit circle.*/</ |
r2r: return arg(1) // (pi() * 2) /*normalize radians ───► a unit circle.*/</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 1,724: | Line 1,724: | ||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
< |
<syntaxhighlight lang=ruby>def mapp(x, min_x, max_x, min_to, max_to) |
||
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to |
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to |
||
end |
end |
||
Line 1,752: | Line 1,752: | ||
end |
end |
||
main()</ |
main()</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 1,768: | Line 1,768: | ||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/DqRNe2A/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/M5Ye6h8ZRkmTCNzexUh3uw Scastie (remote JVM)]. |
{{Out}}Best seen running in your browser either by [https://scalafiddle.io/sf/DqRNe2A/0 ScalaFiddle (ES aka JavaScript, non JVM)] or [https://scastie.scala-lang.org/M5Ye6h8ZRkmTCNzexUh3uw Scastie (remote JVM)]. |
||
< |
<syntaxhighlight lang=Scala>import scala.math.{Pi, cos} |
||
object ChebyshevCoefficients extends App { |
object ChebyshevCoefficients extends App { |
||
Line 1,798: | Line 1,798: | ||
c.foreach(d => println(f"$d%23.16e")) |
c.foreach(d => println(f"$d%23.16e")) |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang=ruby>func chebft (callback, a, b, n) { |
||
var bma = (0.5 * b-a); |
var bma = (0.5 * b-a); |
||
Line 1,816: | Line 1,816: | ||
chebft(func(v){v.cos}, 0, 1, 10).each { |v| |
chebft(func(v){v.cos}, 0, 1, 10).each { |v| |
||
say ("%+.10e" % v); |
say ("%+.10e" % v); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,836: | Line 1,836: | ||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
< |
<syntaxhighlight lang=swift>import Foundation |
||
typealias DFunc = (Double) -> Double |
typealias DFunc = (Double) -> Double |
||
Line 1,893: | Line 1,893: | ||
print(String(format: "%1.3f %1.8f %1.8f % 4.1e", x, f, approx, approx - f)) |
print(String(format: "%1.3f %1.8f %1.8f % 4.1e", x, f, approx, approx - f)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,936: | Line 1,936: | ||
{{trans|Microsoft Small Basic}} |
{{trans|Microsoft Small Basic}} |
||
To run in console mode with cscript. |
To run in console mode with cscript. |
||
< |
<syntaxhighlight lang=vb>' N Chebyshev coefficients for the range 0 to 1 |
||
Dim coef(10),cheby(10) |
Dim coef(10),cheby(10) |
||
pi=4*Atn(1) |
pi=4*Atn(1) |
||
Line 1,951: | Line 1,951: | ||
If cheby(i)<=0 Then t="" Else t=" " |
If cheby(i)<=0 Then t="" Else t=" " |
||
WScript.StdOut.WriteLine i&" : "&t&cheby(i) |
WScript.StdOut.WriteLine i&" : "&t&cheby(i) |
||
Next</ |
Next</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,968: | Line 1,968: | ||
=={{header|Visual Basic .NET}}== |
=={{header|Visual Basic .NET}}== |
||
{{trans|C#}} |
{{trans|C#}} |
||
< |
<syntaxhighlight lang=vbnet>Module Module1 |
||
Structure ChebyshevApprox |
Structure ChebyshevApprox |
||
Line 2,074: | Line 2,074: | ||
End Sub |
End Sub |
||
End Module</ |
End Module</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Coefficients: |
<pre>Coefficients: |
||
Line 2,114: | Line 2,114: | ||
{{trans|Kotlin}} |
{{trans|Kotlin}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang=ecmascript>import "/fmt" for Fmt |
||
var mapRange = Fn.new { |x, min, max, minTo, maxTo| (x - min)/(max - min)*(maxTo - minTo) + minTo } |
var mapRange = Fn.new { |x, min, max, minTo, maxTo| (x - min)/(max - min)*(maxTo - minTo) + minTo } |
||
Line 2,161: | Line 2,161: | ||
diffStr = (diff >= 0) ? " " + diffStr[0..3] : diffStr[0..4] |
diffStr = (diff >= 0) ? " " + diffStr[0..3] : diffStr[0..4] |
||
Fmt.print("$1.3f $1.8f $1.8f $s", x, f, approx, diffStr + e) |
Fmt.print("$1.3f $1.8f $1.8f $s", x, f, approx, diffStr + e) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,204: | Line 2,204: | ||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|C}}{{trans|Perl}} |
{{trans|C}}{{trans|Perl}} |
||
< |
<syntaxhighlight lang=zkl>var [const] PI=(1.0).pi; |
||
fcn chebft(a,b,n,func){ |
fcn chebft(a,b,n,func){ |
||
bma,bpa,fac := 0.5*(b - a), 0.5*(b + a), 2.0/n; |
bma,bpa,fac := 0.5*(b - a), 0.5*(b + a), 2.0/n; |
||
Line 2,212: | Line 2,212: | ||
}) |
}) |
||
} |
} |
||
chebft(0.0,1.0,10,fcn(x){ x.cos() }).enumerate().concat("\n").println();</ |
chebft(0.0,1.0,10,fcn(x){ x.cos() }).enumerate().concat("\n").println();</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Revision as of 19:22, 26 August 2022
Chebyshev coefficients are the basis of polynomial approximations of functions.
- Task
Write a program to generate Chebyshev coefficients.
Calculate coefficients: cosine function, 10 coefficients, interval 0 1
11l
F test_func(Float x)
R cos(x)
F mapper(x, min_x, max_x, min_to, max_to)
R (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to
F cheb_coef(func, n, min, max)
V coef = [0.0] * n
L(i) 0 .< n
V f = func(mapper(cos(math:pi * (i + 0.5) / n), -1, 1, min, max)) * 2 / n
L(j) 0 .< n
coef[j] += f * cos(math:pi * j * (i + 0.5) / n)
R coef
F cheb_approx(=x, n, min, max, coef)
V a = 1.0
V b = mapper(x, min, max, -1, 1)
V res = coef[0] / 2 + coef[1] * b
x = 2 * b
V i = 2
L i < n
V c = x * b - a
res = res + coef[i] * c
(a, b) = (b, c)
i++
R res
V n = 10
V minv = 0
V maxv = 1
V c = cheb_coef(test_func, n, minv, maxv)
print(‘Coefficients:’)
L(i) 0 .< n
print(c[i])
print("\n\nApproximation:\n x func(x) approx diff")
L(i) 20
V x = mapper(i, 0.0, 20.0, minv, maxv)
V f = test_func(x)
V approx = cheb_approx(x, n, minv, maxv, c)
print(‘#.3 #.10 #.10 #.’.format(x, f, approx, format_float_exp(approx - f, 2, 9)))
- Output:
Coefficients: 1.64717 -0.232299 -0.0537151 0.00245824 0.000282119 -7.72223e-06 -5.89856e-07 1.15214e-08 6.5963e-10 -1.00219e-11 Approximation: x func(x) approx diff 0.000 1.0000000000 1.0000000000 4.68e-13 0.050 0.9987502604 0.9987502604 -9.36e-14 0.100 0.9950041653 0.9950041653 4.62e-13 0.150 0.9887710779 0.9887710779 -4.73e-14 0.200 0.9800665778 0.9800665778 -4.60e-13 0.250 0.9689124217 0.9689124217 -2.32e-13 0.300 0.9553364891 0.9553364891 2.62e-13 0.350 0.9393727128 0.9393727128 4.61e-13 0.400 0.9210609940 0.9210609940 1.98e-13 0.450 0.9004471024 0.9004471024 -2.47e-13 0.500 0.8775825619 0.8775825619 -4.58e-13 0.550 0.8525245221 0.8525245221 -2.46e-13 0.600 0.8253356149 0.8253356149 1.96e-13 0.650 0.7960837985 0.7960837985 4.53e-13 0.700 0.7648421873 0.7648421873 2.54e-13 0.750 0.7316888689 0.7316888689 -2.28e-13 0.800 0.6967067093 0.6967067093 -4.47e-13 0.850 0.6599831459 0.6599831459 -4.37e-14 0.900 0.6216099683 0.6216099683 4.46e-13 0.950 0.5816830895 0.5816830895 -8.98e-14
BASIC
BASIC256
Given the limitations of the language, only 8 coefficients are calculated
a = 0: b = 1: n = 8
dim cheby(n)
dim coef(n)
for i = 0 to n-1
coef[i] = cos(cos(pi/n*(i+1/2))*(b-a)/2+(b+a)/2)
next i
for i = 0 to n-1
w = 0
for j = 0 to n-1
w += coef[j] * cos(pi/n*i*(j+1/2))
next j
cheby[i] = w*2/n
print i; " : "; cheby[i]
next i
end
- Output:
0 : 1.64716947539 1 : -0.23229937162 2 : -0.05371511462 3 : 0.00245823527 4 : 0.00028211906 5 : -0.00000772223 6 : -5.89855645106e-07 7 : 1.15214275009e-08
QBasic
pi = 4 * ATN(1)
a = 0: b = 1: n = 10
DIM cheby!(n)
DIM coef!(n)
FOR i = 0 TO n - 1
coef(i) = COS(COS(pi / n * (i + 1 / 2)) * (b - a) / 2 + (b + a) / 2)
NEXT i
FOR i = 0 TO n - 1
w = 0
FOR j = 0 TO n - 1
w = w + coef(j) * COS(pi / n * i * (j + 1 / 2))
NEXT j
cheby(i) = w * 2 / n
PRINT USING " # : ##.#####################"; i; cheby(i)
NEXT i
END
- Output:
0 : 1.647169470787048000000 1 : -0.232299402356147800000 2 : -0.053715050220489500000 3 : 0.002458173315972090000 4 : 0.000282166845863685000 5 : -0.000007787576578266453 6 : -0.000000536595905487047 7 : 0.000000053614126471757 8 : 0.000000079823998078155 9 : -0.000000070922546058227
FreeBASIC
Const pi As Double = 4 * Atn(1)
Dim As Double i, w, j
Dim As Double a = 0, b = 1, n = 10
Dim As Double cheby(10), coef(10)
For i = 0 To n-1
coef(i) = Cos(Cos(pi/n*(i+1/2))*(b-a)/2+(b+a)/2)
Next i
For i = 0 To n-1
w = 0
For j = 0 To n-1
w += coef(j) * Cos(pi/n*i*(j+1/2))
Next j
cheby(i) = w*2/n
Print i; " : "; cheby(i)
Next i
Sleep
- Output:
0 : 1.647169475390314 1 : -0.2322993716151719 2 : -0.05371511462204768 3 : 0.002458235266981634 4 : 0.0002821190574339161 5 : -7.7222291556156e-006 6 : -5.898556451056081e-007 7 : 1.152142750093788e-008 8 : 6.596299062522348e-010 9 : -1.002201654998203e-011
Yabasic
a = 0: b = 1: n = 10
dim cheby(n)
dim coef(n)
for i = 0 to n-1
coef(i) = cos(cos(pi/n*(i+1/2))*(b-a)/2+(b+a)/2)
next i
for i = 0 to n-1
w = 0
for j = 0 to n-1
w = w + coef(j) * cos(pi/n*i*(j+1/2))
next j
cheby(i) = w*2/n
print i, " : ", cheby(i)
next i
end
- Output:
0 : 1.64717 1 : -0.232299 2 : -0.0537151 3 : 0.00245824 4 : 0.000282119 5 : -7.72223e-06 6 : -5.89856e-07 7 : 1.15214e-08 8 : 6.5963e-10 9 : -1.0022e-11
C
C99.
#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
// Note that n >= 2 is assumed; probably should check for that, however silly it is.
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 10
double c[N], min = 0, max = 1;
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;
}
C#
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace Chebyshev {
class Program {
struct ChebyshevApprox {
public readonly List<double> coeffs;
public readonly Tuple<double, double> domain;
public ChebyshevApprox(Func<double, double> func, int n, Tuple<double, double> domain) {
coeffs = ChebCoef(func, n, domain);
this.domain = domain;
}
public double Call(double x) {
return ChebEval(coeffs, domain, x);
}
}
static double AffineRemap(Tuple<double, double> from, double x, Tuple<double, double> to) {
return to.Item1 + (x - from.Item1) * (to.Item2 - to.Item1) / (from.Item2 - from.Item1);
}
static List<double> ChebCoef(List<double> fVals) {
int n = fVals.Count;
double theta = Math.PI / n;
List<double> retval = new List<double>();
for (int i = 0; i < n; i++) {
retval.Add(0.0);
}
for (int ii = 0; ii < n; ii++) {
double f = fVals[ii] * 2.0 / n;
double phi = (ii + 0.5) * theta;
double c1 = Math.Cos(phi);
double s1 = Math.Sin(phi);
double c = 1.0;
double s = 0.0;
for (int j = 0; j < n; j++) {
retval[j] += f * c;
// update c -> cos(j*phi) for next value of j
double cNext = c * c1 - s * s1;
s = c * s1 + s * c1;
c = cNext;
}
}
return retval;
}
static List<double> ChebCoef(Func<double, double> func, int n, Tuple<double, double> domain) {
double remap(double x) {
return AffineRemap(new Tuple<double, double>(-1.0, 1.0), x, domain);
}
double theta = Math.PI / n;
List<double> fVals = new List<double>();
for (int i = 0; i < n; i++) {
fVals.Add(0.0);
}
for (int ii = 0; ii < n; ii++) {
fVals[ii] = func(remap(Math.Cos((ii + 0.5) * theta)));
}
return ChebCoef(fVals);
}
static double ChebEval(List<double> coef, double x) {
double a = 1.0;
double b = x;
double c;
double retval = 0.5 * coef[0] + b * coef[1];
var it = coef.GetEnumerator();
it.MoveNext();
it.MoveNext();
while (it.MoveNext()) {
double pc = it.Current;
c = 2.0 * b * x - a;
retval += pc * c;
a = b;
b = c;
}
return retval;
}
static double ChebEval(List<double> coef, Tuple<double, double> domain, double x) {
return ChebEval(coef, AffineRemap(domain, x, new Tuple<double, double>(-1.0, 1.0)));
}
static void Main() {
const int N = 10;
ChebyshevApprox fApprox = new ChebyshevApprox(Math.Cos, N, new Tuple<double, double>(0.0, 1.0));
Console.WriteLine("Coefficients: ");
foreach (var c in fApprox.coeffs) {
Console.WriteLine("\t{0: 0.00000000000000;-0.00000000000000;zero}", c);
}
Console.WriteLine("\nApproximation:\n x func(x) approx diff");
const int nX = 20;
const int min = 0;
const int max = 1;
for (int i = 0; i < nX; i++) {
double x = AffineRemap(new Tuple<double, double>(0, nX), i, new Tuple<double, double>(min, max));
double f = Math.Cos(x);
double approx = fApprox.Call(x);
Console.WriteLine("{0:0.000} {1:0.00000000000000} {2:0.00000000000000} {3:E}", x, f, approx, approx - f);
}
}
}
}
- Output:
Coefficients: 1.64716947539031 -0.23229937161517 -0.05371511462205 0.00245823526698 0.00028211905743 -0.00000772222916 -0.00000058985565 0.00000001152143 0.00000000065963 -0.00000000001002 Approximation: x func(x) approx diff 0.000 1.00000000000000 1.00000000000047 4.689582E-013 0.050 0.99875026039497 0.99875026039487 -9.370282E-014 0.100 0.99500416527803 0.99500416527849 4.622969E-013 0.150 0.98877107793604 0.98877107793600 -4.662937E-014 0.200 0.98006657784124 0.98006657784078 -4.604095E-013 0.250 0.96891242171065 0.96891242171041 -2.322587E-013 0.300 0.95533648912561 0.95533648912587 2.609024E-013 0.350 0.93937271284738 0.93937271284784 4.606315E-013 0.400 0.92106099400289 0.92106099400308 1.980638E-013 0.450 0.90044710235268 0.90044710235243 -2.473577E-013 0.500 0.87758256189037 0.87758256188991 -4.586331E-013 0.550 0.85252452205951 0.85252452205926 -2.461364E-013 0.600 0.82533561490968 0.82533561490988 1.961764E-013 0.650 0.79608379854906 0.79608379854951 4.536371E-013 0.700 0.76484218728449 0.76484218728474 2.553513E-013 0.750 0.73168886887382 0.73168886887359 -2.267075E-013 0.800 0.69670670934717 0.69670670934672 -4.467537E-013 0.850 0.65998314588498 0.65998314588494 -4.485301E-014 0.900 0.62160996827066 0.62160996827111 4.444223E-013 0.950 0.58168308946388 0.58168308946379 -8.992806E-014
C++
Based on the C99 implementation above. The main improvement is that, because C++ containers handle memory for us, we can use a more functional style.
The two overloads of cheb_coef show a useful idiom for working with C++ templates; the non-template code, which does all the mathematical work, can be placed in a source file so that it is compiled only once (reducing code bloat from repeating substantial blocks of code). The template function is a minimal wrapper to call the non-template implementation.
The wrapper class ChebyshevApprox_ supports very terse user code.
#include <iostream>
#include <iomanip>
#include <string>
#include <cmath>
#include <utility>
#include <vector>
using namespace std;
static const double PI = acos(-1.0);
double affine_remap(const pair<double, double>& from, double x, const pair<double, double>& to)
{
return to.first + (x - from.first) * (to.second - to.first) / (from.second - from.first);
}
vector<double> cheb_coef(const vector<double>& f_vals)
{
const int n = f_vals.size();
const double theta = PI / n;
vector<double> retval(n, 0.0);
for (int ii = 0; ii < n; ++ii)
{
double f = f_vals[ii] * 2.0 / n;
const double phi = (ii + 0.5) * theta;
double c1 = cos(phi), s1 = sin(phi);
double c = 1.0, s = 0.0;
for (int j = 0; j < n; j++)
{
retval[j] += f * c;
// update c -> cos(j*phi) for next value of j
const double cNext = c * c1 - s * s1;
s = c * s1 + s * c1;
c = cNext;
}
}
return retval;
}
template<class F_> vector<double> cheb_coef(const F_& func, int n, const pair<double, double>& domain)
{
auto remap = [&](double x){return affine_remap({ -1.0, 1.0 }, x, domain); };
const double theta = PI / n;
vector<double> fVals(n);
for (int ii = 0; ii < n; ++ii)
fVals[ii] = func(remap(cos((ii + 0.5) * theta)));
return cheb_coef(fVals);
}
double cheb_eval(const vector<double>& coef, double x)
{
double a = 1.0, b = x, c;
double retval = 0.5 * coef[0] + b * coef[1];
for (auto pc = coef.begin() + 2; pc != coef.end(); a = b, b = c, ++pc)
{
c = 2.0 * b * x - a;
retval += (*pc) * c;
}
return retval;
}
double cheb_eval(const vector<double>& coef, const pair<double, double>& domain, double x)
{
return cheb_eval(coef, affine_remap(domain, x, { -1.0, 1.0 }));
}
struct ChebyshevApprox_
{
vector<double> coeffs_;
pair<double, double> domain_;
double operator()(double x) const { return cheb_eval(coeffs_, domain_, x); }
template<class F_> ChebyshevApprox_
(const F_& func,
int n,
const pair<double, double>& domain)
:
coeffs_(cheb_coef(func, n, domain)),
domain_(domain)
{ }
};
int main(void)
{
static const int N = 10;
ChebyshevApprox_ fApprox(cos, N, { 0.0, 1.0 });
cout << "Coefficients: " << setprecision(14);
for (const auto& c : fApprox.coeffs_)
cout << "\t" << c << "\n";
for (;;)
{
cout << "Enter x, or non-numeric value to quit:\n";
double x;
if (!(cin >> x))
return 0;
cout << "True value: \t" << cos(x) << "\n";
cout << "Approximate: \t" << fApprox(x) << "\n";
}
}
D
This imperative code retains some of the style of the original C version.
import std.math: PI, cos;
/// Map x from range [min, max] to [min_to, max_to].
real map(in real x, in real min_x, in real max_x, in real min_to, in real max_to)
pure nothrow @safe @nogc {
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to;
}
void chebyshevCoef(size_t N)(in real function(in real) pure nothrow @safe @nogc func,
in real min, in real max, ref real[N] coef)
pure nothrow @safe @nogc {
coef[] = 0.0;
foreach (immutable i; 0 .. N) {
immutable f = func(map(cos(PI * (i + 0.5f) / N), -1, 1, min, max)) * 2 / N;
foreach (immutable j, ref cj; coef)
cj += f * cos(PI * j * (i + 0.5f) / N);
}
}
/// f(x) = sum_{k=0}^{n - 1} c_k T_k(x) - c_0/2
real chebyshevApprox(size_t N)(in real x, in real min, in real max, in ref real[N] coef)
pure nothrow @safe @nogc if (N >= 2) {
real a = 1.0L,
b = map(x, min, max, -1, 1),
result = coef[0] / 2 + coef[1] * b;
immutable x2 = 2 * b;
foreach (immutable ci; coef[2 .. $]) {
// T_{n+1} = 2x T_n - T_{n-1}
immutable c = x2 * b - a;
result += ci * c;
a = b;
b = c;
}
return result;
}
void main() @safe {
import std.stdio: writeln, writefln;
enum uint N = 10;
real[N] c;
real min = 0, max = 1;
static real test(in real x) pure nothrow @safe @nogc { return x.cos; }
chebyshevCoef(&test, min, max, c);
writefln("Coefficients:\n%( %+2.25g\n%)", c);
enum nX = 20;
writeln("\nApproximation:\n x func(x) approx diff");
foreach (immutable i; 0 .. nX) {
immutable x = map(i, 0, nX, min, max);
immutable f = test(x);
immutable approx = chebyshevApprox(x, min, max, c);
writefln("%1.3f % 10.10f % 10.10f % 4.2e", x, f, approx, approx - f);
}
}
- Output:
Coefficients: +1.6471694753903136868 -0.23229937161517194216 -0.053715114622047555044 +0.0024582352669814797779 +0.00028211905743400579387 -7.7222291558103533853e-06 -5.898556452178771968e-07 +1.1521427332860788728e-08 +6.5963000382704222411e-10 -1.0022591914390921452e-11 Approximation: x func(x) approx diff 0.000 1.00000000000000000000 1.00000000000046961190 4.70e-13 0.050 0.99875026039496624654 0.99875026039487216781 -9.41e-14 0.100 0.99500416527802576609 0.99500416527848803832 4.62e-13 0.150 0.98877107793604228670 0.98877107793599569749 -4.66e-14 0.200 0.98006657784124163110 0.98006657784078136889 -4.60e-13 0.250 0.96891242171064478408 0.96891242171041249593 -2.32e-13 0.300 0.95533648912560601967 0.95533648912586667367 2.61e-13 0.350 0.93937271284737892005 0.93937271284783928305 4.60e-13 0.400 0.92106099400288508277 0.92106099400308274515 1.98e-13 0.450 0.90044710235267692169 0.90044710235242891114 -2.48e-13 0.500 0.87758256189037271615 0.87758256188991362600 -4.59e-13 0.550 0.85252452205950574283 0.85252452205925896211 -2.47e-13 0.600 0.82533561490967829723 0.82533561490987400509 1.96e-13 0.650 0.79608379854905582896 0.79608379854950937939 4.54e-13 0.700 0.76484218728448842626 0.76484218728474395029 2.56e-13 0.750 0.73168886887382088633 0.73168886887359430061 -2.27e-13 0.800 0.69670670934716542091 0.69670670934671868322 -4.47e-13 0.850 0.65998314588498217039 0.65998314588493717370 -4.50e-14 0.900 0.62160996827066445648 0.62160996827110870299 4.44e-13 0.950 0.58168308946388349416 0.58168308946379353278 -9.00e-14
The same code, with N = 16:
Coefficients: +1.6471694753903136868 -0.23229937161517194214 -0.053715114622047555035 +0.0024582352669814797982 +0.00028211905743400571932 -7.722229155810705751e-06 -5.898556452177348953e-07 +1.1521427330794028337e-08 +6.5963022091481034181e-10 -1.0016894235462866363e-11 -4.5865582517937500406e-13 +5.6974586994888026802e-15 +2.1752822525027137867e-16 -2.3140940118987485263e-18 -1.0333801956502464137e-19 +2.5410988417629010172e-20 Approximation: x func(x) approx diff 0.000 1.00000000000000000000 1.00000000000000000030 3.25e-19 0.050 0.99875026039496624654 0.99875026039496624646 -1.08e-19 0.100 0.99500416527802576609 0.99500416527802576557 -5.42e-19 0.150 0.98877107793604228670 0.98877107793604228636 -3.79e-19 0.200 0.98006657784124163110 0.98006657784124163127 1.08e-19 0.250 0.96891242171064478408 0.96891242171064478451 3.79e-19 0.300 0.95533648912560601967 0.95533648912560601967 0.00e+00 0.350 0.93937271284737892005 0.93937271284737891962 -3.79e-19 0.400 0.92106099400288508277 0.92106099400288508260 -2.17e-19 0.450 0.90044710235267692169 0.90044710235267692169 5.42e-20 0.500 0.87758256189037271615 0.87758256189037271632 2.17e-19 0.550 0.85252452205950574283 0.85252452205950574274 -5.42e-20 0.600 0.82533561490967829723 0.82533561490967829697 -2.17e-19 0.650 0.79608379854905582896 0.79608379854905582861 -3.25e-19 0.700 0.76484218728448842626 0.76484218728448842630 5.42e-20 0.750 0.73168886887382088633 0.73168886887382088637 5.42e-20 0.800 0.69670670934716542091 0.69670670934716542087 -5.42e-20 0.850 0.65998314588498217039 0.65998314588498217022 -1.63e-19 0.900 0.62160996827066445648 0.62160996827066445674 2.71e-19 0.950 0.58168308946388349416 0.58168308946388349403 -1.63e-19
EasyLang
numfmt 0 5
a = 0
b = 1
n = 10
len coef[] n
len cheby[] n
for i range n
coef[i] = cos (180 / pi * (cos (180 / n * (i + 1 / 2)) * (b - a) / 2 + (b + a) / 2))
.
for i range n
w = 0
for j range n
w += coef[j] * cos (180 / n * i * (j + 1 / 2))
.
cheby[i] = w * 2 / n
print cheby[i]
.
Go
Wikipedia gives a formula for coefficients in a section "Example 1". Read past the bit about the inner product to where it gives the technique based on the discrete orthogonality condition. The N of the WP formulas is the parameter nNodes in the code here. It is not necessarily the same as n, the number of polynomial coefficients, the parameter nCoeff here.
The evaluation method is the Clenshaw algorithm.
Two variances here from the WP presentation and most mathematical presentations follow other examples on this page and so keep output directly comparable. One variance is that the Kronecker delta factor is dropped, which has the effect of doubling the first coefficient. This simplifies both coefficient generation and polynomial evaluation. A further variance is that there is no scaling for the range of function values. The result is that coefficients are not necessarily bounded by 1 (2 for the first coefficient) but by the maximum function value over the argument range from min to max (or twice that for the first coefficient.)
package main
import (
"fmt"
"math"
)
type cheb struct {
c []float64
min, max float64
}
func main() {
fn := math.Cos
c := newCheb(0, 1, 10, 10, fn)
fmt.Println("coefficients:")
for _, c := range c.c {
fmt.Printf("% .15f\n", c)
}
fmt.Println("\nx computed approximated computed-approx")
const n = 10
for i := 0.; i <= n; i++ {
x := (c.min*(n-i) + c.max*i) / n
computed := fn(x)
approx := c.eval(x)
fmt.Printf("%.1f %12.8f %12.8f % .3e\n",
x, computed, approx, computed-approx)
}
}
func newCheb(min, max float64, nCoeff, nNodes int, fn func(float64) float64) *cheb {
c := &cheb{
c: make([]float64, nCoeff),
min: min,
max: max,
}
f := make([]float64, nNodes)
p := make([]float64, nNodes)
z := .5 * (max + min)
r := .5 * (max - min)
for k := 0; k < nNodes; k++ {
p[k] = math.Pi * (float64(k) + .5) / float64(nNodes)
f[k] = fn(z + math.Cos(p[k])*r)
}
n2 := 2 / float64(nNodes)
for j := 0; j < nCoeff; j++ {
sum := 0.
for k := 0; k < nNodes; k++ {
sum += f[k] * math.Cos(float64(j)*p[k])
}
c.c[j] = sum * n2
}
return c
}
func (c *cheb) eval(x float64) float64 {
x1 := (2*x - c.min - c.max) / (c.max - c.min)
x2 := 2 * x1
var s, t float64
for j := len(c.c) - 1; j >= 1; j-- {
t, s = x2*t-s+c.c[j], t
}
return x1*t - s + .5*c.c[0]
}
- Output:
coefficients: 1.647169475390314 -0.232299371615172 -0.053715114622048 0.002458235266982 0.000282119057434 -0.000007722229156 -0.000000589855645 0.000000011521427 0.000000000659630 -0.000000000010022 x computed approximated computed-approx 0.0 1.00000000 1.00000000 -4.685e-13 0.1 0.99500417 0.99500417 -4.620e-13 0.2 0.98006658 0.98006658 4.601e-13 0.3 0.95533649 0.95533649 -2.607e-13 0.4 0.92106099 0.92106099 -1.972e-13 0.5 0.87758256 0.87758256 4.587e-13 0.6 0.82533561 0.82533561 -1.965e-13 0.7 0.76484219 0.76484219 -2.552e-13 0.8 0.69670671 0.69670671 4.470e-13 0.9 0.62160997 0.62160997 -4.449e-13 1.0 0.54030231 0.54030231 -4.476e-13
Groovy
class ChebyshevCoefficients {
static 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
}
static void chebyshevCoef(Closure<Double> func, double min, double max, double[] coef) {
final int N = coef.length
for (int i = 0; i < N; i++) {
double m = map(Math.cos(Math.PI * (i + 0.5f) / N), -1, 1, min, max)
double f = func(m) * 2 / N
for (int j = 0; j < N; j++) {
coef[j] += f * Math.cos(Math.PI * j * (i + 0.5f) / N)
}
}
}
static void main(String[] args) {
final int N = 10
double[] c = new double[N]
double min = 0, max = 1
chebyshevCoef(Math.&cos, min, max, c)
println("Coefficients:")
for (double d : c) {
println(d)
}
}
}
- Output:
Coefficients: 1.6471694753903139 -0.23229937161517178 -0.0537151146220477 0.002458235266981773 2.8211905743405485E-4 -7.722229156320592E-6 -5.898556456745974E-7 1.1521427770166959E-8 6.59630183807991E-10 -1.0021913854352249E-11
J
From 'J for C Programmers: Calculating Chebyshev Coefficients [[1]]
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
)
Calculate coefficients:
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
Java
Partial translation of C via D
import static java.lang.Math.*;
import java.util.function.Function;
public class ChebyshevCoefficients {
static 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;
}
static void chebyshevCoef(Function<Double, Double> func, double min,
double max, double[] coef) {
int N = coef.length;
for (int i = 0; i < N; i++) {
double m = map(cos(PI * (i + 0.5f) / N), -1, 1, min, max);
double f = func.apply(m) * 2 / N;
for (int j = 0; j < N; j++) {
coef[j] += f * cos(PI * j * (i + 0.5f) / N);
}
}
}
public static void main(String[] args) {
final int N = 10;
double[] c = new double[N];
double min = 0, max = 1;
chebyshevCoef(x -> cos(x), min, max, c);
System.out.println("Coefficients:");
for (double d : c)
System.out.println(d);
}
}
Coefficients: 1.6471694753903139 -0.23229937161517178 -0.0537151146220477 0.002458235266981773 2.8211905743405485E-4 -7.722229156320592E-6 -5.898556456745974E-7 1.1521427770166959E-8 6.59630183807991E-10 -1.0021913854352249E-11
jq
Adapted from Wren
Works with gojq, the Go implementation of jq
Preliminaries
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
def rpad($len; $fill): tostring | ($len - length) as $l | . + ($fill * $l)[:$l];
# Format a decimal number so that there are at least `left` characters
# to the left of the decimal point, and at most `right` characters to its right.
# No left-truncation occurs, so `left` can be specified as 0 to prevent left-padding.
# If tostring has an "e" then eparse as defined below is used.
def pp(left; right):
def lpad: if (left > length) then ((left - length) * " ") + . else . end;
def eparse: index("e") as $ix | (.[:$ix]|pp(left;right)) + .[$ix:];
tostring as $s
| $s
| if test("e") then eparse
else index(".") as $ix
| ((if $ix then $s[0:$ix] else $s end) | lpad) + "." +
(if $ix then $s[$ix+1:] | .[0:right] else "" end)
end;
Chebyshev Coefficients
def mapRange($x; $min; $max; $minTo; $maxTo):
(($x - $min)/($max - $min))*($maxTo - $minTo) + $minTo;
def chebCoeffs(func; n; min; max):
(1 | atan * 4) as $pi
| reduce range(0;n) as $i ([]; # coeffs
((mapRange( ($pi * ($i + 0.5) / n)|cos; -1; 1; min; max) | func) * 2 / n) as $f
| reduce range(0;n) as $j (.;
.[$j] += $f * ($pi * $j * (($i + 0.5) / n)|cos)) );
def chebApprox(x; n; min; max; coeffs):
if n < 2 or (coeffs|length) < 2 then "'n' can't be less than 2." | error
else { a: 1,
b: mapRange(x; min; max; -1; 1) }
| .res = coeffs[0]/2 + coeffs[1]*.b
| .xx = 2 * .b
| reduce range(2;n) as $i (.;
(.xx * .b - .a) as $c
| .res += coeffs[$i]*$c)
| .a = .b
| .b = $c)
| .res
end ;
def task:
[10, 0, 1] as [$n, $min, $max]
| chebCoeffs(cos; $n; $min; $max) as $coeffs
| "Coefficients:",
($coeffs[]|pp(2;14)),
"\nApproximations:\n x func(x) approx diff",
(range(0;21) as $i
| mapRange($i; 0; 20; $min; $max) as $x
| ($x|cos) as $f
| chebApprox($x; $n; $min; $max; $coeffs) as $approx
| ($approx - $f) as $diff
| [ ($x|pp(0;3)|rpad( 4;"0")),
($f|pp(0;8)|rpad(10;"0")),
($approx|pp(0;8)),
($diff |pp(2;2)) ]
| join(" ") );
task
- Output:
Coefficients: 1.64716947539031 -0.23229937161517 -0.05371511462204 0.00245823526698 0.00028211905743 -7.72222915562670e-06 -5.89855645688475e-07 1.15214280338449e-08 6.59629580124221e-10 -1.00220526322303e-11 Approximations: x func(x) approx diff 0.00 1.00000000 1.00000000 4.66e-13 0.05 0.99875026 0.99875026 -9.21e-14 0.10 0.99500416 0.99500416 4.62e-13 0.15 0.98877107 0.98877107 -4.74e-14 0.20 0.98006657 0.98006657 -4.60e-13 0.25 0.96891242 0.96891242 -2.32e-13 0.30 0.95533648 0.95533648 2.61e-13 0.35 0.93937271 0.93937271 4.60e-13 0.40 0.92106099 0.92106099 1.98e-13 0.45 0.90044710 0.90044710 -2.47e-13 0.50 0.87758256 0.87758256 -4.59e-13 0.55 0.85252452 0.85252452 -2.46e-13 0.60 0.82533561 0.82533561 1.95e-13 0.65 0.79608379 0.79608379 4.53e-13 0.70 0.76484218 0.76484218 2.55e-13 0.75 0.73168886 0.73168886 -2.26e-13 0.80 0.69670670 0.69670670 -4.46e-13 0.85 0.65998314 0.65998314 -4.45e-14 0.90 0.62160996 0.62160996 4.44e-13 0.95 0.58168308 0.58168308 -9.01e-14 1.00 0.54030230 0.54030230 4.47e-13
Julia
mutable struct Cheb
c::Vector{Float64}
min::Float64
max::Float64
end
function Cheb(min::Float64, max::Float64, ncoeff::Int, nnodes::Int, fn::Function)::Cheb
c = Cheb(Vector{Float64}(ncoeff), min, max)
f = Vector{Float64}(nnodes)
p = Vector{Float64}(nnodes)
z = (max + min) / 2
r = (max - min) / 2
for k in 0:nnodes-1
p[k+1] = π * (k + 0.5) / nnodes
f[k+1] = fn(z + cos(p[k+1]) * r)
end
n2 = 2 / nnodes
for j in 0:nnodes-1
s = sum(fk * cos(j * pk) for (fk, pk) in zip(f, p))
c.c[j+1] = s * n2
end
return c
end
function evaluate(c::Cheb, x::Float64)::Float64
x1 = (2x - c.max - c.min) / (c.max - c.min)
x2 = 2x1
t = s = 0
for j in length(c.c):-1:2
t, s = x2 * t - s + c.c[j], t
end
return x1 * t - s + c.c[1] / 2
end
fn = cos
c = Cheb(0.0, 1.0, 10, 10, fn)
# coefs
println("Coefficients:")
for x in c.c
@printf("% .15f\n", x)
end
# values
println("\nx computed approximated computed-approx")
const n = 10
for i in 0.0:n
x = (c.min * (n - i) + c.max * i) / n
computed = fn(x)
approx = evaluate(c, x)
@printf("%.1f %12.8f %12.8f % .3e\n", x, computed, approx, computed - approx)
end
- Output:
Coefficients: 1.647169475390314 -0.232299371615172 -0.053715114622048 0.002458235266981 0.000282119057434 -0.000007722229156 -0.000000589855645 0.000000011521427 0.000000000659630 -0.000000000010022 x computed approximated computed-approx 0.0 1.00000000 1.00000000 -4.685e-13 0.1 0.99500417 0.99500417 -4.620e-13 0.2 0.98006658 0.98006658 4.601e-13 0.3 0.95533649 0.95533649 -2.605e-13 0.4 0.92106099 0.92106099 -1.970e-13 0.5 0.87758256 0.87758256 4.586e-13 0.6 0.82533561 0.82533561 -1.967e-13 0.7 0.76484219 0.76484219 -2.551e-13 0.8 0.69670671 0.69670671 4.470e-13 0.9 0.62160997 0.62160997 -4.449e-13 1.0 0.54030231 0.54030231 -4.476e-13
Kotlin
// version 1.1.2
typealias DFunc = (Double) -> Double
fun mapRange(x: Double, min: Double, max: Double, minTo: Double, maxTo:Double): Double {
return (x - min) / (max - min) * (maxTo - minTo) + minTo
}
fun chebCoeffs(func: DFunc, n: Int, min: Double, max: Double): DoubleArray {
val coeffs = DoubleArray(n)
for (i in 0 until n) {
val f = func(mapRange(Math.cos(Math.PI * (i + 0.5) / n), -1.0, 1.0, min, max)) * 2.0 / n
for (j in 0 until n) coeffs[j] += f * Math.cos(Math.PI * j * (i + 0.5) / n)
}
return coeffs
}
fun chebApprox(x: Double, n: Int, min: Double, max: Double, coeffs: DoubleArray): Double {
require(n >= 2 && coeffs.size >= 2)
var a = 1.0
var b = mapRange(x, min, max, -1.0, 1.0)
var res = coeffs[0] / 2.0 + coeffs[1] * b
val xx = 2 * b
var i = 2
while (i < n) {
val c = xx * b - a
res += coeffs[i] * c
a = b
b = c
i++
}
return res
}
fun main(args: Array<String>) {
val n = 10
val min = 0.0
val max = 1.0
val coeffs = chebCoeffs(Math::cos, n, min, max)
println("Coefficients:")
for (coeff in coeffs) println("%+1.15g".format(coeff))
println("\nApproximations:\n x func(x) approx diff")
for (i in 0..20) {
val x = mapRange(i.toDouble(), 0.0, 20.0, min, max)
val f = Math.cos(x)
val approx = chebApprox(x, n, min, max, coeffs)
System.out.printf("%1.3f %1.8f %1.8f % 4.1e\n", x, f, approx, approx - f)
}
}
- Output:
Coefficients: +1.64716947539031 -0.232299371615172 -0.0537151146220477 +0.00245823526698177 +0.000282119057434055 -7.72222915632059e-06 -5.89855645674597e-07 +1.15214277701670e-08 +6.59630183807991e-10 -1.00219138543522e-11 Approximations: x func(x) approx diff 0.000 1.00000000 1.00000000 4.7e-13 0.050 0.99875026 0.99875026 -9.4e-14 0.100 0.99500417 0.99500417 4.6e-13 0.150 0.98877108 0.98877108 -4.7e-14 0.200 0.98006658 0.98006658 -4.6e-13 0.250 0.96891242 0.96891242 -2.3e-13 0.300 0.95533649 0.95533649 2.6e-13 0.350 0.93937271 0.93937271 4.6e-13 0.400 0.92106099 0.92106099 2.0e-13 0.450 0.90044710 0.90044710 -2.5e-13 0.500 0.87758256 0.87758256 -4.6e-13 0.550 0.85252452 0.85252452 -2.5e-13 0.600 0.82533561 0.82533561 2.0e-13 0.650 0.79608380 0.79608380 4.5e-13 0.700 0.76484219 0.76484219 2.5e-13 0.750 0.73168887 0.73168887 -2.3e-13 0.800 0.69670671 0.69670671 -4.5e-13 0.850 0.65998315 0.65998315 -4.4e-14 0.900 0.62160997 0.62160997 4.5e-13 0.950 0.58168309 0.58168309 -9.0e-14 1.000 0.54030231 0.54030231 4.5e-13
Lua
function map(x, min_x, max_x, min_to, max_to)
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to
end
function chebyshevCoef(func, minn, maxx, coef)
local N = table.getn(coef)
for j=1,N do
local i = j - 1
local m = map(math.cos(math.pi * (i + 0.5) / N), -1, 1, minn, maxx)
local f = func(m) * 2 / N
for k=1,N do
local p = k -1
coef[k] = coef[k] + f * math.cos(math.pi * p * (i + 0.5) / N)
end
end
end
function main()
local N = 10
local c = {}
local minn = 0.0
local maxx = 1.0
for i=1,N do
table.insert(c, 0)
end
chebyshevCoef(function (x) return math.cos(x) end, minn, maxx, c)
print("Coefficients:")
for i,d in pairs(c) do
print(d)
end
end
main()
- Output:
Coefficients: 1.6471694753903 -0.23229937161517 -0.053715114622048 0.0024582352669818 0.00028211905743405 -7.7222291563483e-006 -5.898556456746e-007 1.1521427756289e-008 6.5963018380799e-010 -1.0021913854352e-011
Microsoft Small Basic
' N Chebyshev coefficients for the range 0 to 1 - 18/07/2018
pi=Math.pi
a=0
b=1
n=10
For i=0 To n-1
coef[i]=Math.cos(Math.cos(pi/n*(i+1/2))*(b-a)/2+(b+a)/2)
EndFor
For i=0 To n-1
w=0
For j=0 To n-1
w=w+coef[j]*Math.cos(pi/n*i*(j+1/2))
EndFor
cheby[i]=w*2/n
t=" "
If cheby[i]<=0 Then
t=""
EndIf
TextWindow.WriteLine(i+" : "+t+cheby[i])
EndFor
- Output:
0 : 1,6471694753903136 1 : -0,2322993716151700684187787635 2 : -0,0537151146220494010749946688 3 : 0,0024582352669837594966069584 4 : 0,0002821190574317389206759282 5 : -0,0000077222291539069653168878 6 : -0,0000005898556481086082412444 7 : 0,0000000115214300591398939205 8 : 0,0000000006596278553822696656 9 : -0,0000000000100189955816952521
Nim
import lenientops, math, strformat, sugar
type Cheb = object
c: seq[float]
min, max: float
func initCheb(min, max: float; nCoeff, nNodes: int; fn: float -> float): Cheb =
result = Cheb(c: newSeq[float](nCoeff), min: min, max: max)
var f, p = newSeq[float](nNodes)
let z = 0.5 * (max + min)
let r = 0.5 * (max - min)
for k in 0..<nNodes:
p[k] = PI * (k + 0.5) / nNodes
f[k] = fn(z + cos(p[k]) * r)
let n2 = 2 / nNodes
for j in 0..<nCoeff:
var sum = 0.0
for k in 0..<nNodes:
sum += f[k] * cos(j * p[k])
result.c[j] = sum * n2
func eval(cheb: Cheb; x: float): float =
let x1 = (2 * x - cheb.min - cheb.max) / (cheb.max - cheb.min)
let x2 = 2 * x1
var s, t: float
for j in countdown(cheb.c.high, 1):
s = x2 * t - s + cheb.c[j]
swap s, t
result = x1 * t - s + 0.5 * cheb.c[0]
when isMainModule:
let fn: float -> float = cos
let cheb = initCheb(0, 1, 10, 10, fn)
echo "Coefficients:"
for c in cheb.c:
echo &"{c: .15f}"
echo "\n x computed approximated computed-approx"
const N = 10
for i in 0..N:
let x = (cheb.min * (N - i) + cheb.max * i) / N
let computed = fn(x)
let approx = cheb.eval(x)
echo &"{x:.1f} {computed:12.8f} {approx:12.8f} {computed-approx: .3e}"
- Output:
Coefficients: 1.647169475390314 -0.232299371615172 -0.053715114622048 0.002458235266981 0.000282119057434 -0.000007722229156 -0.000000589855645 0.000000011521427 0.000000000659630 -0.000000000010022 x computed approximated computed-approx 0.0 1.00000000 1.00000000 -4.685e-13 0.1 0.99500417 0.99500417 -4.620e-13 0.2 0.98006658 0.98006658 4.601e-13 0.3 0.95533649 0.95533649 -2.605e-13 0.4 0.92106099 0.92106099 -1.970e-13 0.5 0.87758256 0.87758256 4.586e-13 0.6 0.82533561 0.82533561 -1.967e-13 0.7 0.76484219 0.76484219 -2.551e-13 0.8 0.69670671 0.69670671 4.470e-13 0.9 0.62160997 0.62160997 -4.450e-13 1.0 0.54030231 0.54030231 -4.476e-13
Perl
use constant PI => 2 * atan2(1, 0);
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
}
printf "%+13.7e\n", $_ for chebft(sub{cos($_[0])}, 0, 1, 10);
- 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
Phix
function Cheb(atom cmin, cmax, integer ncoeff, nnodes) sequence c = repeat(0,ncoeff), f = repeat(0,nnodes), p = repeat(0,nnodes) atom z = (cmax + cmin) / 2, r = (cmax - cmin) / 2 for k=1 to nnodes do p[k] = PI * ((k-1) + 0.5) / nnodes f[k] = cos(z + cos(p[k]) * r) end for atom n2 = 2 / nnodes for j=1 to nnodes do atom s := 0 for k=1 to nnodes do s += f[k] * cos((j-1)*p[k]) end for c[j] = s * n2 end for return c end function function evaluate(sequence c, atom cmin, cmax, x) atom x1 = (2*x - cmax - cmin) / (cmax - cmin), x2 = 2*x1, t = 0, s = 0 for j=length(c) to 2 by -1 do {t, s} = {x2 * t - s + c[j], t} end for return x1 * t - s + c[1] / 2 end function atom cmin = 0.0, cmax = 1.0 sequence c = Cheb(cmin, cmax, 10, 10) printf(1, "Coefficients:\n") pp(c,{pp_Nest,1,pp_FltFmt,"%18.15f"}) printf(1,"\nx computed approximated computed-approx\n") constant n = 10 for i=0 to 10 do atom x = (cmin * (n - i) + cmax * i) / n, calc = cos(x), est = evaluate(c, cmin, cmax, x) printf(1,"%.1f %12.8f %12.8f %10.3e\n", {x, calc, est, calc-est}) end for
- Output:
Coefficients: { 1.647169475390314, -0.232299371615172, -0.053715114622048, 0.002458235266981, 0.000282119057434, -0.000007722229156, -0.000000589855645, 0.000000011521427, 0.000000000659630, -0.000000000010022} x computed approximated computed-approx 0.0 1.00000000 1.00000000 -4.686e-13 0.1 0.99500417 0.99500417 -4.620e-13 0.2 0.98006658 0.98006658 4.600e-13 0.3 0.95533649 0.95533649 -2.604e-13 0.4 0.92106099 0.92106099 -1.970e-13 0.5 0.87758256 0.87758256 4.587e-13 0.6 0.82533561 0.82533561 -1.968e-13 0.7 0.76484219 0.76484219 -2.551e-13 0.8 0.69670671 0.69670671 4.470e-13 0.9 0.62160997 0.62160997 -4.450e-13 1.0 0.54030231 0.54030231 -4.477e-13
Python
import math
def test_func(x):
return math.cos(x)
def mapper(x, min_x, max_x, min_to, max_to):
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to
def cheb_coef(func, n, min, max):
coef = [0.0] * n
for i in xrange(n):
f = func(mapper(math.cos(math.pi * (i + 0.5) / n), -1, 1, min, max)) * 2 / n
for j in xrange(n):
coef[j] += f * math.cos(math.pi * j * (i + 0.5) / n)
return coef
def cheb_approx(x, n, min, max, coef):
a = 1
b = mapper(x, min, max, -1, 1)
c = float('nan')
res = coef[0] / 2 + coef[1] * b
x = 2 * b
i = 2
while i < n:
c = x * b - a
res = res + coef[i] * c
(a, b) = (b, c)
i += 1
return res
def main():
N = 10
min = 0
max = 1
c = cheb_coef(test_func, N, min, max)
print "Coefficients:"
for i in xrange(N):
print " % lg" % c[i]
print "\n\nApproximation:\n x func(x) approx diff"
for i in xrange(20):
x = mapper(i, 0.0, 20.0, min, max)
f = test_func(x)
approx = cheb_approx(x, N, min, max, c)
print "%1.3f %10.10f %10.10f % 4.2e" % (x, f, approx, approx - f)
return None
main()
- Output:
Coefficients: 1.64717 -0.232299 -0.0537151 0.00245824 0.000282119 -7.72223e-06 -5.89856e-07 1.15214e-08 6.5963e-10 -1.00219e-11 Approximation: x func(x) approx diff 0.000 1.0000000000 1.0000000000 4.68e-13 0.050 0.9987502604 0.9987502604 -9.36e-14 0.100 0.9950041653 0.9950041653 4.62e-13 0.150 0.9887710779 0.9887710779 -4.73e-14 0.200 0.9800665778 0.9800665778 -4.60e-13 0.250 0.9689124217 0.9689124217 -2.32e-13 0.300 0.9553364891 0.9553364891 2.62e-13 0.350 0.9393727128 0.9393727128 4.61e-13 0.400 0.9210609940 0.9210609940 1.98e-13 0.450 0.9004471024 0.9004471024 -2.47e-13 0.500 0.8775825619 0.8775825619 -4.58e-13 0.550 0.8525245221 0.8525245221 -2.46e-13 0.600 0.8253356149 0.8253356149 1.96e-13 0.650 0.7960837985 0.7960837985 4.53e-13 0.700 0.7648421873 0.7648421873 2.54e-13 0.750 0.7316888689 0.7316888689 -2.28e-13 0.800 0.6967067093 0.6967067093 -4.47e-13 0.850 0.6599831459 0.6599831459 -4.37e-14 0.900 0.6216099683 0.6216099683 4.46e-13 0.950 0.5816830895 0.5816830895 -8.99e-14
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
- 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)
Raku
(formerly Perl 6)
sub chebft ( Code $func, Real \a, Real \b, Int \n ) {
my \bma = ½ × (b - a);
my \bpa = ½ × (b + a);
my @pi-n = ( ^n »+» ½ ) »×» (π/n);
my @f = ( @pi_n».cos »×» bma »+» bpa )».&$func;
my @sums = (^n).map: { [+] @f »×« ( @pi-n »×» $_ )».cos };
@sums »×» (2/n)
}
say chebft(&cos, 0, 1, 10)».fmt: '%+13.7e';
- 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
REXX
This REXX program is a translation of the C program plus added optimizations.
Pafnuty Lvovich Chebysheff: Chebyshev [English transliteration] Chebysheff [ " " ] Chebyshov [ " " ] Tchebychev [French " ] Tchebysheff [ " " ] Tschebyschow [German " ] Tschebyschev [ " " ] Tschebyschef [ " " ] Tschebyscheff [ " " ]
The numeric precision is dependent on the number of decimal digits specified in the value of pi.
/*REXX program calculates N Chebyshev coefficients for the range 0 ──► 1 (inclusive)*/
numeric digits length( pi() ) - length(.) /*DIGITS default is nine, but use 71. */
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 " " " " " */
fac= 2 / N; pin= pi / N /*calculate a couple handy─dandy values*/
Dma= (b-a) / 2 /*calculate one─half of the difference.*/
Dpa= (b+a) / 2 /* " " " " sum. */
do k=0 for N; f.k= cos( cos( pin * (k + .5) ) * Dma + Dpa)
end /*k*/
do j=0 for N; z= pin * j /*The LEFT('', ···) ────────►──────┐ */
$= 0 /* clause is used to align │ */
do m=0 for N /* the non─negative values with ↓ */
$= $ + f.m * cos(z*(m +.5)) /* the negative values. │ */
end /*m*/ /* ┌─────◄──────┘ */
cheby.j= fac * $ /* ↓ */
say right(j, length(N) +3) " Chebyshev coefficient is:" left('', cheby.j >= 0),
format(cheby.j, , 30) /*only show 30 decimal digits of coeff.*/
end /*j*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure; parse arg x; numeric digits digits()+9; x=r2r(x); a=abs(x); numeric fuzz 5
if a=pi then return -1; if a=pi*.5 | a=pi*2 then return 0; pit= pi/3; z= 1
if a=pit then return .5; if a=pit*2 then return -.5; q= x*x; _= 1
do k=2 by 2 until p=z; p=z; _= -_ * q/(k*k - k); z= z+_; end; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164;return pi
r2r: return arg(1) // (pi() * 2) /*normalize radians ───► a unit circle.*/
- 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
- output when using the following input of: , , 20
0 Chebyshev coefficient is: 1.647169475390313686961473816799 1 Chebyshev coefficient is: -0.232299371615171942121038341150 2 Chebyshev coefficient is: -0.053715114622047555071596207909 3 Chebyshev coefficient is: 0.002458235266981479866768726383 4 Chebyshev coefficient is: 0.000282119057434005702429677244 5 Chebyshev coefficient is: -0.000007722229155810577212604038 6 Chebyshev coefficient is: -5.898556452177850238987693546709E-7 7 Chebyshev coefficient is: 1.152142733081886533841160480101E-8 8 Chebyshev coefficient is: 6.596302208686010678189261798322E-10 9 Chebyshev coefficient is: -1.001689435637395512060196156843E-11 10 Chebyshev coefficient is: -4.586557765969596848147502951921E-13 11 Chebyshev coefficient is: 5.697353072301630964243748212466E-15 12 Chebyshev coefficient is: 2.173565878297512401879760404343E-16 13 Chebyshev coefficient is: -2.284293234863639106096540267786E-18 14 Chebyshev coefficient is: -7.468956910165861862760811388638E-20 15 Chebyshev coefficient is: 6.802288097339388765485830636223E-22 16 Chebyshev coefficient is: 1.945994872442404773393679283660E-23 17 Chebyshev coefficient is: -1.563704507245591241161562138364E-25 18 Chebyshev coefficient is: -3.976201538410589537318561880598E-27 19 Chebyshev coefficient is: 2.859065292763079576513213370136E-29
Ruby
def mapp(x, min_x, max_x, min_to, max_to)
return (x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to
end
def chebyshevCoef(func, min, max, coef)
n = coef.length
for i in 0 .. n-1 do
m = mapp(Math.cos(Math::PI * (i + 0.5) / n), -1, 1, min, max)
f = func.call(m) * 2 / n
for j in 0 .. n-1 do
coef[j] = coef[j] + f * Math.cos(Math::PI * j * (i + 0.5) / n)
end
end
end
N = 10
def main
c = Array.new(N, 0)
min = 0
max = 1
chebyshevCoef(lambda { |x| Math.cos(x) }, min, max, c)
puts "Coefficients:"
puts c
end
main()
- Output:
Coefficients: 1.6471694753903139 -0.23229937161517178 -0.0537151146220477 0.002458235266981773 0.00028211905743405485 -7.722229156348348e-06 -5.898556456745974e-07 1.1521427756289171e-08 6.59630183807991e-10 -1.0021913854352249e-11
Scala
- Output:
Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).
import scala.math.{Pi, cos}
object ChebyshevCoefficients extends App {
final val N = 10
final val (min, max) = (0, 1)
val c = new Array[Double](N)
def chebyshevCoef(func: Double => Double,
min: Double,
max: Double,
coef: Array[Double]): Unit =
for (i <- coef.indices) {
def map(x: Double,
min_x: Double,
max_x: Double,
min_to: Double,
max_to: Double): Double =
(x - min_x) / (max_x - min_x) * (max_to - min_to) + min_to
val m = map(cos(Pi * (i + 0.5f) / N), -1, 1, min, max)
def f = func.apply(m) * 2 / N
for (j <- coef.indices) coef(j) += f * cos(Pi * j * (i + 0.5f) / N)
}
chebyshevCoef((x: Double) => cos(x), min, max, c)
println("Coefficients:")
c.foreach(d => println(f"$d%23.16e"))
}
Sidef
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) »*» (Number.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);
}
- 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
Swift
import Foundation
typealias DFunc = (Double) -> Double
func mapRange(x: Double, min: Double, max: Double, minTo: Double, maxTo: Double) -> Double {
return (x - min) / (max - min) * (maxTo - minTo) + minTo
}
func chebCoeffs(fun: DFunc, n: Int, min: Double, max: Double) -> [Double] {
var res = [Double](repeating: 0, count: n)
for i in 0..<n {
let dI = Double(i)
let dN = Double(n)
let f = fun(mapRange(x: cos(.pi * (dI + 0.5) / dN), min: -1, max: 1, minTo: min, maxTo: max)) * 2.0 / dN
for j in 0..<n {
res[j] += f * cos(.pi * Double(j) * (dI + 0.5) / dN)
}
}
return res
}
func chebApprox(x: Double, n: Int, min: Double, max: Double, coeffs: [Double]) -> Double {
var a = 1.0
var b = mapRange(x: x, min: min, max: max, minTo: -1, maxTo: 1)
var res = coeffs[0] / 2.0 + coeffs[1] * b
let xx = 2 * b
var i = 2
while i < n {
let c = xx * b - a
res += coeffs[i] * c
(a, b) = (b, c)
i += 1
}
return res
}
let coeffs = chebCoeffs(fun: cos, n: 10, min: 0, max: 1)
print("Coefficients")
for coeff in coeffs {
print(String(format: "%+1.15g", coeff))
}
print("\nApproximations:\n x func(x) approx diff")
for i in stride(from: 0.0, through: 20, by: 1) {
let x = mapRange(x: i, min: 0, max: 20, minTo: 0, maxTo: 1)
let f = cos(x)
let approx = chebApprox(x: x, n: 10, min: 0, max: 1, coeffs: coeffs)
print(String(format: "%1.3f %1.8f %1.8f % 4.1e", x, f, approx, approx - f))
}
- Output:
Coefficients +1.64716947539031 -0.232299371615172 -0.0537151146220476 +0.00245823526698177 +0.000282119057434055 -7.72222915632059e-06 -5.89855645688475e-07 +1.15214277562892e-08 +6.59630204624673e-10 -1.0021858343201e-11 Approximations: x func(x) approx diff 0.000 1.00000000 1.00000000 4.7e-13 0.050 0.99875026 0.99875026 -9.3e-14 0.100 0.99500417 0.99500417 4.6e-13 0.150 0.98877108 0.98877108 -4.7e-14 0.200 0.98006658 0.98006658 -4.6e-13 0.250 0.96891242 0.96891242 -2.3e-13 0.300 0.95533649 0.95533649 2.6e-13 0.350 0.93937271 0.93937271 4.6e-13 0.400 0.92106099 0.92106099 2.0e-13 0.450 0.90044710 0.90044710 -2.5e-13 0.500 0.87758256 0.87758256 -4.6e-13 0.550 0.85252452 0.85252452 -2.5e-13 0.600 0.82533561 0.82533561 2.0e-13 0.650 0.79608380 0.79608380 4.5e-13 0.700 0.76484219 0.76484219 2.5e-13 0.750 0.73168887 0.73168887 -2.3e-13 0.800 0.69670671 0.69670671 -4.5e-13 0.850 0.65998315 0.65998315 -4.4e-14 0.900 0.62160997 0.62160997 4.5e-13 0.950 0.58168309 0.58168309 -9.0e-14 1.000 0.54030231 0.54030231 4.5e-13
VBScript
To run in console mode with cscript.
' N Chebyshev coefficients for the range 0 to 1
Dim coef(10),cheby(10)
pi=4*Atn(1)
a=0: b=1: n=10
For i=0 To n-1
coef(i)=Cos(Cos(pi/n*(i+1/2))*(b-a)/2+(b+a)/2)
Next
For i=0 To n-1
w=0
For j=0 To n-1
w=w+coef(j)*Cos(pi/n*i*(j+1/2))
Next
cheby(i)=w*2/n
If cheby(i)<=0 Then t="" Else t=" "
WScript.StdOut.WriteLine i&" : "&t&cheby(i)
Next
- Output:
0 : 1,64716947539031 1 : -0,232299371615172 2 : -5,37151146220477E-02 3 : 2,45823526698163E-03 4 : 2,82119057433916E-04 5 : -7,72222915563781E-06 6 : -5,89855645105608E-07 7 : 1,15214275009379E-08 8 : 6,59629917354465E-10 9 : -1,0022016549982E-11
Visual Basic .NET
Module Module1
Structure ChebyshevApprox
Public ReadOnly coeffs As List(Of Double)
Public ReadOnly domain As Tuple(Of Double, Double)
Public Sub New(func As Func(Of Double, Double), n As Integer, domain As Tuple(Of Double, Double))
coeffs = ChebCoef(func, n, domain)
Me.domain = domain
End Sub
Public Function Eval(x As Double) As Double
Return ChebEval(coeffs, domain, x)
End Function
End Structure
Function AffineRemap(from As Tuple(Of Double, Double), x As Double, t0 As Tuple(Of Double, Double)) As Double
Return t0.Item1 + (x - from.Item1) * (t0.Item2 - t0.Item1) / (from.Item2 - from.Item1)
End Function
Function ChebCoef(fVals As List(Of Double)) As List(Of Double)
Dim n = fVals.Count
Dim theta = Math.PI / n
Dim retval As New List(Of Double)
For i = 1 To n
retval.Add(0.0)
Next
For i = 1 To n
Dim ii = i - 1
Dim f = fVals(ii) * 2.0 / n
Dim phi = (ii + 0.5) * theta
Dim c1 = Math.Cos(phi)
Dim s1 = Math.Sin(phi)
Dim c = 1.0
Dim s = 0.0
For j = 1 To n
Dim jj = j - 1
retval(jj) += f * c
' update c -> cos(j*phi) for next value of j
Dim cNext = c * c1 - s * s1
s = c * s1 + s * c1
c = cNext
Next
Next
Return retval
End Function
Function ChebCoef(func As Func(Of Double, Double), n As Integer, domain As Tuple(Of Double, Double)) As List(Of Double)
Dim Remap As Func(Of Double, Double)
Remap = Function(x As Double)
Return AffineRemap(Tuple.Create(-1.0, 1.0), x, domain)
End Function
Dim theta = Math.PI / n
Dim fVals As New List(Of Double)
For i = 1 To n
fVals.Add(0.0)
Next
For i = 1 To n
Dim ii = i - 1
fVals(ii) = func(Remap(Math.Cos((ii + 0.5) * theta)))
Next
Return ChebCoef(fVals)
End Function
Function ChebEval(coef As List(Of Double), x As Double) As Double
Dim a = 1.0
Dim b = x
Dim c As Double
Dim retval = 0.5 * coef(0) + b * coef(1)
Dim it = coef.GetEnumerator
it.MoveNext()
it.MoveNext()
While it.MoveNext
Dim pc = it.Current
c = 2.0 * b * x - a
retval += pc * c
a = b
b = c
End While
Return retval
End Function
Function ChebEval(coef As List(Of Double), domain As Tuple(Of Double, Double), x As Double) As Double
Return ChebEval(coef, AffineRemap(domain, x, Tuple.Create(-1.0, 1.0)))
End Function
Sub Main()
Dim N = 10
Dim fApprox As New ChebyshevApprox(AddressOf Math.Cos, N, Tuple.Create(0.0, 1.0))
Console.WriteLine("Coefficients: ")
For Each c In fApprox.coeffs
Console.WriteLine(vbTab + "{0: 0.00000000000000;-0.00000000000000;zero}", c)
Next
Console.WriteLine(vbNewLine + "Approximation:" + vbNewLine + " x func(x) approx diff")
Dim nX = 20.0
Dim min = 0.0
Dim max = 1.0
For i = 1 To nX
Dim x = AffineRemap(Tuple.Create(0.0, nX), i, Tuple.Create(min, max))
Dim f = Math.Cos(x)
Dim approx = fApprox.Eval(x)
Console.WriteLine("{0:0.000} {1:0.00000000000000} {2:0.00000000000000} {3:E}", x, f, approx, approx - f)
Next
End Sub
End Module
- Output:
Coefficients: 1.64716947539031 -0.23229937161517 -0.05371511462205 0.00245823526698 0.00028211905743 -0.00000772222916 -0.00000058985565 0.00000001152143 0.00000000065963 -0.00000000001002 Approximation: x func(x) approx diff 0.050 0.99875026039497 0.99875026039487 -9.370282E-014 0.100 0.99500416527803 0.99500416527849 4.622969E-013 0.150 0.98877107793604 0.98877107793600 -4.662937E-014 0.200 0.98006657784124 0.98006657784078 -4.604095E-013 0.250 0.96891242171065 0.96891242171041 -2.322587E-013 0.300 0.95533648912561 0.95533648912587 2.609024E-013 0.350 0.93937271284738 0.93937271284784 4.606315E-013 0.400 0.92106099400289 0.92106099400308 1.980638E-013 0.450 0.90044710235268 0.90044710235243 -2.473577E-013 0.500 0.87758256189037 0.87758256188991 -4.586331E-013 0.550 0.85252452205951 0.85252452205926 -2.461364E-013 0.600 0.82533561490968 0.82533561490988 1.961764E-013 0.650 0.79608379854906 0.79608379854951 4.536371E-013 0.700 0.76484218728449 0.76484218728474 2.553513E-013 0.750 0.73168886887382 0.73168886887359 -2.267075E-013 0.800 0.69670670934717 0.69670670934672 -4.467537E-013 0.850 0.65998314588498 0.65998314588494 -4.485301E-014 0.900 0.62160996827066 0.62160996827111 4.444223E-013 0.950 0.58168308946388 0.58168308946379 -8.992806E-014 1.000 0.54030230586814 0.54030230586859 4.468648E-013
Wren
import "/fmt" for Fmt
var mapRange = Fn.new { |x, min, max, minTo, maxTo| (x - min)/(max - min)*(maxTo - minTo) + minTo }
var chebCoeffs = Fn.new { |func, n, min, max|
var coeffs = List.filled(n, 0)
for (i in 0...n) {
var f = func.call(mapRange.call((Num.pi * (i + 0.5) / n).cos, -1, 1, min, max)) * 2 / n
for (j in 0...n) coeffs[j] = coeffs[j] + f * (Num.pi * j * (i + 0.5) / n).cos
}
return coeffs
}
var chebApprox = Fn.new { |x, n, min, max, coeffs|
if (n < 2 || coeffs.count < 2) Fiber.abort("'n' can't be less than 2.")
var a = 1
var b = mapRange.call(x, min, max, -1, 1)
var res = coeffs[0]/2 + coeffs[1]*b
var xx = 2 * b
var i = 2
while (i < n) {
var c = xx*b - a
res = res + coeffs[i]*c
a = b
b = c
i = i + 1
}
return res
}
var n = 10
var min = 0
var max = 1
var coeffs = chebCoeffs.call(Fn.new { |x| x.cos }, n, min, max)
System.print("Coefficients:")
for (coeff in coeffs) Fmt.print("$0s$1.15f", (coeff >= 0) ? " " : "", coeff)
System.print("\nApproximations:\n x func(x) approx diff")
for (i in 0..20) {
var x = mapRange.call(i, 0, 20, min, max)
var f = x.cos
var approx = chebApprox.call(x, n, min, max, coeffs)
var diff = approx - f
var diffStr = diff.toString
var e = diffStr[-4..-1]
diffStr = diffStr[0..-5]
diffStr = (diff >= 0) ? " " + diffStr[0..3] : diffStr[0..4]
Fmt.print("$1.3f $1.8f $1.8f $s", x, f, approx, diffStr + e)
}
- Output:
Coefficients: 1.64716947539031 -0.23229937161517 -0.05371511462205 0.00245823526698 0.00028211905743 -0.00000772222916 -0.00000058985565 0.00000001152143 0.00000000065963 -0.00000000001002 Approximations: x func(x) approx diff 0.000 1.00000000 1.00000000 4.68e-13 0.050 0.99875026 0.99875026 -9.35e-14 0.100 0.99500417 0.99500417 4.61e-13 0.150 0.98877108 0.98877108 -4.72e-14 0.200 0.98006658 0.98006658 -4.60e-13 0.250 0.96891242 0.96891242 -2.31e-13 0.300 0.95533649 0.95533649 2.61e-13 0.350 0.93937271 0.93937271 4.61e-13 0.400 0.92106099 0.92106099 1.98e-13 0.450 0.90044710 0.90044710 -2.47e-13 0.500 0.87758256 0.87758256 -4.58e-13 0.550 0.85252452 0.85252452 -2.46e-13 0.600 0.82533561 0.82533561 1.95e-13 0.650 0.79608380 0.79608380 4.52e-13 0.700 0.76484219 0.76484219 2.54e-13 0.750 0.73168887 0.73168887 -2.27e-13 0.800 0.69670671 0.69670671 -4.47e-13 0.850 0.65998315 0.65998315 -4.37e-14 0.900 0.62160997 0.62160997 4.45e-13 0.950 0.58168309 0.58168309 -8.99e-14 1.000 0.54030231 0.54030231 4.47e-13
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();
- 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)