You can always use ApproxFun.jl
or you can write your own fully automated (textual) Chebyshev function generator for sin function using Symbolics.jl.
using Symbolics
@syms Cos(var1) func(var1)
@variables a b pi x new_x u
###############################
# targetfunc = sin(u) #
# targetlowerbound = 0 #
# targetupperbound = pi/2 #
###############################
funcname = "sine"
targetfunc = sin(u)
targetlowerbound = 0
targetupperbound = pi//2
transformedfunc = substitute( targetfunc, Dict(u => (1//2)*(a+b-a*x+b*x)) )
transformedfunc = substitute( transformedfunc, Dict(a => targetlowerbound, b => targetupperbound) )
# myfunc = simplify(myfunc,expand=true)
println("Stage 1:")
println("transformedfunc is ",transformedfunc)
println()
# Now use Symbolics.jl to automatically built myfunc
myfunc_expr=build_function(transformedfunc,x, expression = Val{false})
myfunc=eval(myfunc_expr)
X(k,n) = Cos((2*k+1)//(2*n+2)*pi)
function c(j,n)
if j == 0
return (1//(n+1)) * sum(k->func(X(k,n)),0:1:n)
else
return (2//(n+1)) * sum(k->func(X(k,n))*Cos(j * (2*k+1)//(2*n+2) * pi),0:1:n)
end
end
function T(j,z)
if j == 0
return 1
elseif j == 1
return z
else
return 2 * z * T(j-1,z) - T(j-2,z)
end
end
P(n) = sum(j->c(j,n)*T(j,x),0:1:n)
# Get the 13th Degree Chebyshev Approximation
soln = P(13)
soln = substitute( soln, Dict(Cos => cos,func => myfunc, pi => Base.pi) )
soln = simplify(soln,expand=true)
println("Stage 2: Intermediate Solution")
println("soln is ",soln)
println()
# Normalise soln to final_soln by scaling the x-axis
final_soln = substitute( soln, Dict( x => (-1 + (new_x-a)//(((a+b)//2)-a)) ) )
final_soln = substitute( final_soln, Dict(a => targetlowerbound, b => targetupperbound) )
final_soln = substitute( final_soln, Dict( pi => Base.pi ) )
final_soln = simplify(final_soln,expand=true)
# replace new_x variable back to plain old x variable
final_soln = substitute( final_soln, Dict(new_x => x) )
println("Stage 3: Final Solution")
println("final_soln is ",final_soln)
println()
final_soln_str = ""
firsttime = true
println("Summary of terms:")
for k=0:1:Symbolics.degree(final_soln)
global final_soln_str
# global final_solnrat
global firsttime
if !firsttime
final_soln_str = final_soln_str * " + "
else
firsttime = false
end
println("The $k term is ",Float64(Symbolics.coeff(final_soln,x^k))," * x^$k")
final_soln_str = final_soln_str * string(Float64(Symbolics.coeff(final_soln,x^k))) * " * x^$k"
end
println()
println("Explicit $(funcname) function is")
println("$(funcname)(x) = ",final_soln_str)
println()
println("Horner method $(funcname) function is")
println("HornerArray = zeros($(Symbolics.degree(final_soln) + 1))")
for k=Symbolics.degree(final_soln):-1:0
println("HornerArray[$(Symbolics.degree(final_soln)+1-k)]= ",Symbolics.coeff(final_soln,x^k))
end
println()
print(
"""
function HornerFunction(x)
global HornerArray
result = 0.0
for value in HornerArray
result = result * x + value
end
return result
end
"""
)
with the output
Stage 1:
transformedfunc is sin((1//2)*((1//2)*pi + (1//2)*pi*x))
Stage 2: Intermediate Solution
soln is 0.707106781186548 + 0.5553603672697928x - 0.21808950623873247(x^2) - 0.057095699218664704(x^3) + 0.011210714326173552(x^4) + 0.0017609748880502549(x^5) - 0.00023051107399396642(x^6) - 2.5863280733824062e-5(x^7) + 2.5391214367890106e-6(x^8) + 2.2157795116178087e-7(x^9) - 1.7401841237837546e-8(x^10) - 1.2406837462159664e-9(x^11) + 8.040907622281727e-11(x^12) + 4.287617879786662e-12(x^13)
Stage 3: Final Solution
final_soln is 2.1312365392454902e-16 + 0.9999999999999833x + 9.802328219954415e-13(x^2) - 0.16666666668431543(x^3) + 1.47824161765564e-10(x^4) + 0.00833333263088191(x^5) + 2.101024760311869e-9(x^6) - 0.00019841689126759302(x^7) + 5.761480511332656e-9(x^8) + 2.7502155501620215e-6(x^9) + 3.6549238210666234e-9(x^10) - 2.667607157014283e-8(x^11) + 4.4780748032530374e-10(x^12) + 9.909405845154378e-11(x^13)
Summary of terms:
The 0 term is 2.1312365392454902e-16 * x^0
The 1 term is 0.9999999999999833 * x^1
The 2 term is 9.802328219954415e-13 * x^2
The 3 term is -0.16666666668431543 * x^3
The 4 term is 1.47824161765564e-10 * x^4
The 5 term is 0.00833333263088191 * x^5
The 6 term is 2.101024760311869e-9 * x^6
The 7 term is -0.00019841689126759302 * x^7
The 8 term is 5.761480511332656e-9 * x^8
The 9 term is 2.7502155501620215e-6 * x^9
The 10 term is 3.6549238210666234e-9 * x^10
The 11 term is -2.667607157014283e-8 * x^11
The 12 term is 4.4780748032530374e-10 * x^12
The 13 term is 9.909405845154378e-11 * x^13
Explicit sine function is
sine(x) = 2.1312365392454902e-16 * x^0 + 0.9999999999999833 * x^1 + 9.802328219954415e-13 * x^2 + -0.16666666668431543 * x^3 + 1.47824161765564e-10 * x^4 + 0.00833333263088191 * x^5 + 2.101024760311869e-9 * x^6 + -0.00019841689126759302 * x^7 + 5.761480511332656e-9 * x^8 + 2.7502155501620215e-6 * x^9 + 3.6549238210666234e-9 * x^10 + -2.667607157014283e-8 * x^11 + 4.4780748032530374e-10 * x^12 + 9.909405845154378e-11 * x^13
Horner method sine function is
HornerArray = zeros(14)
HornerArray[1]= 9.909405845154378e-11
HornerArray[2]= 4.4780748032530374e-10
HornerArray[3]= -2.667607157014283e-8
HornerArray[4]= 3.6549238210666234e-9
HornerArray[5]= 2.7502155501620215e-6
HornerArray[6]= 5.761480511332656e-9
HornerArray[7]= -0.00019841689126759302
HornerArray[8]= 2.101024760311869e-9
HornerArray[9]= 0.00833333263088191
HornerArray[10]= 1.47824161765564e-10
HornerArray[11]= -0.16666666668431543
HornerArray[12]= 9.802328219954415e-13
HornerArray[13]= 0.9999999999999833
HornerArray[14]= 2.1312365392454902e-16
function HornerFunction(x)
global HornerArray
result = 0.0
for value in HornerArray
result = result * x + value
end
return result
end
julia>
using the formula
and the transformation
