diff --git a/Project.toml b/Project.toml index 516f264..1dcf6ac 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,13 @@ uuid = "ce1ef771-b552-5ca6-80e4-7358b8c578e2" authors = ["Florian Oswald "] version = "0.1.0" +[deps] +ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" + [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/q1.png b/q1.png new file mode 100644 index 0000000..3df5c8e Binary files /dev/null and b/q1.png differ diff --git a/q2.png b/q2.png new file mode 100644 index 0000000..aa697e5 Binary files /dev/null and b/q2.png differ diff --git a/q3.png b/q3.png new file mode 100644 index 0000000..3d0cdba Binary files /dev/null and b/q3.png differ diff --git a/src/HWfuncapp.jl b/src/HWfuncapp.jl index 732fe3f..db98281 100644 --- a/src/HWfuncapp.jl +++ b/src/HWfuncapp.jl @@ -4,7 +4,7 @@ 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 +using PyPlot import ApproXD: getBasis, BSpline using Distributions using ApproxFun @@ -13,37 +13,77 @@ using Plots ChebyT(x,deg) = cos(acos(x)*deg) unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1 #[a,b] -> [-1,1] -function q1(n) +function q1(n = 15) f(x) = x .+ 2x.^2 - exp.(-x) + deg = n-1 + lb = -3 + ub = 3 + node = gausschebyshev(n)[1] + chebnode = 0.5 .* (lb .+ ub) .+ 0.5 .* (ub .- lb) .* node + y = map(f, chebnode) + basismatrix = [cos((n - i + 0.5) * (j - 1)* π / n) for i = 1:n, j = 1:n] + c = basismatrix \ y + + n_new = 100 + x_2 = range(lb, stop = ub, length = n_new) + base2 = [ChebyT(unitmap(x_2[i], lb, ub), j) for i in 1:n_new, j = 0:deg] + y2 = base2 * c + error = y2 - f(collect(-3:6/99:3)) + Plots.plot(x_2, [f(collect(-3:6/99:3)) error], labels = ["true_val" "error"], title = "Q1", layout = 2) + scatter!(x_2, y2, labels = ["approx"]) # without using PyPlot, just erase the `PyPlot.` part - PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png")) - return Dict(:error=>maximum(abs,err)) + Plots.savefig(joinpath(dirname(@__FILE__),"..","q1.png")) + return Dict(:error=>maximum(abs,error)) end function q2(b::Number) @assert b > 0 # use ApproxFun.jl to do the same: - - Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q2.png")) + b = 4 + n_new = 100 + f(x) = x .+ 2x.^2 - exp.(-x) + f2 = Fun(f,-b..b) + x = range(-b, b, length = n_new) + y2 = f2.(x) + y = f.(x) + error = y - y2 + Plots.plot(x, [y error], labels = ["true_val" "error"], title = "Q2", layout = 2) + scatter!(x, y2, labels = ["approx"]) + + Plots.savefig(joinpath(dirname(@__FILE__),"..","q2.png")) end - +using LinearAlgebra #for cumsum function q3(b::Number) - + b = 10 + x = Fun(identity, -b..b) + f = sin(x^2) + g = cos(x) + h = f - g + r = roots(h) + p = Plots.plot(h, label = "h", title = "Q3") + scatter!(r,h.(r),labels="roots h") + Plots.savefig(joinpath(dirname(@__FILE__),"..","q3.png")) + Plots.savefig(joinpath(p,dirname(@__FILE__),"..","q3.png")) + k1 = cumsum(h) + k = k1 + h(-b) + integral = norm(h-k) # p is your plot return (p,integral) end # optinal -function q4() - +function q4(n) + x = range(0, 1, length = n) + x2 = collect(0:0.01:1) + fig = Plots.plot(x, [ChebyT(x2, deg) for deg in 1:9]) #acos not working + Plots.savefig(joinpath(p,dirname(@__FILE__),"..","q4.png")) return fig 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 +119,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 @@ -113,6 +153,4 @@ function runall() q7() end - - end # module diff --git a/test/runtests.jl b/test/runtests.jl index 6453446..4c9e8d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,8 +3,8 @@ using Test @testset "HWfuncapp.jl" begin # test that q1 with 15 nodes has an error smaller than 1e-9 - + @test q1(15)[:error] < 1e-9 # test that the integral of function h [-10,0] ≈ 1.46039789878568 - + @test q3(10) ≈ 1.46039789878568 @test false # by default fail end