using Plots
using CurveFit
using Polynomials
using Interpolations
using LaTeXStrings
xgiven = collect(0:10)
ygiven = [1,2,3,1,4,6,1,5,10,2,2]
x = LinRange(0,10,1000)
p = poly_fit(xgiven,ygiven,length(xgiven)-1)
fp = Polynomial(p)
ypoly = fp.(x)
f = LinearInterpolation(xgiven, ygiven)
yinterp = f.(x)
plot(x, zeros(length(x)), fillrange=yinterp, color="lightgray")
plot!(x,yinterp, color="red", lw=2)
scatter!(xgiven,ygiven, color="black")
plot!(legend=nothing)
plot!(ylim=[0,13])
p1=plot!(title="linear interpolation")
plot(x, zeros(length(x)), fillrange=ypoly, color="lightgray")
plot!(x,ypoly, color="blue", lw=2)
scatter!(xgiven,ygiven, color="black")
plot!(legend=nothing)
plot!(ylim=[0,13])
p2=plot!(title="polynomial fit")
plot(x, min(yinterp, ypoly), fillrange=max(yinterp, ypoly), color="lightgray")
plot!(x,ypoly, color="blue", lw=2)
plot!(x,yinterp, color="red", lw=2)
scatter!(xgiven,ygiven, color="black")
plot!(legend=nothing)
plot!(ylim=[0,13])
p3=plot!(title="difference")
plot(p1, p2, p3, layout=(3,1), size=(400,900))
x = 0:5
fpts = [2, 1, 2.5, 3.5, 1, 1.75]
Fi = LinearInterpolation(x, fpts)
Fs = CubicSplineInterpolation(x, fpts)
xx = LinRange(0,5,1000)
fi = Fi(xx)
fs = Fs(xx)
plot(xx, zeros(length(xx)), fillrange=fi, color="lightgray", label=nothing)
plot!(xx,fi, color="red", lw=2, label="trapazoid")
plot!(xx,fs, color="black", lw=2, label="function")
plot!(x,fpts, line=:stem, marker=:circle, lw=2, color="black", label=nothing)
plot!(legend_foreground_color=nothing)
plot!(ylim=[0,4.1])
annotate!(x[1]+0.55,0.5, L"I_0")
annotate!(x[2]+0.55,0.5, L"I_1")
annotate!(x[3]+0.55,0.5, L"I_2")
annotate!(x[4]+0.55,0.5, L"I_3")
annotate!(x[5]+0.55,0.5, L"I_4")
annotate!(x[1],fpts[1]+0.35, L"f_0")
annotate!(x[2],fpts[2]+0.35, L"f_1")
annotate!(x[3],fpts[3]+0.35, L"f_2")
annotate!(x[4],fpts[4]+0.35, L"f_3")
annotate!(x[5],fpts[5]+0.35, L"f_4")
annotate!(x[6],fpts[6]+0.35, L"f_5")
n = 7
x = 0:n-1
fpts = [2, 1, 2.5, 3.5, 1, 1.75, 3]
dx = x[2]-x[1]
Fs = CubicSplineInterpolation(x, fpts)
xx = LinRange(0,n-1,1000)
fs = Fs(xx)
nI = Integer((n-1)/2)
p = plot()
for i in 1:nI
a = 1/2/dx^2*(fpts[1+2*(i-1)] -2*fpts[2+2*(i-1)] + fpts[3+2*(i-1)])
b = 1/dx*(-1.5*fpts[1+2*(i-1)] +2*fpts[2+2*(i-1)] - 0.5*fpts[3+2*(i-1)])
c = fpts[1+2*(i-1)]
xi = LinRange(0+2*(i-1),2+2*(i-1),1000)
fi = a*(xi.-xi[1]).^2 .+ b*(xi.-xi[1]) .+ c
plot!(xi,fi, color="red", lw=4, label= i==1 ? "Simpson" : nothing)
plot!(xi,zeros(length(xi)), fillrange=fi, color="lightgray", label=nothing)
end
plot!(xx,fs, color="black", lw=2, label="function")
scatter!(x,fpts, marker=:circle, lw=2, color="black", label=nothing)
plot!(x[1:2:end],fpts[1:2:end], line=:stem, marker=:circle, lw=2, color="black", label=nothing)
plot!(legend_foreground_color=nothing)
plot!(ylim=[0,5])
annotate!(x[1]+1,0.5, L"I_0")
annotate!(x[1]+3,0.5, L"I_1")
annotate!(x[1]+5,0.5, L"I_2")
annotate!(x[1],fpts[1]+0.35, L"f_0")
annotate!(x[2],fpts[2]+0.35, L"f_1")
annotate!(x[3],fpts[3]+0.55, L"f_2")
annotate!(x[4],fpts[4]+0.35, L"f_3")
annotate!(x[5]+0.1,fpts[5]+0.45, L"f_4")
annotate!(x[6],fpts[6]+0.45, L"f_5")
annotate!(x[7],fpts[7]+0.35, L"f_6")
display(p)
Here, we have done $\pm\frac{\Delta x}{3}f_0$ and $\pm\frac{\Delta x}{3}f_{n-1}$ to get 2 as the coefficient of $f_0$ and $f_{n-1}$ so that the even sum can go over all even points; this then results in the final correction term.
Simpson's rule requires an odd number of points: $n$ is odd.
The global integration error is $\propto\Delta x^3\rightarrow \mathcal{O}(\Delta x^3)$.
Could do Simpson's rule for $n-1$ points and Trapazoid rule for the last interval.
Could do Simpson's rule for $n-1$ points, then integrate a parabola through the last three points over the last sub interval.
Could use Simpson's 3/8 rule.
$$\int_{x_0}^{x_3}f(x)dx = \frac{3}{8}\Delta x(f_0 + 3f_1 + 3f_2 + f_3).$$
Hence,
\begin{align} a&=\frac{f_2-f_1}{\Delta x}, \\ b&= 2f_1 - f_2 \end{align}$$I_0 = \int_{x_0}^{x_1}f(x)dx = \frac{a\Delta x^2}{2} + b\Delta x,$$Insert $a$ and $b$ to get:
$$I_0 = \Delta x\left(\frac{3}{2}f_1 - \frac{1}{2}f_2\right).$$