Inverse quadratic interpolation¶
Here we look for a root of \(x+\cos(10x)\) that is close to 1.
using FundamentalsNumericalComputation
f = x -> x + cos(10*x);
interval = [0.5,1.5];
plot(f,interval...,grid=:y,label="Function",legend=:bottomright,xlabel="x",ylabel="y")
r = nlsolve(x->f(x[1]),[1.]).zero
1-element Array{Float64,1}:
0.9678884021293104
We choose three values to get the iteration started.
x = [0.8,1.2,1]
y = @. f(x)
scatter!(x,y,label="initial points")
If we were using “forward” interpolation, we would ask for the polynomial interpolant of \(y\) as a function of \(x\). But that parabola has no real roots.
using Polynomials
q = fit(x,y,2); # interpolating polynomial
plot!(x->q(x),interval...,l=:dash,label="interpolant",ylim=[-.1,3])
To do inverse interpolation, we swap the roles of \(x\) and \(y\) in the interpolation.
plot(f,interval...,grid=:y,label="Function",legend=:bottomright,xlabel="x",ylabel="y")
scatter!(x,y,label="initial points")
q = fit(y,x,2); # interpolating polynomial
plot!(y->q(y),y->y,-.1,2.6,l=:dash,label="inverse interpolant")
We seek the value of \(x\) that makes \(y\) zero. This means evaluating \(q\) at zero.
@show x = [x; q(0)];
@show y = [y; f(x[end])];
x = [x; q(0)] = [0.8, 1.2, 1.0, 1.103981385440472]
y = [y; f(x[end])] = [0.6544999661913865, 2.043853958732492, 0.16092847092354756, 1.1482065231940184]
We repeat the process a few more times.
for k = 4:8
q = fit(y[k-2:k],x[k-2:k],2)
x = [x; q(0)]
y = [y; f(x[k+1])]
end
Here is the sequence of errors.
err = @. x - r
9-element Array{Float64,1}:
-0.16788840212931033
0.23211159787068958
0.03211159787068962
0.1360929833111617
0.015347343251992385
0.0032683144541431064
0.0004617433331599585
6.295567151926029e-6
3.1584921345029215e-9
The error seems to be superlinear, but subquadratic.
logerr = @. log(abs(err));
ratios = @. logerr[2:end] / logerr[1:end-1]
8-element Array{Float64,1}:
0.8184775454303989
2.3542970942956387
0.5800188686031708
2.094252632650976
1.3702986042756682
1.341928370566436
1.5592295620791041
1.6344120701878158