diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..7a3fe23 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,10 @@ +language: julia +os: + - linux +julia: + - 0.6 +notifications: + email: false +#script: # the default script is equivalent to the following +# - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi +# - julia -e 'Pkg.clone(pwd()); Pkg.build("Example"); Pkg.test("Example"; coverage=true)'; diff --git a/src/funcapp.jl b/src/funcapp.jl index a85eb99..6a70ad2 100644 --- a/src/funcapp.jl +++ b/src/funcapp.jl @@ -1,32 +1,103 @@ - -module funcapp +module funcapp +using FastGaussQuadrature +using ApproxFun +using PyPlot +using DataFrames +#Pkg.clone("https://github.com/floswald/ApproXD.jl") +using ApproXD # use chebyshev to interpolate this: function q1(n) + # Define the function + f(x) = x + 2x^2 - exp(-x) + # Specify degree of approximation + deg = n-1 + # Define bounds of domain + ub = 3.0 + lb = -3.0 + # Get interpolation nodes + nodes = gausschebyshev(n)[1] + tnodes = (nodes+1)*(ub-lb)/2 + lb + # Get function value at x, y=f(x) + ytnodes = map(f, tnodes) + # Evaluate the chebyshev basis matrix at Z + basmat = [cos.(acos.(nodes[i])*j) for i=1:n,j=0:(deg)] + #Compute approx coeffs c by matrix inversion + coef = basmat \ ytnodes - return Dict(:err => 1.0) - - end + # Evaluation at new points + + n_new = 100 + x_new = linspace(lb,ub,n_new) + ## Evaluate the basis matrix at x_new + ## We need to map x_new to -1,1 to find the base matrix + x_new_t = 2(x_new-lb)/(ub-lb) -1 + basmat_new = [cos.(acos.(x_new_t[i]) * j) for i=1:n_new,j=0:(deg)] + y_new = basmat_new * coef + ytrue = map(f, x_new) + error = y_new - ytrue + figure("Question 1", figsize=(10, 8)) + subplot(121) + plot(x_new,ytrue,color="blue",label="True Values") + scatter(x_new,y_new,label="Approximation") + title("Plain Vanilla Chebyshev approximation") + subplot(122) + plot(x_new, error, color="red", label="Approximation error") + return(DataFrame(error = abs.(error))) +end function q2(n) - + ff = x -> x + 2*x.^2 - exp.(-x) + F = Fun(ff,Chebyshev(Interval(-3.0,3.0))) + x = collect(linspace(-3.0, 3.0, n)) + y = [] + for i in 1:n + push!(y, F(x[i])) + end + error = ff(x) - y + + figure("Question 2", figsize = (10,8)) + subplot(121) + scatter(x, y, label = "Function Approximation", s = 50, c = "red", marker = "+", alpha=0.5) + plot(x, ff(x), label = "True Function") + legend() + subplot(122) + plot(x, error) + title("Approximation Error") + end + ChebyT(x,deg) = cos.(acos.(x)*deg) # plot the first 9 basis Chebyshev Polynomial Basisi Fnctions function q3() + x = collect(linspace(0,1.0,500)) + + data = [x] + for i in 0:8 + data = hcat(data, [ChebyT(x, i)]) + end + + figure("Question 3", figsize = (8, 10)) + for i in 331:339 + subplot(i) + plot(data[1], data[i-329]) + title("Polynomial of order $(i-331)", fontsize = 10) + xticks(size = 6) + yticks(size = 6) + end end - ChebyT(x,deg) = cos(acos(x)*deg) + unitmap(x,lb,ub) = 2.*(x.-lb)/(ub.-lb) - 1 #[a,b] -> [-1,1] type ChebyType - f::Function # fuction to approximate + f::Function # fuction to approximate nodes::Union{Vector,LinSpace} # evaluation points basis::Matrix # basis evaluated at nodes coefs::Vector # estimated coefficients @@ -45,7 +116,7 @@ module funcapp new(_f,_nodes,_basis,_coefs,_deg,_lb,_ub) end end - + # function to predict points using info stored in ChebyType function predict(Ch::ChebyType,x_new) @@ -60,17 +131,141 @@ module funcapp function q4a(deg=(5,9,15),lb=-1.0,ub=1.0) + h(x) = 1 ./ (1+ 25 .* x .^2) + x_new = linspace(lb, ub, 100) + y_new = h(x_new) + #Chebyshev nodes + + nodes = gausschebyshev(deg[1]+1)[1] + d5 = predict(ChebyType(nodes, deg[1], lb, ub, h), + x_new) + nodes = gausschebyshev(deg[2]+1)[1] + d9 = predict(ChebyType(nodes, deg[2], lb, ub, h), + x_new) + nodes = gausschebyshev(deg[3]+1)[1] + d16 = predict(ChebyType(nodes, deg[3], lb, ub, h), + x_new) + + # Uniform nodes + + grid = collect(linspace(lb, ub, (deg[1]+1))) + u5 = predict(ChebyType(grid, deg[1], lb, ub, h), + x_new) + grid = collect(linspace(lb, ub, (deg[2]+1)) ) + u9 = predict(ChebyType(grid, deg[2], lb, ub, h), + x_new) + grid = collect(linspace(lb, ub, (deg[3]+1))) + u16 = predict(ChebyType(grid, deg[3], lb, ub, h), + x_new) + + #Plots + + figure("Question 4-a", figsize=(10, 8)) + subplot(121) + plot(x_new, y_new, label="True function") + title("Function Approximation\nwith Chebyshev nodes") + plot(d5[:"x"], d5[:"preds"], label = "k=5") + plot(d9[:"x"], d9[:"preds"], label = "k=9") + plot(d16[:"x"], d16[:"preds"], label = "k=16") + subplot(122) + plot(x_new, y_new, label="True function") + title("Function Approximation\nwith Uniform nodes") + plot(u5[:"x"], u5[:"preds"], label = "k=5") + plot(u9[:"x"], u9[:"preds"], label = "k=9") + plot(u16[:"x"], u16[:"preds"], label = "k=16") end function q4b() + #Spaced knots + nknots = 13 + deg = 5 + lb = -5.0 + ub = 5.0 + h(x) = 1 ./ (1 + 25 .* x .^2) + b1 = BSpline(nknots,deg,lb,ub) + nevals = 5 * nknots + evals = collect(linspace(lb, ub, nevals)) + mat1 = full(getBasis(evals,b1)) + coef1 = mat1 \ h(evals) + y1 = mat1 * coef1 + error1 = h(evals) - y1 + + #Knots around zero + + t = gausschebyshev(12)[1] + tt = (t+1)*(ub-lb)/2 + lb + t1 = tt[1:6] + 1 + t2 = tt[7:12] - 1 + T = sort(vcat(t1, 0.0, t2)) + b2 = BSpline(T, deg) + mat2 = full(getBasis(evals, b2)) + coef2 = mat2 \ h(evals) + y2 = mat2 * coef2 + error2 = h(evals) - y2 + + # Plots - + x_new = collect(linspace(lb, ub, 200)) + ytrue = h(x_new) + figure("Question 4-b", figsize=(10, 8)) + subplot(121) + plot(x_new, ytrue) + title("Runge's function") + subplot(122) + plot(evals, error1, label = "spaced knots") + plot(evals, error2, label = "concentrated knots") + scatter(T, zeros(13), label = "knots location", s = 50, c = "green", marker = "x", alpha=0.5 ) + legend() + title("Error in Runge's function") end function q5() - + v(x) = sqrt.(abs.(x)) + lb = -1.0 + ub = 1.0 + nknots = 13 + deg=3 + nevals = 5 * nknots + evals = collect(linspace(lb, ub, nevals)) + + # Uniform Knots + b1 = BSpline(nknots,deg,lb,ub) + mat1 = full(getBasis(evals,b1)) + coef1 = mat1 \ v(evals) + y1 = mat1 * coef1 + error1 = v(evals) - y1 + + # Knots with multiplicity at 0 + + k = vcat(collect(linspace(-1, -0.1, 5)), 0, 0, 0, collect(linspace(0.1, 1, 5)) ) + b2 = BSpline(k, deg) + mat2 = full(getBasis(evals, b2)) + coef2 = mat2 \ v(evals) + y2 = mat2 * coef2 + error2 = v(evals) - y2 + + x_new = collect(linspace(-1.0, 1.0, 200)) + y1_new = full(getBasis(x_new, b2)) * coef2 + y2_new = full(getBasis(x_new, b1)) * coef1 + ytrue = v(x_new) + + figure("Question 5", figsize = (10,10)) + subplot(221) + plot(x_new, ytrue) + title("True function") + subplot(222) + plot(x_new, y1_new, label = "Uniform knots") + plot(x_new, y2_new, label = "Knots with multiplicity at 0") + legend() + title("Function Approximations") + subplot(223) + plot(evals, error1, label = "Uniform knots") + plot(evals, error2, label = "Knots with multiplicity at 0") + legend() + title("Approximation errors") + end @@ -87,4 +282,3 @@ module funcapp end - diff --git a/test/runtests.jl b/test/runtests.jl index e76e903..02c45bb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,8 @@ using Base.Test @testset "basics" begin - @test funcapp.q1(15)[:error] < 1e-9 -end \ No newline at end of file + n = 15 + for i in 1:length(funcapp.q1(n)[:error]) + @test funcapp.q1(n)[:error][i] < 1e-6 + end +end