my homework_Xiaojie

.gitignore

@@ -2,4 +2,4 @@
 *.jl.cov
 *.jl.mem
 .DS_Store
-/Manifest.toml
+

Project.toml

@@ -4,7 +4,17 @@ authors = ["Florian Oswald <[email protected]>"]
 version = "0.1.0"
 
 [extras]
+[deps]
+ApproXD = "f0e97480-066f-5960-9ba5-e027352bc8af"
+ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
+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"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test"]
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

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

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