'**************************************************** '* Program to demonstrate the general integration * '* subroutine. Example is the integral of SIN(X) * '* from X1 to X2 in double precision. * '* ------------------------------------------------ * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* '* * '* ------------------------------------------------ * '* SAMPLE RUN: * '* * '* Integration of SIN(X) * '* * '* Input start point end point: 0, 0.5 * '* Integral from 0.000 to 0.500 equals 0.1224231 * '* Error check: 1 * '* * '**************************************************** defint i-n defdbl a-h,o-z iv=14 'Number of pooints in table dim X(iv),Y(iv),XM(iv+3),Z(iv) cls print print " Integration of SIN(X)" print 'Input table for i=1 to iv read X(i), Y(i) next i '------------------------------------------------------------ 'Sine table values from Handbook of mathematical 'functions by M. Abramowitz and I.A. Stegun, 'NBS, june 1964 data 0.000#,0.00000000#,0.125#,0.12467473# data 0.217#,0.21530095#,0.299#,0.29456472#,0.376#,0.36720285# data 0.450#,0.43496553#,0.520#,0.49688014#,0.589#,0.55552980# data 0.656#,0.60995199#,0.721#,0.66013615# data 0.7853981634#,0.7071067812# data 0.849#,0.75062005#,0.911#,0.79011709# data 0.972#,0.82601466# '------------------------------------------------------------ Input " Input start point, end point : ",x1, x2 print gosub 2000 print using " Integral from ##.### to ##.### equals ##.#######"; x1; x2; zz print using " Error check: #"; z1 print end '******************************************************** '* Akima spline fitting subroutine * '* ---------------------------------------------------- * '* The input table is X(i), Y(i), where Y(i) is the * '* dependant variable. The interpolation point is xx, * '* which is assumed to be in the interval of the table * '* with at least one table value to the left, and three * '* to the right. The interpolated returned value is yy. * '* n is returned as an error check (n=0 implies error). * '* It is also assumed that the X(i) are in ascending * '* order. * '******************************************************** 1000 'Check to see if interpolation point is correct n=1 if xx>=X(1) then goto 1100 n=0 : return 1100 if xx<=X(iv-3) then goto 1200 n=0 : return 1200 X(0)=2#*X(1)-X(2) 'Calculate Akima coefficients for i=1 to iv-1 'Shift i to i+2 XM(i+2)=(Y(i+1)-Y(i))/(X(i+1)-X(i)) next i XM(iv+2)=2#*XM(iv+1)-XM(iv) XM(iv+3)=2#*XM(iv+2)-XM(iv+1) XM(2)=2#*Xm(3)-XM(4) XM(1)=2#*XM(2)-XM(3) for i=1 to iv a=ABS(XM(i+3)-XM(i+2)) b=ABS(XM(i+1)-XM(i)) if a+b<>0 then goto 1300 Z(i)=(XM(i+2)+XM(i+1))/2# goto 1350 1300 Z(i)=(a*XM(i+1)+b*XM(i+2))/(a+b) 1350 next i 'Find relevant table interval i=0 1400 i=i+1 if xx>X(i) then goto 1400 i=i-1 'Begin interpolation b=X(i+1)-X(i) a=xx-X(i) yy=Y(i)+Z(i)*a+(3#*XM(i+2)-2#*Z(i)-Z(i+1))*a*a/b yy=yy+(Z(i)+Z(i+1)-2#*XM(i+2))*a*a*a/(b*b) return '******************************************************** '* General integration subroutine (ITEG) * '* Interpolation by Akima (or other). Integration by * '* enhanced trapezoidal rule with Richardson extrapola- * '* tion to give cubic accuracy. The integration range * '* is (x1,x2). It is assumed that x1error). * '******************************************************** 2000 zz=0 : z1=0 'Check to see if end points are correct if x1X(iv-3) then return 'If x1>x2 then switch and set flag if x1X(i+1) then goto 3100 'If not, integral is simple n1=n1+1 : d=yy : xx=x2 'Find end point yy value gosub 1000 xi1=(yy+d)*(x2-x1)/2# return 'At least one table interval must be summed over 3100 j1=i xi1=xi1+(yy+Y(i+1))*(X(i+1)-xx)/2# 'Any more intervals? If not, finish integral with end point 3150 if x2X(i+1) then goto 3600 xx=x1+(x2-x1)/2# gosub 1000 xi2=xi2+(d+yy)*(x2-x1)/4# return 3600 xx=x1+(X(i+1)-x1)/2# j1=i gosub 1000 xi2=xi2+(yy+d)*(xx-x1)/2# xi2=xi2+(yy+T(j1+1))*(X(j1+1)-xx)/2# 3650 if x2