Given that f(x)= 1/(1+x^2) construct the clamped cubic spline that interpolates the function at x= [-2,2]

Question

It's a numerical analysis question

Answer 1

Divide the x interval [α,β]=[-2,2] into three ordered points (-2,0.2), (1,0.5), (2,0.2). So n=3. Refer to these as (x₀,y₀), (x₁,y₁), (x_2,y₂). The question doesn't specify the number of splines (interpolants), so as an illustration I'll just use 2 (that is, n-1), which are splines between (-2,0.2) and (1,0.5) (knot), and (1,0.5) and (2,0.2). (In general, for n data points there would be n-1 splines and n-2 knots.)

A cubic spline: S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)²+d_i(x-x_i)³ where i=0, 1, or 2 in this example. The splines must pass through the 3 given points. Therefore we have:

S₀(x)=a₀+b₀(x-x₀)+c₀(x-x₀)²+d₀(x-x₀)³=f(x₀)=y₀ for x∈[x₀,x₁];

S₁(x)=a₁+b₁(x-x₁)+c₁(x-x₁)²+d₁(x-x₁)³=f(x₁)=y₁ for x∈[x₁,x₂];

S₂(x)=a₂+b₂(x-x₂)+c₂(x-x₂)²+d₂(x-x₂)³=f(x₂)=y₂ for x≥x₂.

These are the pieces of the piecewise function S(x). We have 12 (=4n) unknown coefficients here, so we need to find more conditions to satisfy the equations. For the required condition of continuity, the end-point of the first interval is also the start-point of the second interval:

S₀(x₁)=S₁(x₁)=f(x₁)=f(1)=y₁=0.5 and, similarly for the end-point of the second interval and (β,f(β)):  S₁(x₂)=S₂(x₂)=f(x₂)=f(2)=y₂=0.2=f(β).

Also: S₀(x₀)=S₀(α)=S₀(-2)=f(-2)=y₀=0.2.

We can now partly solve for these unknowns:

S₀(x)=a₀+b₀(x-x₀)+c₀(x-x₀)²+d₀(x-x₀)³ at (-2,0.2) becomes:

S₀(-2)=a₀=f(-2)=y₀=0.2, so a₀=0.2, because x₀-x₀=0. Also:

S₂(2)=a₂=f(2)=y₂=0.2, so a₂=0.2.

S₁(1)=a₁=f(1)=y₁=0.5, so a₁=0.5.

S₀(1)=S₁(1)=0.5⇒S₀(1)=0.2+3b₀+9c₀+27d₀=0.5.

S₁(2)=S₂(2)=0.2⇒S₁(2)=0.5+b₁+c₂+d₂=0.2.

Now we turn to the "smoothness" of the knot. This involves another condition, that the first derivatives are equal at common points on the splines. In this example there's only one common point at (1,0.5).

S₀'(1)=S₁'(1) and S₁'(2)=S₂'(2).

S₀'(x)=b₀+2c₀(x-x₀)+3d₀(x-x₀)²

S₁'(x)=b₁+2c₁(x-x₁)+3d₁(x-x₁)²

S₂'(x)=b₂+2c₂(x-x₂)+3d₂(x-x₂)²

The second derivative is used to equate the curvature at the knot.

S₀"(1)=S₁"(1) and S₁"(2)=S₂"(2).

S₀"(x₀)=2c₀+6d₀(x-x₀)

S₁"(x₁)=2c₁+6d₁(x-x₁)

S₂"(x₂)=2c₂+6d₂(x-x₂)

Another condition distinguishes natural splines from clamped splines and enables us to find all the coefficients. This condition applies specifically to the end-points of the encompassing interval [α,β].

For clamped cubic splines we need the first derivatives at the end-points (α,f(α)) and (β,f(β)). We can only work out all the coefficients if we know f'(α) and f'(β). f'(x)=-2x/(1+x²)²; f'(-2)=4/25=0.16; f'(2)=-4/25=-0.16, because we know f(x)=1/(1+x²), just as we knew f(x) for x=-2, 1, 2=0.2, 0.5, 0.2. For natural cubic splines we only have to know that the second derivatives are zero at the end-points (α,f(α)) and (β,f(β)).

S₀'(-2)=b₀=f'(-2)=0.16, so b₀=0.16; and S₂'(2)=b₂=f'(2)=-0.16. So we now have all a_i, b₀ and b₂. We need b₁, all c_i and all d_i (7 coefficients). Using all the conditions and making substitutions for known coefficients, we have:

S₀(1)=S₁(1): 0.2+0.48+9c₀+27d₀=0.5;

S₁(2)=S₂(2): 0.5+b₁+c₁+d₁=0.2;

S₀'(1)=S₁'(1): 0.16+2c₀(1+2)+3d₀(1+2)²=0.16+6c₀+27d₀=b₁;

S₁'(2)=S₂'(2): b₁+2c₁+3d₁=-0.16;

S₀"(1)=S₁"(1): 2c₀+6d₀(1+2)=2c₀+18d₀=2c₁;

S₁"(2)=S₂"(2): 2c₁+6d₁=2c₂ (to be omitted).

But we have no use for c₂ or d₂ because we don't need a spline for anything beyond x=β, a₂ was only used in S₁(2)=S₂(2), and b₂ was only used in S₁'(2)=S₂'(2). So the last second derivative equation (containing c₂) can be omitted. That leaves 5 equations and 5 unknown coefficients: b₁, c₀, c₁, d₀, d₁. We have all we need. We can solve using a 5×5 matrix A or elimination and substitution.

Matrix A=

⎛ 0  9   0  27  0 ⎞

⎜ 1   0  1    0 1 ⎟

⎟ 1 -6   0 -27 0 ⎟

⎜ 1  0   2   0  3 ⎟

⎝ 0  2 -2  18  0 ⎠

Matrix X=

⎛ b₁ ⎞

⎜c₀ ⎜

⎜c₁ ⎟

⎜ d₀ ⎜

⎝ d₁ ⎠

and B=

⎛ -0.18 ⎞

⎜ -0.30 ⎟

⎜  0.16 ⎟

⎜-0.16 ⎟

⎝    0    ⎠

We solve for X in AX=B.

d₀=-(0.18+9c₀)/27=-(0.02+c₀)/3,

b₁-6c₀+0.18+9c₀=0.16, b₁=-0.02-3c₀,

c₁=c₀-(0.18+9c₀)/3=-2c₀-0.6,

d₁=-0.3+0.02+3c₀+2c₀+0.6=5c₀+0.32,

-0.02-3c₀-4c₀-1.2+15c₀+0.96=-0.16, 8c₀=-0.16+0.02+1.2-0.96=0.1, c₀=0.0125.

Therefore: d₀=-13/1200=-0.01083 approx; b₁=-0.0575; c₁=-0.625; d₁=0.3825.

S₀(x)=0.2+0.16(x+2)+0.0125(x+2)²-(13/1200)(x+2)³,

S₁(x)=0.5-0.0575(x-1)-0.625(x-1)²+0.3825(x-1)³.

Comments

I don't have access to MATLAB or other computational apps. The example I used is just an illustration of the method. In practice there would be multiple data points and, without an appropriate app, computations would be very laborious and time-consuming, as well as being prone to error. There may be errors in my solution despite the relative simplicity of the example.