LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= restart; with(Student[NumericalAnalysis]): # We define the function to interpolate f:= x -> x^2; # Nodes Xlist:= [-1, 0, 1,3]; # Number of points( nodes) k:= nops(Xlist): #Compute cubic polynomial on ith segment. ` for i to k - 1 do S[i]:= a[i] + b[i]*(x - Xlist[i]) + c[i]*(x - Xlist[i])^2 + d[i]*(x - Xlist[i])^3; # Compute its first derivative. dS[i]:= diff(S[i], x); # Compute its second derivative. ddS[i]:= diff(dS[i], x); end do: # We must find the 4(k-1) unkown coefficients of the cubics. # We will use a system 4(k-1) equations in these unkown coefficients which #we obtain by constraints on the cubics and their derivatives. for i to k - 1 do # Each cubic must equal to f at left endpoint. EQ[i]:= subs(x = Xlist[i], S[i]) = f(Xlist[i]); end do: for i to k - 1 do # Each cubic must agree with f at right endpoint. EQ[k - 1 + i]:= subs(x = Xlist[i + 1], S[i]) = f(Xlist[i + 1]); end do: for i to k - 2 do #1st derivatives must agree at nodes. EQ[2*k - 2 + i]:= subs(x = Xlist[i + 1], dS[i]) = subs(x = Xlist[i + 1], dS[i + 1]); end do: for i to k - 2 do #2nd derivatives must agree at nodes. EQ[3*k - 4 + i]:= subs(x = Xlist[i + 1], ddS[i]) = subs(x = Xlist[i + 1], ddS[i + 1]); end do: # We now have 4k-6 equations. We need 2 more to have a total of 4k-4, the number of coefficients. # We get these 2 equations from the endpoints. There are two common ways of doing this, called #Clamped Splines and Natural (Free) Splines. # Require 1st derivative to agree with function Clamp:= eval( dS,x = Xlist ) = D(f)(Xlist): # at the end points. Clamp:= eval( dS[k - 1], x = Xlist[k] ) = D(f)(Xlist[k]): # The natural (free) end conditions are: # Require 2nd derivative to be zero at end points Natural:= eval( ddS, x = Xlist ) = 0: Natural:= eval(ddS[k - 1], x = Xlist[k] ) = 0: #Collect all 3k-4 equations. ClampedEqs:= {Clamp, Clamp, seq(EQ[i], i = 1 .. 4*k - 6)}: NaturalEqs:= {Natural, Natural, seq(EQ[i], i = 1 .. 4*k - 6)}; # Solve to obtain the unknown coefficients ClampedCoeffs:= solve ( ClampedEqs, [ seq ( [ a[i], b[i], c[i],d[i]][], i = 1 .. k - 1 ) ] )[]; # HOW CAN I make a general unapply available for any number of splines functions, Here, By hand I try to #modify the code but with three splines there is a problem ClampedSpline:= unapply ( piecewise ( x < eval( Xlist, ClampedCoeffs), expand( eval( S, ClampedCoeffs) ), ( x < eval( Xlist, ClampedCoeffs), expand( eval( S, ClampedCoeffs) ), expand( eval( S, ClampedCoeffs) ) ), x ); plot([f(x), ClampedSpline(x)], x = -1 .. 1); NaturalCoeffs:= solve ( NaturalEqs, [ seq ( [ a[i], b[i], c[i],d[i]][], i = 1 .. k - 1 ) ] )[]: NaturalSpline:= unapply ( piecewise ( x < eval( Xlist, NaturalCoeffs), expand( eval( S, NaturalCoeffs) ), ( x < eval( Xlist, NaturalCoeffs), expand( eval( S, NaturalCoeffs) ), expand( eval( S, NaturalCoeffs) ) ), x ); #Manually edit this line when using a different number of data points. plot([f(x), NaturalSpline(x)], x = -1 .. 1); Zio2I0kieEc2IkYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlKiQ5JCIiI0YlRiVGJQ== NyYhIiIiIiEiIiIiIiQ= PC4vJkkiYUc2IjYjIiIiRigvJkYlNiMiIiMiIiEvJkYlNiMiIiRGKC8sJCZJImNHRiZGJ0YsRi0vLCZGNEYsJkkiZEdGJkYnIiInLCQmRjVGK0YsLywmRjxGLCZGOUYrRjosJCZGNUYwRiwvLCZGQUYsJkY5RjAiIzdGLS8sKCZJImJHRiZGJ0YoRjRGLEY4RjEmRklGKy8sKEZKRihGPEYsRj9GMSZGSUYwLywqRiRGKEZIRihGNEYoRjhGKEYtLywqRipGKEZKRihGPEYoRj9GKEYoLywqRi9GKEZNRixGQSIiJUZEIiIpIiIq Ny4vJkkiYUc2IjYjIiIiRigvJkkiYkdGJkYnISIjLyZJImNHRiZGJ0YoLyZJImRHRiZGJyIiIS8mRiU2IyIiI0YzLyZGK0Y2RjMvJkYvRjZGKC8mRjJGNkYzLyZGJTYjIiIkRigvJkYrRkBGNy8mRi9GQEYoLyZGMkZARjM= Error, `;` unexpected #Compare our spline with the one given by Maple's built in package. with(CurveFitting): Ylist:= map(f, Xlist); xy := [[Xlist, Ylist], [Xlist, Ylist], [Xlist, Ylist], [Xlist, Ylist]; # MapleNatural:= CubicSpline(Xlist, Ylist, x, endpoints = 'natural'); P1:= CubicSpline(xy, independentvar=x, endpoints = 'natural'): P2:=expand(Interpolant(P1)); plot( [P2, NaturalSpline(x)], x = -1 .. 1); # Draw(P1); NyYiIiIiIiFGIyIiKg== Error, `;` unexpected LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=