"File ADD-UTIL.MTH, copyright (c) 1997-1999 by R&D Publishing." "See add-summary.mth for a summary of new functions defined below:" "Substitute x=a into the expression u." "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,DIM(data),1),VECTOR(DPOLY(ddata SUB i SUB 1,a,n)=ddata SUB i SUB 2,i,1~ ,DIM(ddata),1)) UNK(n):=VECTOR((SOLUTIONS(upsilon=upsilon,upsilon)) SUB 1,i,1,n,1) ANS(x,data,ddata,a,n,v):=PROG(v:=SOLUTIONS(EQNS(data,ddata,a,n),a),IF(DIM(v)=~ 0,[],POLY(x,v SUB 1,n))) "Finds the polynomial of the right degree through the nx2-data matrix." CURVEFIT1(x,data):=ANS(x,data,[],UNK(DIM(data)),DIM(data)-1) CURVEFIT2(x,data,ddata):=ANS(x,data,ddata,UNK(DIM(data)+DIM(ddata)),DIM(data)~ +DIM(ddata)-1) CURVEFIT(x,data,ddata):=IF(DIM(ddata)>0,CURVEFIT2(x,data,ddata),CURVEFIT1(x,d~ ata),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) SLOPE1(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],DIM(data)-1) "Note that SLOPE1 returns the matrix [[1,m1], [2,m2], ...]." SPLINE(x,data,m1):=SPLINE_AUX(x,data,SLOPE1(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