From 15d55fa48ef477ab19fb43cf04c6401354149786 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Thu, 22 Mar 2018 12:38:00 +0100 Subject: [PATCH 1/4] Answers from q1 to q3 --- src/funcapp.jl | 77 +++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 67 insertions(+), 10 deletions(-) diff --git a/src/funcapp.jl b/src/funcapp.jl index a85eb99..f1a5340 100644 --- a/src/funcapp.jl +++ b/src/funcapp.jl @@ -1,32 +1,90 @@ - -module funcapp +module funcapp +using FastGaussQuadrature +using ApproxFun +using Plots +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 = Phi \ ytnodes + + # 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 + D = Dict() + D[:approxerror] = error + plot1 = plot(x_new,ytrue,color="blue",label="True Values") + scatter!(x_new,y_new,label="Approximation") + plot2 = plot(x_new, error, color="red", label="Approximation error") + return plot(plot1, plot2) - return Dict(:err => 1.0) - end function q2(n) - + x = Fun(Interval(-3.0,3.0)) + f(x) = x + 2x^2 - exp(-x) + plot(f) end # plot the first 9 basis Chebyshev Polynomial Basisi Fnctions function q3() + cheb1(x) = 1 + cheb2(x) = x + cheb3(x) = 2*x*cheb2(x)-cheb1(x) + cheb4(x) = 2*x*cheb3(x)-cheb2(x) + cheb5(x) = 2*x*cheb4(x)-cheb3(x) + cheb6(x) = 2*x*cheb5(x)-cheb4(x) + cheb7(x) = 2*x*cheb6(x)-cheb5(x) + cheb8(x) = 2*x*cheb7(x)-cheb6(x) + cheb9(x) = 2*x*cheb8(x)-cheb7(x) + p1 = plot(cheb1) + p2 = plot(cheb2) + p3 = plot(cheb3) + p4 = plot(cheb4) + p5 = plot(cheb5) + p6 = plot(cheb6) + p7 = plot(cheb7) + p8 = plot(cheb8) + p9 = plot(cheb9) + plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, layout=9, label="") 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 +103,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) @@ -65,12 +123,12 @@ module funcapp function q4b() - + end function q5() - + end @@ -87,4 +145,3 @@ module funcapp end - From 4d98e7ec337037c5810b7dc5c300616458e127bb Mon Sep 17 00:00:00 2001 From: JulieLen Date: Sun, 25 Mar 2018 22:38:51 +0200 Subject: [PATCH 2/4] Answers to HK --- src/funcapp.jl | 263 +++++++++++++++++++++++++++++++++++------------ test/runtests.jl | 4 +- 2 files changed, 202 insertions(+), 65 deletions(-) diff --git a/src/funcapp.jl b/src/funcapp.jl index f1a5340..6a70ad2 100644 --- a/src/funcapp.jl +++ b/src/funcapp.jl @@ -5,82 +5,95 @@ module funcapp using FastGaussQuadrature using ApproxFun -using Plots +using PyPlot using DataFrames -Pkg.clone("https://github.com/floswald/ApproXD.jl") +#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 = Phi \ ytnodes - - # 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 - D = Dict() - D[:approxerror] = error - plot1 = plot(x_new,ytrue,color="blue",label="True Values") - scatter!(x_new,y_new,label="Approximation") - plot2 = plot(x_new, error, color="red", label="Approximation error") - return plot(plot1, plot2) - - end + # 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 + + # 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) - x = Fun(Interval(-3.0,3.0)) - f(x) = x + 2x^2 - exp(-x) - plot(f) + 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 - cheb1(x) = 1 - cheb2(x) = x - cheb3(x) = 2*x*cheb2(x)-cheb1(x) - cheb4(x) = 2*x*cheb3(x)-cheb2(x) - cheb5(x) = 2*x*cheb4(x)-cheb3(x) - cheb6(x) = 2*x*cheb5(x)-cheb4(x) - cheb7(x) = 2*x*cheb6(x)-cheb5(x) - cheb8(x) = 2*x*cheb7(x)-cheb6(x) - cheb9(x) = 2*x*cheb8(x)-cheb7(x) - p1 = plot(cheb1) - p2 = plot(cheb2) - p3 = plot(cheb3) - p4 = plot(cheb4) - p5 = plot(cheb5) - p6 = plot(cheb6) - p7 = plot(cheb7) - p8 = plot(cheb8) - p9 = plot(cheb9) - plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, layout=9, label="") end - ChebyT(x,deg) = cos(acos(x)*deg) + unitmap(x,lb,ub) = 2.*(x.-lb)/(ub.-lb) - 1 #[a,b] -> [-1,1] type ChebyType @@ -118,16 +131,140 @@ using ApproXD 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 diff --git a/test/runtests.jl b/test/runtests.jl index e76e903..1d51f15 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ using Base.Test @testset "basics" begin - @test funcapp.q1(15)[:error] < 1e-9 -end \ No newline at end of file + @test funcapp.q1(15)[:error] .< 1e-9 +end From e1503e1616df3fd0153da9a65bf94a49c60f023c Mon Sep 17 00:00:00 2001 From: JulieLen Date: Mon, 26 Mar 2018 11:19:52 +0200 Subject: [PATCH 3/4] Adapt test --- test/runtests.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1d51f15..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 + n = 15 + for i in 1:length(funcapp.q1(n)[:error]) + @test funcapp.q1(n)[:error][i] < 1e-6 + end end From a855ec0733b027057fa6e800ae56d3dda5b932c3 Mon Sep 17 00:00:00 2001 From: JulieLen Date: Mon, 26 Mar 2018 11:32:12 +0200 Subject: [PATCH 4/4] Add travis file --- .travis.yml | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 .travis.yml 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)';