#restart; read `W:/hw3.txt`;

Help := proc();

print(`NewtonDD, HermiteInter, FreeCubic`);
print(`ClampedCubic, TwoPoint`);
end:

# Section 3.2 
# Newton's Divided Difference Method
# Program from the book page 121.

# Input: the list of point x , the list  
# of point y and variable x;
#
# Output: print the table and 
# return Interpolating polynomial 
# of degree nops(x)-1

# Example: Table 3.8 in the book
#
# Try: 
# X := [1.0,1.3,1.6,1.9,2.2]:
# Y := [0.7651977,0.6200860,
#  0.4554022,0.2818186,0.1103623]:
# NewtonDD(X,Y,x); 

# Another example: problem 1a)
# Check the answer in the back to the book
#
# Try: 
# X := [8.1,8.3,8.6,8.7]:
# Y := [16.94410,17.56492,
#  18.50515,18.82091]:
# NewtonDD(X,Y,x); 

NewtonDD := proc(X,Y,x) local i, j, n, F,P;

n := nops(X);
F := [[0$n]$n];

for i from 1 to n do
	F[i][1] := Y[i];
od:

for i from 2 to n do
	
	for j from 2 to i do

F[i][j] := (F[i][j-1]-F[i-1][j-1])/(X[i]-X[i-j+1]);

	od:
od:

#Write interpolating polynomial

P := 0; 
for i from 1 to n do

  P := P + F[i][i]*mul(x-X[k],k=1..i-1);

od:

for i from 1 to nops(F) do
	print(F[i]);
od:

P;
end:


# Section 3.3 
# Hermite polynomial
# Program from the book page 135.

# Input: the list of point x , the list  
# of point y , the list of y' and variable x;
#
# Output: print the table and 
# return Interpolating polynomial 
# of degree 2*nops(x)-1

# Example: Table 3.14 in the book
#
# Try: 
# X := [1.3,1.6,1.9]:
# Y := [0.6200860,0.4554022,0.2818186]:
# Yprime := [-0.5220232,-0.5698959,-0.5811571]:
# HermiteInter(X,Y,Yprime,x); 

# Another example: problem 1a)
# Check the answer in the back to the book
#
# Try: 
# X := [8.3,8.6]:
# Y := [17.56492,18.50515]:
# Yprime := [3.116256,3.151762]:
# HermiteInter(X,Y,Yprime,x); 


HermiteInter := proc(X,Y,Yprime,x) local n, i, j, z, Q, A;

n := nops(X);	
z := [0$2*n];
Q := [[0$(2*n)]$(2*n)];

for i from 1 to n do

	z[2*i-1] := X[i];	
	z[2*i] := X[i];
	Q[2*i-1][1] := Y[i];
	Q[2*i][1] := Y[i];
	Q[2*i][2] := Yprime[i];

	if i <> 1 then
		Q[2*i-1][2] := 
(Q[2*i-1][1]-Q[2*i-2][1])/(z[2*i-1]-z[2*i-2]); 
	fi:

od:

for i from 3 to 2*n do
	for j from 3 to i do
		Q[i][j] := (Q[i][j-1]-Q[i-1][j-1])/(z[i]-z[i-j+1]);
	od:
od:

for i from 1 to nops(Q) do
	print(Q[i]);
od:

A := [Q[1][1], seq( Q[i][i]*mul(x-X[floor(j/2)],j=2..i), i = 2..2*n)];

add(A[i], i =1..2*n);

end:

## Section 3.4 
## Free Cubic Spline

# Input: the list of point x , the list  
# of point y , the list of y' and variable x;
#
# Output: Cublic Spline Interpolation 

# Example: problem 3a)
#
# Try: 
# X := [8.3,8.6]:
# Y := [17.56492,18.50515]:
# FreeCubic(X,Y,x); 

# Another example: problem 3c)
# Check the answer in the back to the book
#
# Try: 
# X := [-0.5,-0.25,0]:
# Y := [-0.02475,0.3349375,1.101]:
# FreeCubic(X,Y,x); 

