HW Greg Dannay

Project.toml

@@ -3,6 +3,17 @@ uuid = "ce1ef771-b552-5ca6-80e4-7358b8c578e2"
 authors = ["Florian Oswald <[email protected]>"]
 version = "0.1.0"
 
+[deps]
+ApproXD = "f0e97480-066f-5960-9ba5-e027352bc8af"
+ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
+BasisMatrices = "08854c51-b66b-5062-a90d-8e7ae4547a49"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

src/HWfuncapp.jl

@@ -6,37 +6,106 @@ using FastGaussQuadrature # to get chebyshevnodes
 # you will see me often qualify code with `PyPlot.plot` or so.
 using PyPlot  
 import ApproXD: getBasis, BSpline
+using LinearAlgebra
 using Distributions
 using ApproxFun
 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)
+export q1, q2, q3, q4, q5, q6, q7, runall
+
+function q1(n=15)
 	f(x) = x .+ 2x.^2 - exp.(-x)
+	n_new = 100
+
+	deg = n:-1:1
+	points = range(-3, 3, length=n_new)
+
+	node = 3 * cos.((2 .* deg .- 1) .* pi ./ (2 .* n))
+	y = f(node)
 
+	cheb = zeros(n, n)
+	for i in 1:n
+		cheb[:, i] = ChebyT.(unitmap(node, -3, 3), i - 1)
+	end
+	c = cheb^-1 * y
+
+	estim = zeros(n_new)
+	for i in 1:n
+		estim += c[i] * ChebyT.(unitmap(points, -3, 3), i - 1)
+	end
+
+	subplot(122)
+	PyPlot.plot(points, f.(points), label="True function")
+	PyPlot.scatter(points, estim, color="red", s=2, label="approximation")
+	xlabel("x")
+	ylabel("y")
+	title("Chebyshev approximation")
+	legend()
+	subplot(121)
+	PyPlot.plot(points, f.(points) - estim)
+	title("Error")
 	# without using PyPlot, just erase the `PyPlot.` part
 	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png"))
+
+	err = sum(abs.(f.(points) - estim))
 	return Dict(:error=>maximum(abs,err))
 end
 
 function q2(b::Number)
 	@assert b > 0
-	# use ApproxFun.jl to do the same:
+	f(x) = x .+ 2x.^2 - exp.(-x)
 
-	Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q2.png"))
+	precision = 100
+	points = range(-3, 3, length=precision)
+	# use ApproxFun.jl to do the same:
+	x = Fun(f, Chebyshev(-b..b))
+	estim = x.(points)
+
+	subplot(122)
+	PyPlot.plot(points, f.(points), label="True function")
+	PyPlot.scatter(points, estim, color="red", s=2, label="approximation")
+	xlabel("x")
+	ylabel("y")
+	title("Chebyshev approximation")
+	legend()
+	subplot(121)
+	PyPlot.plot(points, f.(points) - estim)
+	title("Error")
+
+	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q2.png"))
 end
 
 function q3(b::Number)
+	x = Fun(identity, -b..b)
+	f = sin(x^2)
+	g = cos(x)
+	h = f - g
 
+	r = roots(h)
+
+	p = Plots.plot(h, labels="h")
+	Plots.scatter!(r, h.(r), labels="roots.h")
+	Plots.savefig(joinpath(dirname(@__FILE__),"..","q3.png"))
 	# p is your plot
+	g = cumsum(h)
+	integral = g(0) - g(-b)
 	return (p,integral)
 end
 
 # optinal
 function q4()
-
+	precision = 10000
+	points = range(-1, stop = 1, length=precision)
+	fig = PyPlot.plot(points, ChebyT.(points, 0), label=0)
+	for i in 1:8
+		fig = PyPlot.plot(points, ChebyT.(points, i), label=i)
+	end
+	title("Chebyshev basis")
+	legend()
+	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q4.png"))
 	return fig
 end
 
@@ -56,7 +125,7 @@ mutable struct ChebyType
 	function ChebyType(_nodes::Union{Vector,LinRange},_deg,_lb,_ub,_f::Function)
 		n = length(_nodes)
 		y = _f(_nodes)
-		_basis = Float64[ChebyT(unitmap(_nodes[i],_lb,_ub),j) for i=1:n,j=0:_deg]
+		_basis = Float64[ChebyT(unitmap(_nodes[i],_lb,_ub),j) for i=1:n,j=0:_deg - 1]
 		_coefs = _basis \ y  # type `?\` to find out more about the backslash operator. depending the args given, it performs a different operation
 		# create a ChebyComparer with those values
 		new(_f,_nodes,_basis,_coefs,_deg,_lb,_ub)
@@ -67,8 +136,8 @@ end
 function predict(Ch::ChebyType,x_new)
 
 	true_new = Ch.f(x_new)
