"File F-ADD-UTIL-VER4.MTH, copyright (c) 1997-2000 by R&D Publishing."

"See add-summary.mth for a summary of new functions defined below:"

"Substitute x=a into the expression u."

SUBST(u,x,a):=LIM(u,x,a)

"The secant line of u(x) through x = a and x = a + h."

SECANT(u,x,a,h):=y=(SUBST(u,x,a+h)-SUBST(u,x,a))/h*(x-a)+SUBST(u,x,a)

"The tangent line of u(x) at x = a."

TANGENT(u,x,a):=y=SUBST(u,x,a)+SUBST(DIF(u,x),x,a)*(x-a)

"Helper functions for CURVEFIT."

POLY(x,a,n):=SUM(a SUB (i+1)*x^i,i,0,n)

DPOLY(x,a,n):=SUM(i*a SUB (i+1)*x^(i-1),i,1,n)

EQNS(data,ddata,a,n):=APPEND(VECTOR(POLY(data SUB i SUB 1,a,n)=data SUB i SUB~
 2,i,1,DIMENSION(data),1),VECTOR(DPOLY(ddata SUB i SUB 1,a,n)=ddata SUB i SUB~
 2,i,1,DIMENSION(ddata),1))

UNK(n):=RHS(VECTOR(SOLVE(upsilon=upsilon,upsilon),i,1,n,1)` SUB 1)

ANS(x,data,ddata,a,n):=IF(DIMENSION(ans_:=SOLVE(EQNS(data,ddata,a,n),[a]))=0,~
[],POLY(x,(RHS(ans_)) SUB 1,n))

"Finds the polynomial of the right degree through the nx2-data matrix."

CURVEFIT1(x,data):=ANS(x,data,[],UNK(DIMENSION(data)),DIMENSION(data)-1)

CURVEFIT2(x,data,ddata):=ANS(x,data,ddata,UNK(DIMENSION(data)+DIMENSION(ddata~
)),DIMENSION(data)+DIMENSION(ddata)-1)

CURVEFIT(x,data,ddata):=IF(DIMENSION(ddata)>0,CURVEFIT2(x,data,ddata),CURVEFI~
T1(x,data),CURVEFIT1(x,data))

"Quadratic spline function interpolating data points with initial slope m1."

SPLINE_AUX(x,data,m):=SUM(CURVEFIT(x,[data SUB k,data SUB (k+1)],[[data SUB k~
 SUB 1,m SUB k SUB 2]])*CHI(data SUB k SUB 1,x,data SUB (k+1) SUB 1),k,1,DIME~
NSION(data)-1)

SLOPE(data,m1):=ITERATES([v SUB 1+1,2*(data SUB (v SUB 1+1) SUB 2-data SUB (v~
 SUB 1) SUB 2)/(data SUB (v SUB 1+1) SUB 1-data SUB (v SUB 1) SUB 1)-v SUB 2]~
,v,[1,m1],DIMENSION(data)-1)

"Note that SLOPE returns the matrix [[1,m1], [2,m2], ...]."

SPLINE(x,data,m1):=SPLINE_AUX(x,data,SLOPE(data,m1))

"Newton algorithm"

NEWT_ITERATES(u,x,x0,k):=ITERATES(x-u/DIF(u,x),x,x0,k)

NEWT(u,x,x0,k):=IF(k>0,NEWT_ITERATES(u,x,x0,k),?,SUBST(x-u/DIF(u,x),x,x0))

"This produces a vector which when plotted demonstrates Newton's method."

DRAW_NEWT(u,x,x0,k):=VECTOR([[v,0],[v,SUBST(u,x,v)],[NEWT(u,x,v),0]],v,ITERAT~
ES(NEWT(u,x,w),w,x0,k))

DRAW_COMPLEX(v):=VECTOR([RE(z),IM(z)],z,v)

"Bisection method helper"

BIS_AUX(u,x,a,b):=IF(SUBST(u,x,a)*SUBST(u,x,(a+b)/2)<0,[a,(a+b)/2],[(a+b)/2,b~
])

"Bisection method"

BISECT(u,x,v0,k):=IF(k>0,ITERATES(BIS_AUX(u,x,v SUB 1,v SUB 2),v,v0,k),?,BIS_~
AUX(u,x,v0 SUB 1,v0 SUB 2))

"Formula for the left-endpoint method for integrating u(x) over [a,b] with n ~
subdivisions."

LEFT(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+k*(b-a)/n),k,0,n-1)

"Formula for the midpoint method for integrating u(x) over [a,b] with n subdi~
visions."

MID(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+(k+0.5)*(b-a)/n),k,0,n-1)

"Formula for the right-endpoint method for integrating u(x) over [a,b] with n~
 subdivisions."

RIGHT(u,x,n,a,b):=(b-a)/n*SUM(SUBST(u,x,a+k*(b-a)/n),k,1,n)

"Formula for the trapezoid method for integrating u(x) over [a,b] with n subd~
ivisions."

TRAP(u,x,n,a,b):=(b-a)/(2*n)*(SUBST(u,x,a)+SUBST(u,x,b)+2*SUM(SUBST(u,x,a+k*(~
b-a)/n),k,1,n-1))

"Formula for Simpson's rule for integrating u(x) over [a,b] with n subdivisio~
ns."

SIMP(u,x,n,a,b):=(b-a)/(6*n)*(SUBST(u,x,a)+SUBST(u,x,b)+2*SUM(SUBST(u,x,a+k*(~
b-a)/n),k,1,n-1)+4*SUM(SUBST(u,x,a+(k+1/2)*(b-a)/n),k,0,n-1))

"The box and trapezoid drawing functions used in the graphical demonstrations~
 of numerical integration techniques."

D_BOX(x1,y1,x2,y2):=[[x1,y1],[x2,y1],[x2,y2],[x1,y2],[x1,y1]]

D_TRAP(x1,y1,x2,y2,x3,y3,x4,y4):=[[x1,y1],[x2,y2],[x3,y3],[x4,y4],[x1,y1]]

DRAW_LEFT(u,x,n,a,b):=VECTOR(D_BOX(t,0,t+(b-a)/n,SUBST(u,x,t)),t,a,b-(b-a)/n,~
(b-a)/n)

DRAW_RIGHT(u,x,n,a,b):=VECTOR(D_BOX(t,0,t+(b-a)/n,SUBST(u,x,t+(b-a)/n)),t,a,b~
-(b-a)/n,(b-a)/n)

DRAW_TRAP(u,x,n,a,b):=VECTOR(D_TRAP(t,0,t+(b-a)/n,0,t+(b-a)/n,SUBST(u,x,t+(b-~
a)/n),t,SUBST(u,x,t)),t,a,b-(b-a)/n,(b-a)/n)

"The solution to the differential equation (DE) y'+p(x)y=q(x) with y(x0)=y0."

DE(p,q,x,y,x0,y0):=y=(y0+INT(q*#e^INT(p,x,x0,x),x,x0,x))/#e^INT(p,x,x0,x)

"Direction field helper function."

SEG(rc,x,y,ëx,ëy):=IF(ABS(rc)>1 AND ëy<ABS(rc)*ëx,[[x-ëy/rc,y-ëy],[x+ëy/rc,y+~
ëy]],[[x-ëx,y-rc*ëx],[x+ëx,y+rc*ëx]])

"The direction field (DF) for the differential equation: y'=r(x,y) with a  de~
termined by x0<x<xm with m subdivisions and y0<y<ym with n subdivisions."

DF(r,x,x0,xm,m,y,y0,yn,n):=VECTOR(VECTOR(SEG(LIM(r,[x,y],[x0+j*(xm-x0)/m,y0+k~
*(yn-y0)/n]),x0+j*(xm-x0)/m,y0+k*(yn-y0)/n,(xm-x0)/(4*m),(yn-y0)/(4*n)),j,0,m~
),k,0,n)

"The EULER function gives an approximate solution to: y'=r(x,y) with y(x0)=y0~
."

EULER(r,x,y,x0,y0,xn,n):=ITERATES(v+(xn-x0)/n*[1,LIM(r,[x,y],v)],v,[x0,y0],n)