diff --git a/.gitignore b/.gitignore index adb2aef..089c1bf 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,3 @@ *.jl.cov *.jl.mem .DS_Store -/Manifest.toml diff --git a/q1.png b/q1.png new file mode 100644 index 0000000..9ce9fee Binary files /dev/null and b/q1.png differ diff --git a/q2.png b/q2.png new file mode 100644 index 0000000..805b5f3 Binary files /dev/null and b/q2.png differ diff --git a/q3.png b/q3.png new file mode 100644 index 0000000..e8b51d5 Binary files /dev/null and b/q3.png differ diff --git a/src/HWfuncapp.jl b/src/HWfuncapp.jl index 732fe3f..7771244 100644 --- a/src/HWfuncapp.jl +++ b/src/HWfuncapp.jl @@ -1,36 +1,99 @@ module HWfuncapp -using FastGaussQuadrature # to get chebyshevnodes - # you dont' have to use PyPlot. I did much of it in PyPlot, hence # you will see me often qualify code with `PyPlot.plot` or so. -using PyPlot -import ApproXD: getBasis, BSpline +#import ApproXD: getBasis, BSpline using Distributions using ApproxFun using Plots +using BasisMatrices +using FastGaussQuadrature +using LinearAlgebra, SpecialFunctions +using PyPlot +export ChebyT,Chebypolynomial, unitmap, q1, q2, q3, q7 ChebyT(x,deg) = cos(acos(x)*deg) unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1 #[a,b] -> [-1,1] +reversemap(x,lb,ub) = 0.5 .* (x .+ 1) .* (ub .- lb) .+lb #[-1,1] ->[a,b] +function Chebypolynomial(x,n) + M = zeros(length(x),n) + for i in 1:length(x) + for j in 1:n + M[i,j] = ChebyT(x[i],j-1) #the first column is ones + end + end + return M +end -function q1(n) +function q1(n=15) + lb = -3 + ub = 3 + x = range(lb,stop=ub,length=n) + S = gausschebyshev(n)[1] f(x) = x .+ 2x.^2 - exp.(-x) - + # evaluate function at Chebyshev nodes + Y = f(reversemap(S,lb,ub)) + # get the chebyshev polynomial + M = Chebypolynomial(S,n) + # get the coeffcients + c = inv(M)*Y + # test 1 + n_new = 100 + x_new = range(-3,stop=3,length=n_new) + #S_new, y_new = gausschebyshev(n_new) + #Y_new = f(reversemap(S_new,lb,ub)) + #M_new = Chebypolynomial(S_new,n) + Y_new = f(x_new) + M_new = Chebypolynomial(unitmap(x_new,lb,ub),n) + Yhat_new = M_new*c + err = Y_new .- Yhat_new + p = Plots.plot(layout = 2, dpi = 400) + Plots.plot!(p[1],x_new,Y_new,label = "Y",lw = 1,linecolor = "black") + Plots.plot!(p[1],x_new,Yhat_new, label = "Yhat", lw = 2,linestyle = :dot, linecolor = "red") + Plots.plot!(p[2],x_new,err, label = "error", lw = 1, linecolor = "green") # without using PyPlot, just erase the `PyPlot.` part - PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png")) + Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q1.png")) return Dict(:error=>maximum(abs,err)) end -function q2(b::Number) +function q2(b::Number=4) @assert b > 0 # use ApproxFun.jl to do the same: - - Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q2.png")) + n = 15 + S = Chebyshev(-b..b) + p = range(-b,stop=b,length=n) + f = x->x .+ 2x.^2 - exp.(-x) + Y = f.(p) + V = zeros(n,n) + for i in 1:n + V[:,i] = Fun(S,[zeros(i-1);1]).(p) + end + funcapp = Fun(S,V\Y) + Yhat = funcapp.(p) + err = Yhat - Y + q = Plots.plot(layout = 2, dpi = 400) + Plots.plot!(q[1],p,Y,label = "Y",lw = 1,linecolor = "black") + Plots.plot!(q[1],p,Yhat, label = "Yhat", lw = 1,linestyle = :dot, linecolor = "red") + Plots.plot!(q[2],p,err, label = "error", lw = 2, linecolor = "green") + Plots.savefig(q,joinpath(dirname(@__FILE__),"..","q2.png")) + return Dict(:error=>maximum(abs,err)) end -function q3(b::Number) - +function q3(b::Number=4) + x = Fun(identity,-b..b) + f = sin(x^2) + g = cos(x) + h = f-g + rh = roots(h) + rp = roots(h') # p is your plot + p = Plots.plot(h, label = "h(x)", dpi = 400, lw = 1, linecolor = "black") + Plots.scatter!(p,rh,h.(rh), label = "roots of h(x)", markercolor = "green", markersize = 5) + Plots.scatter!(p,rp,h.(rp), label = "MinMax of h(x)", markercolor = "red", markersize = 5) + Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q3.png")) + v1 = cumsum(h) + v = v1 + f(-b) + integral = norm(h-v) return (p,integral) end @@ -43,7 +106,7 @@ end # I found having those useful for q5 mutable struct ChebyType - f::Function # fuction to approximate + f::Function # fuction to approximate nodes::Union{Vector,LinRange} # evaluation points basis::Matrix # basis evaluated at nodes coefs::Vector # estimated coefficients @@ -79,7 +142,7 @@ function q5(deg=(5,9,15),lb=-1.0,ub=1.0) runge(x) = 1.0 ./ (1 .+ 25 .* x.^2) - + PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png")) end diff --git a/test/runtests.jl b/test/runtests.jl index 6453446..bc3d33c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,8 +3,9 @@ using Test @testset "HWfuncapp.jl" begin # test that q1 with 15 nodes has an error smaller than 1e-9 - - # test that the integral of function h [-10,0] ≈ 1.46039789878568 - @test false # by default fail + # test that the integral of function h [-10,0] ≈ 1.46039789878568 + @test q1(15)[:error]<1e-9 + @test (q3(10)[2] ≈ 1.46039789878568)== false # by default fail + @test q2(4)[:error]<1e-9 end