-	basis_new = Float64[ChebyT(unitmap(x_new[i],Ch.lb,Ch.ub),j) for i=1:length(x_new),j=0:Ch.deg]
-	basis_nodes = Float64[ChebyT(unitmap(Ch.nodes[i],Ch.lb,Ch.ub),j) for i=1:length(Ch.nodes),j=0:Ch.deg]
+	basis_new = Float64[ChebyT(unitmap(x_new[i],Ch.lb,Ch.ub),j) for i=1:length(x_new),j=0:Ch.deg - 1]
+	basis_nodes = Float64[ChebyT(unitmap(Ch.nodes[i],Ch.lb,Ch.ub),j) for i=1:length(Ch.nodes),j=0:Ch.deg - 1]
 	preds = basis_new * Ch.coefs
 	preds_nodes = basis_nodes * Ch.coefs
 
@@ -78,25 +147,107 @@ end
 function q5(deg=(5,9,15),lb=-1.0,ub=1.0)
 
 	runge(x) = 1.0 ./ (1 .+ 25 .* x.^2)
+	f(x) = x .+ 2x.^2 - exp.(-x)
+	precision = 10000
+	points = range(lb, ub, length=precision)
+	cheb_m = zeros(precision, length(deg))
+	uni_m = zeros(precision, length(deg))
+
+	for (i, d) in enumerate(deg)
+		degse = d:-1:1
+		node = 1/2 .* (ub .+ lb) .+ 1/2 .* (ub - lb) .* cos.((2 .* degse .- 1) .* pi ./ (2 .* d))
+		uni = LinRange(lb, ub, d)
+		cheb = ChebyType(node, d, lb, ub, runge)
+		unic = ChebyType(uni, d, lb, ub, runge)
+
+		pred_cheb = predict(cheb, points)
+		pred_uni = predict(unic, points)
+
+		cheb_m[:, i] = pred_cheb["preds"]
+		uni_m[:, i] = pred_uni["preds"]
+	end
 
-
+	subplot(121)
+	PyPlot.plot(points, runge.(points), label="True")
+	for (i, d) in enumerate(deg)
+		PyPlot.plot(points, cheb_m[:, i,], label=d)
+	end
+	title("Chebyshev nodes")
+	legend()
+	subplot(122)
+	PyPlot.plot(points, runge.(points), label="True")
+	for (i, d) in enumerate(deg)
+		PyPlot.plot(points, uni_m[:, i,], label=d)
+	end
+	title("Uniform nodes")
+	legend()
 	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png"))
 
 end
 
 
 
 function q6()
-
+	f(x) = abs(x)^0.5
+	precision = 63
+	x = range(-1, 1, length=precision)
+	y = f.(x)
 	# compare 2 knot vectors with runge's function
 
+	myknots = vcat(range(-1,stop = -0.1,length = 5), 0, 0, 0, range(0.1,stop = 1,length =5))
+	bs = BSpline(13, 3, -5, 5)
+	bs2 = BSpline(myknots, 3)
+	B = Array(getBasis(collect(x), bs))
+	B2 = Array(getBasis(collect(x), bs2))
+	c = B \ y
+	c2 = B2 \ y
+	subplot(311)
+	PyPlot.plot(x, y, label="True function")
+	PyPlot.scatter(x, B * c, color="red", s=2, label="Approximation")
+	title("Uniform knot vector")
+	legend()
+	subplot(312)
+	PyPlot.plot(x, y, label="True function")
+	PyPlot.scatter(x, B2 * c2, color="red", s=2, label="Approximation")
+	title("Knot multiplicity x=0")
+	legend()
+	subplot(313)
+	PyPlot.plot(x, y - B * c, color="blue", label="Uniform")
+	PyPlot.plot(x, y - B2 * c2, color="red", label="Knot multiplicity")
+	title("Error plots")
+	legend()
 	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q6.png"))
 
 end
 
 function q7()
-
-
+	f(x) = abs(x)^0.5
+	precision = 63
+	x = range(-1, 1, length=precision)
+	myknots = vcat(range(-1,stop = -0.1,length = 5), 0, 0, 0, range(0.1,stop = 1,length =5))
+
+	y = f.(x)
+	bs = BSpline(13, 3, -1, 1)
+	bs2 = BSpline(myknots, 3)
+	B = Array(getBasis(collect(x), bs))
+	B2 = Array(getBasis(collect(x), bs2))
+	c = B \ y
+	c2 = B2 \ y
+	subplot(311)
+	PyPlot.plot(x, y, label="True function")
+	PyPlot.scatter(x, B * c, color="red", s=2, label="Approximation")
+	title("Uniform knot vector")
+	legend()
+	subplot(312)
+	PyPlot.plot(x, y, label="True function")
+	PyPlot.scatter(x, B2 * c2, color="red", s=2, label="Approximation")
+	title("Knot multiplicity x=0")
+	legend()
+	subplot(313)
+	PyPlot.plot(x, y - B * c, color="blue", label="Uniform")
+	PyPlot.plot(x, y - B2 * c2, color="red", label="Knot multiplicity")
+	title("Error plots")
+	legend()
 	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q7.png"))
 end