(* AME 521 Numerical Methods--J. M. Powers, Fall 1999
	
	This program calculates the cubic splines for the example given in AME 521 lecture.    *)

(* First define the discrete points.  Here we take N=3, so
			we have 4 points.  We also choose them to be equally spaced
			in x ϵ [0,1] *)

x0 =0.;
x1 = 1./3.;
x2 = 2./3.;
x3 = 1.;

(* h is the constant spacing between points *)

h = x3-x2;

(* λ is the weighting function for the cantilever end conditions; 0 < λ < 1  *)

λ = 0;

(* Define the function to be interpolated.  Here we choose a simple exponential. It could be anything, including raw data. *)

f[x_]:=Sin[2*x]

(* Define the equations and variables we need to solve to get the second derivatives at each interior point.  Here we have already incorporated the cantilevered end condition. *)

eqns = {
f0pp==λ*f1pp,
       f0pp + 4*f1pp + f2pp==(6/h^2)(f[x0]-2 *f[x1]+f[x2]),
       f1pp +4*f2pp+f3pp==(6/h^2)(f[x1]-2*f[x2]+f[x3]),
       f3pp==λ*f2pp
};
vars={f0pp,f1pp,f2pp,f3pp};

(*  Solve the equations. *)

{soln}=Solve[eqns,vars];

(* Knowing the second derivatives, we can now get each interpolating function. *)

f0 =(x1-x)^3 f0pp/(6h) +
			(x-x0)^3 f1pp /(6h)+
			(x1-x)(f[x0]/h - (h/6)f0pp)+
			(x-x0)(f[x1]/h-(h/6)f1pp)/.soln;
			
f1 =(x2-x)^3 f1pp/(6h) +
			(x-x1)^3 f2pp /(6h)+
			(x2-x)(f[x1]/h - (h/6)f1pp)+
			(x-x1)(f[x2]/h-(h/6)f2pp)/.soln;

f2 =(x3-x)^3 f2pp/(6h) +
			 (x-x2)^3 f3pp /(6h)+
			(x3-x)(f[x2]/h - (h/6)f2pp)+
			(x-x2)(f[x3]/h-(h/6)f3pp)/.soln;

(* Combine each interpolating function to form a global function. *)


ff = If[x0<x<x1,f0,If[x1<x<x2,f1,If[x2<x<x3,f2,0]]];

(* Plot the spline function, the original function and then  the error in the spline approximation *)


Plot[{ff,f[x]},{x,0,1},AxesLabel->{"x","y"}]
Plot[ff-f[x],{x,0,1},AxesLabel->{"x","Error: y(x)-f(x)"}]