FreeCubic := proc(X,Y,x) 

local j, n, P, dP, ddP, var, eq, a, b, c, d, i;

n := nops(X);

#set up equations
P   := {seq(a[i]+b[i]*(x-X[i])+
c[i]*(x-X[i])^2+d[i]*(x-X[i])^3, i =1..n-1)};
dP  := diff(P,x);
ddP := diff(diff(P,x),x);

#all the conditions
eq :={  seq(subs(x=X[j],P[j])=Y[j],j=1..n-1),
		seq(subs(x=X[j+1],P[j])=Y[j+1],j=1..n-1),
seq(subs(x=X[j+1],dP[j])=subs(x=X[j+1],dP[j+1]),j=1..n-2),
seq(subs(x=X[j+1],ddP[j])=subs(x=X[j+1],ddP[j+1]),j=1..n-2),
subs(x=X[1],ddP[1])=0, subs(x=X[n],ddP[n-1])=0};

#variables
var := {seq(op([a[j],b[j],c[j],d[j]]),j=1..n-1)};

#solve for the solution
var:=solve(eq,var):

return(
{seq( S[j-1](x) = subs(var,P[j]), j= 1..n-1)});

end:


## Section 3.4 
## Clamped Cubic Spline

# Input: the list of point x , the list  
# of point y , the list of y' and variable x;
#
# Output: Clamped Spline Interpolation 

# Example: problem 7a) (there are typos in the book)
#
# Try: 
# X := [8.3,8.6]:
# Y := [17.56492,18.50515]:
# YPrime := [3.116256,3.151762]:
# ClampedCubic(X, Y, YPrime, x); 

# Another example: problem 7c)
# Check the answer in the back to the book
#
# Try: 
# X := [-0.5,-0.25,0]:
# Y := [-0.02475,0.3349375,1.101]:
# YPrime := [0.751,4.002]:
# ClampedCubic(X, Y, YPrime, x); 

ClampedCubic := proc(X,Y,YPrime,x) 

local j, n, P, dP, ddP, var, eq, a, b, c, d, i;

n := nops(X);

#set up equations
P   := {seq(a[i]+b[i]*(x-X[i])+
c[i]*(x-X[i])^2+d[i]*(x-X[i])^3, i =1..n-1)};
dP  := diff(P,x);
ddP := diff(diff(P,x),x);

#all the conditions
eq :={  seq(subs(x=X[j],P[j])=Y[j],j=1..n-1),
		seq(subs(x=X[j+1],P[j])=Y[j+1],j=1..n-1),
seq(subs(x=X[j+1],dP[j])=subs(x=X[j+1],dP[j+1]),j=1..n-2),
seq(subs(x=X[j+1],ddP[j])=subs(x=X[j+1],ddP[j+1]),j=1..n-2),
subs(x=X[1],dP[1])=YPrime[1], subs(x=X[n],dP[n-1])=YPrime[2]};

#variables
var := {seq(op([a[j],b[j],c[j],d[j]]),j=1..n-1)};

#solve for the solution
var:=solve(eq,var):

return(
{seq( S[j-1](x) = subs(var,P[j]), j= 1..n-1)});

end:

# Section 4.1
# Two-point formula

# Input: [x_0,x_1], [y_0,y_1];
#
# Output: the estimation of y'_0 using
# Two-point formula

# Example: row 1 of table 4.1
#
# Try: 
# X := [1.8,1.9]:
# Y := [ln(1.8),0.64185389]:
# TwoPoint(X,Y); 

# Another example: problem 1a)
# Check the answer in the back to the book
#
# Try: 
# X := [0.5,0.6]:
# Y := [0.4794,0.5646]:
# TwoPoint(X,Y); 

TwoPoint := proc(X,Y)

	return((Y[2]-Y[1])/(X[2]-X[1]));

end: