restart: Digits := 20: # CUBIC SPLINE RAYLEIGH-RITZ ALGORITHM 11.6 # # To approximate the solution to the boundary-value problem # # -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1, Y(0)=Y(1)=0 # # with a sum of cubic splines: # # INPUT: Integer n # # OUTPUT: Coefficients C(0),...,C(n+1) of the basis functions INTE := proc(J,JJ) local inte: inte := JJ-J+3: RETURN(inte): end: XINT := proc(XU,XL,A1,B1,C1,D1,A2,B2,C2,D2,A3,B3,C3,D3) local AA,BB,CC,DD,EE,FF,GG,XHIGH,XLOW,I1,xint,C: AA := A1*A2: BB := A1*B2+A2*B1: CC := A1*C2+B1*B2+C1*A2: DD := A1*D2+B1*C2+C1*B2+D1*A2: EE := B1*D2+C1*C2+D1*B2: FF := C1*D2+D1*C2: GG := D1*D2: C[9] := AA*A3: C[8] := (AA*B3+BB*A3)/2: C[7] := (AA*C3+BB*B3+CC*A3)/3: C[6] := (AA*D3+BB*C3+CC*B3+DD*A3)/4: C[5] := (BB*D3+CC*C3+DD*B3+EE*A3)/5: C[4] := (CC*D3+DD*C3+EE*B3+FF*A3)/6: C[3] := (DD*D3+EE*C3+FF*B3+GG*A3)/7: C[2] := (EE*D3+FF*C3+GG*B3)/8: C[1] := (FF*D3+GG*C3)/9: C[0] := (GG*D3)/10: XHIGH := 0: XLOW := 0: for I1 from 1 to 10 do XHIGH := (XHIGH+C[I1-1])*XU: XLOW := (XLOW+C[I1-1])*XL: od: xint := XHIGH-XLOW: RETURN(xint): end: print(`This is the Cubic Spline Rayleigh-Ritz Method.`): OK := FALSE: print(`Input functions P(X), Q(X), and F(X) in terms of x separated by a space.`): print(`For example: `): print(`1 3.141592654^2 2*3.141592654^2*sin(3.1415926548*x)`): print(`separated by at least one space.`): Pf := scanf(`%a`)[1]: Qf := scanf(`%a`)[1]: Ff := scanf(`%a`)[1]: print(`P(x) = `):print(Pf):print(`Q(x) = `):print(Qf):print(`R(x) = `):print(Ff): FPL := evalf(subs(x=0,diff(Ff,x))): FPR := evalf(subs(x=1,diff(Ff,x))): QPL := evalf(subs(x=0,diff(Qf,x))): QPR := evalf(subs(x=1,diff(Qf,x))): PPL := evalf(subs(x=0,diff(Pf,x))): PPR := evalf(subs(x=1,diff(Pf,x))): F := unapply(Ff,x,y,z): Q := unapply(Qf,x,y,z): P := unapply(Pf,x,y,z): while OK = FALSE do print(`Input a positive integer n, where x(0) = 0.0 and x(n+1) = 1.0`): N := scanf(`%d`)[1]: print(`n = `):print(N): if N <= 0 then print(`Number must be a positive integer.`): else OK := TRUE: fi: od: if OK = TRUE then print(`Choice of output method:`): print(`1. Output to screen`): print(`2. Output to text File`): print(`Please enter 1 or 2.`): FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG): if FLAG = 2 then print(`Input the file name in the form - drive:\134\134name.ext`): print(`for example A:\134\134OUTPUT.DTA`): NAME := scanf(`%s`)[1]: print(`Output file is `):print(NAME): OUP := fopen(NAME,WRITE,TEXT): else OUP := default: fi: fprintf(OUP, `CUBIC SPLINE RAYLEIGH-RITZ METHOD\134n\134n`): # Step 1 H := 1/(N+1): N1 := N+1: N2 := N+2: N3 := N+3: # Initialize matrix A at 0, note that A[I,N+3] = B[I] for I3 from 1 to N2 do for J from 1 to N3 do A[I3-1,J-1] := 0: od: od: # Step 2 # X[1] = 0, ... , X[I] = (I-1)*H, ... , X[N+1] = 1-H, X[N+2] = 1 for I3 from 1 to N2 do X[I3-1] := (I3-1)*H: od: # STEPS 3 and 4 are implemented in what follows. Initialize coefficients # CO[I,J,K], DCO[I,J,K] for I3 from 1 to N2 do for J from 1 to 4 do # JJ corresponds the coefficients of phi and phi' to the proper interval # involving J JJ := I3+J-3: CO[I3-1,J-1,0] := 0: CO[I3-1,J-1,1] := 0: CO[I3-1,J-1,2] := 0: CO[I3-1,J-1,3] := 0: E := I3-1: OK := TRUE: if JJ < I3-2 or JJ >= I3+2 then OK := FALSE: fi: if I3 = 1 and JJ < I3 then OK := FALSE: fi: if I3 = 2 and JJ < I3-1 then OK := FALSE: fi: if I3 = N+1 and JJ > N+1 then OK := FALSE: fi: if I3 = N+2 and JJ >= N+2 then OK := FALSE: fi: if OK = TRUE then if JJ <= I3-2 then CO[I3-1,J-1,0] := (((-E+6)*E-12)*E+8)/24: CO[I3-1,J-1,1] := ((E-4)*E+4)/(8*H): CO[I3-1,J-1,2] := (-E+2)/(8*H^2): CO[I3-1,J-1,3] := 1/(24*H^3): OK := FALSE: else if JJ > I3 then CO[I3-1,J-1,0] := (((E+6)*E+12)*E+8)/24: CO[I3-1,J-1,1] := ((-E-4)*E-4)/(8*H): CO[I3-1,J-1,2] := (E+2)/(8*H^2): CO[I3-1,J-1,3] := -1/(24*H^3): OK := FALSE: else if JJ > I3-1 then CO[I3-1,J-1,0] := ((-3*E-6)*E*E+4)/24: CO[I3-1,J-1,1] := (3*E+4)*E/(8*H): CO[I3-1,J-1,2] := (-3*E-2)/(8*H^2): CO[I3-1,J-1,3] := 1/(8*H^3): if I3 <> 1 and I3 <> N+1 then OK := FALSE: fi: else CO[I3-1,J-1,0] := ((3*E-6)*E*E+4)/24: CO[I3-1,J-1,1] := (-3*E+4)*E/(8*H): CO[I3-1,J-1,2] := (3*E-2)/(8*H^2): CO[I3-1,J-1,3] := -1/(8*H^3): if I3 <> 2 and I3 <> N+2 then OK := FALSE: fi: fi: fi: fi: fi: if OK = TRUE then if I3 <= 2 then AA := 1/24: BB := -1/(8*H): CC := 1/(8*H^2): DD := -1/(24*H^3): if I3 = 2 then CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-AA: CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-BB: CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-CC: CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-DD: else CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-4*AA: CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-4*BB: CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-4*CC: CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-4*DD: fi: else EE := N+2: AA := (((-EE+6)*EE-12)*EE+8)/24: BB := ((EE-4)*EE+4)/(8*H): CC := (-EE+2)/(8*H^2): DD := 1/(24*H^3): if I3 = N+1 then CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-AA: CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-BB: CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-CC: CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-DD: else CO[I3-1,J-1,0] := CO[I3-1,J-1,0]-4*AA: CO[I3-1,J-1,1] := CO[I3-1,J-1,1]-4*BB: CO[I3-1,J-1,2] := CO[I3-1,J-1,2]-4*CC: CO[I3-1,J-1,3] := CO[I3-1,J-1,3]-4*DD: fi: fi: fi: DCO[I3-1,J-1,0] := 0: DCO[I3-1,J-1,1] := 0: DCO[I3-1,J-1,2] := 0: E := I3-1: OK := TRUE: if JJ < I3-2 or JJ >= I3+2 then OK := FALSE: fi: if I3 = 1 and JJ < I3 then OK := FALSE: fi: if I3 = 2 and JJ < I3-1 then OK := FALSE: fi: if I3 = N+1 and JJ > N+1 then OK := FALSE: fi: if I3 = N+2 and JJ >= N+2 then OK := FALSE: fi: if OK = TRUE then if JJ <= I3-2 then DCO[I3-1,J-1,0] := ((E-4)*E+4)/(8*H): DCO[I3-1,J-1,1] := (-E+2)/(4*H^2): DCO[I3-1,J-1,2] := 1/(8*H^3): OK := FALSE: else if JJ > I3 then DCO[I3-1,J-1,0] := ((-E-4)*E-4)/(8*H): DCO[I3-1,J-1,1] := (E+2)/(4*H^2): DCO[I3-1,J-1,2] := -1/(8*H^3): OK := FALSE: else if JJ > I3-1 then DCO[I3-1,J-1,0] := (3*E+4)*E/(8*H): DCO[I3-1,J-1,1] := (-3.0*E-2.0)/(4.0*H^2): DCO[I3-1,J-1,2] := 3/(8*H^3): if I3 <> 1 and I3 <> N+1 then OK := FALSE: fi: else DCO[I3-1,J-1,0] := (-3*E+4)*E/(8*H): DCO[I3-1,J-1,1] := (3*E-2)/(4*H^2): DCO[I3-1,J-1,2] := -3/(8*H^3): if I3 <> 2 and I3 <> N+2 then OK := FALSE: fi: fi: fi: fi: fi: if OK = TRUE then if I3 <= 2 then AA := -1/(8*H): BB := 1/(4*H^2): CC := -1/(8*H^3): if I3 = 2 then DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-AA: DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-BB: DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-CC: else DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-4*AA: DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-4*BB: DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-4*CC: fi: else EE := N+2: AA := ((EE-4)*EE+4)/(8*H): BB := (-EE+2)/(4*H^2): CC := 1/(8*H^3): if I3 = N+1 then DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-AA: DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-BB: DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-CC: else DCO[I3-1,J-1,0] := DCO[I3-1,J-1,0]-4*AA: DCO[I3-1,J-1,1] := DCO[I3-1,J-1,1]-4*BB: DCO[I3-1,J-1,2] := DCO[I3-1,J-1,2]-4*CC: fi: fi: fi: od: od: # Output the basis functions. fprintf(OUP, `Basis Function: A + B*X + C*X**2 + D*X**3\134n\134n`): fprintf(OUP, ` A B C D\134n\134n`): for I3 from 1 to N2 do fprintf(OUP, `phi( %d ) \134n\134n`, I3): for J from 1 to 4 do if I3 <> 1 or (J <> 1 and J <> 2) then if I3 <> 2 or J <> 2 then if I3 <> N1 or J <> 4 then if I3 <> N2 or (J <> 3 and J <> 4) then JJ1 := I3+J-3: JJ2 := I3+J-2: fprintf(OUP, `On (X( %d ), X( %d )) `, JJ1, JJ2): for K from 1 to 4 do fprintf(OUP, ` %12.8f `, CO[I3-1,J-1,K-1]): od: fprintf(OUP, `\134n\134n`): fi: fi: fi: fi: od: od: # Obtain coefficients for F, P, Q for I3 from 1 to N2 do AA[I3-1] := evalf(F(X[I3-1])): od: XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*FPL: XA[N2-1] := 3.0*FPR-3.0*(AA[N2-1]-AA[N2-2])/H: XL[0] := 2.0*H: XU[0] := 0.5: XZ[0] := XA[0]/XL[0]: for I3 from 2 to N1 do XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H: XL[I3-1] := H*(4.0-XU[I3-2]): XU[I3-1] := H/XL[I3-1]: XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]: od: XL[N2-1] := H*(2.0-XU[N2-2]): XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]: CC[N2-1] := XZ[N2-1]: for I3 from 1 to N1 do J := N2-I3: CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]: BB[J-1] := (AA[J]-AA[J-1])/H-H*(CC[J]+2.0*CC[J-1])/3.0: DD[J-1] := (CC[J]-CC[J-1])/(3.0*H): od: for I3 from 1 to N1 do AF[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]: BF[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]: CF[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]: DF[I3-1] := DD[I3-1]: od: for I3 from 1 to N2 do AA[I3-1] := evalf(P(X[I3-1])): od: XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*PPL: XA[N2-1] := 3.0*PPR-3.0*(AA[N2-1]-AA[N2-2])/H: XL[0] := 2.0*H: XU[0] := 0.5: XZ[0] := XA[0]/XL[0]: for I3 from 2 to N1 do XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H: XL[I3-1] := H*(4.0-XU[I3-2]): XU[I3-1] := H/XL[I3-1]: XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]: od: XL[N2-1] := H*(2.0-XU[N2-2]): XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]: CC[N2-1] := XZ[N2-1]: for I3 from 1 to N1 do J := N2-I3: CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]: BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0: DD[J-1] := (CC[J]-CC[J-1])/(3.0*H): od: for I3 from 1 to N1 do AP[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]: BP[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]: CP[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]: DP[I3-1] := DD[I3-1]: od: for I3 from 1 to N2 do AA[I3-1] := evalf(Q(X[I3-1])): od: XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*QPL: XA[N2-1] := 3.0*QPR-3.0*(AA[N2-1]-AA[N2-2])/H: XL[0] := 2.0*H: XU[0] := 0.5: XZ[0] := XA[0]/XL[0]: for I3 from 2 to N1 do XA[I3-1] := 3.0*(AA[I3]-2.0*AA[I3-1]+AA[I3-2])/H: XL[I3-1] := H*(4.0-XU[I3-2]): XU[I3-1] := H/XL[I3-1]: XZ[I3-1] := (XA[I3-1]-H*XZ[I3-2])/XL[I3-1]: od: XL[N2-1] := H*(2.0-XU[N2-2]): XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1]: CC[N2-1] := XZ[N2-1]: for I3 from 1 to N1 do J := N2-I3: CC[J-1] := XZ[J-1]-XU[J-1]*CC[J]: BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0: DD[J-1] := (CC[J]-CC[J-1])/(3.0*H): od: for I3 from 1 to N1 do AQ[I3-1] := ((-DD[I3-1]*X[I3-1]+CC[I3-1])*X[I3-1]-BB[I3-1])*X[I3-1]+AA[I3-1]: BQ[I3-1] := (3.0*DD[I3-1]*X[I3-1]-2.0*CC[I3-1])*X[I3-1]+BB[I3-1]: CQ[I3-1] := CC[I3-1]-3.0*DD[I3-1]*X[I3-1]: DQ[I3-1] := DD[I3-1]: od: # Steps 5-9 are implemented in what follows for I3 from 1 to N2 do # indices for limits of integration for A[I,J] and B[I] J1 := min(I3+2,N+2): J0 := max(I3-2,1): J2 := J1-1: # Integrate over each subinterval where phi(I) is nonzero for JJ from J0 to J2 do # limits of integration for each call XU := X[JJ]: XL := X[JJ-1]: # coefficients of bases K := INTE(I3,JJ): A1 := DCO[I3-1,K-1,0]: B1 := DCO[I3-1,K-1,1]: C1 := DCO[I3-1,K-1,2]: D1 := 0: A2 := CO[I3-1,K-1,0]: B2 := CO[I3-1,K-1,1]: C2 := CO[I3-1,K-1,2]: D2 := CO[I3-1,K-1,3]: # call subprogram for integrations A[I3-1,I3-1] := A[I3-1,I3-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ-1],A1,B1,C1,D1,A1,B1,C1,D1)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ-1],A2,B2,C2,D2,A2,B2,C2,D2): A[I3-1,N+2]:=A[I3-1,N+2]+XINT(XU,XL,AF[JJ-1],BF[JJ-1],CF[JJ-1],DF[JJ-1],A2,B2,C2,D2,1,0,0,0): od: # compute A[I,J] for J = I+1, ..., Min(I+3,N+2) K3 := I3+1: if K3 <= N2 then K2 := min(I3+3,N+2): for J from K3 to K2 do J0 := max(J-2,1): for JJ from J0 to J2 do XU := X[JJ]: XL := X[JJ-1]: K := INTE(I3,JJ): A1 := DCO[I3-1,K-1,0]: B1 := DCO[I3-1,K-1,1]: C1 := DCO[I3-1,K-1,2]: D1 := 0: A2 := CO[I3-1,K-1,0]: B2 := CO[I3-1,K-1,1]: C2 := CO[I3-1,K-1,2]: D2 := CO[I3-1,K-1,3]: K := INTE(J,JJ): A3 := DCO[J-1,K-1,0]: B3 := DCO[J-1,K-1,1]: C3 := DCO[J-1,K-1,2]: D3 := 0: A4 := CO[J-1,K-1,0]: B4 := CO[J-1,K-1,1]: C4 := CO[J-1,K-1,2]: D4 := CO[J-1,K-1,3]: A[I3-1,J-1] := A[I3-1,J-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ- 1],A1,B1,C1,D1,A3,B3,C3,D3)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ- 1],A2,B2,C2,D2,A4,B4,C4,D4): od: A[J-1,I3-1] := A[I3-1,J-1]: od: fi: od: # Step 10 for I3 from 1 to N1 do for J from I3+1 to N2 do CC := A[J-1,I3-1]/A[I3-1,I3-1]: for K from I3+1 to N3 do A[J-1,K-1] := A[J-1,K-1]-CC*A[I3-1,K-1]: od: A[J-1,I3-1] := 0: od: od: C[N2-1] := A[N2-1,N3-1]/A[N2-1,N2-1]: for I3 from 1 to N1 do J := N1-I3+1: C[J-1] := A[J-1,N3-1]: for KK from J+1 to N2 do C[J-1] := C[J-1]-A[J-1,KK-1]*C[KK-1]: od: C[J-1] := C[J-1]/A[J-1,J-1]: od: # Step 14 # Output coefficients fprintf(OUP, `\134nCoefficients: c(1), c(2), ... , c(n+1)\134n\134n`): for I3 from 1 to N1 do fprintf(OUP, ` %12.6e \134n`, C[I3-1]): od: fprintf(OUP, `\134n`): # compute and output value of the approximation at the nodes fprintf(OUP, `The approximation evaluated at the nodes:\134n\134n`): fprintf(OUP, ` Node Value\134n\134n`): for I3 from 1 to N2 do S := 0: for J from 1 to N2 do J0 := max(J-2,1): J1 := min(J+2,N+2): SS := 0: if I3 < J0 or I3 >= J1 then S := S + C[J-1]*SS: else K := INTE(J,I3): SS := ((CO[J-1,K-1,3]*X[I3-1]+CO[J-1,K-1,2])*X[I3-1]+CO[J-1,K-1,1])*X[I3-1]+CO[J-1, K-1,0]: S := S + C[J-1]*SS: fi: od: fprintf(OUP, `%12.8f %12.8f\134n`, X[I3-1], S): od: fi: if OUP <> default then fclose(OUP): print(`Output file `,NAME,` created successfully`): fi: SU9UaGlzfmlzfnRoZX5DdWJpY35TcGxpbmV+UmF5bGVpZ2gtUml0en5NZXRob2QuRzYi SWRvSW5wdXR+ZnVuY3Rpb25zflAoWCksflEoWCksfmFuZH5GKFgpfmlufnRlcm1zfm9mfnh+fCtzZXBhcmF0ZWR+Ynl+YX5zcGFjZS5HNiI= SS5Gb3J+ZXhhbXBsZTp+RzYi SVQxfjMuMTQxNTkyNjU0XjJ+MiozLjE0MTU5MjY1NF4yKnNpbigzLjE0MTU5MjY1NDgqeClHNiI= SUFzZXBhcmF0ZWR+Ynl+YXR+bGVhc3R+b25lfnNwYWNlLkc2Ig== SShQKHgpfj1+RzYi IiIi SShRKHgpfj1+RzYi JCI0O1B3bU9TLydwKSohIz0= SShSKHgpfj1+RzYi LCQtSSRzaW5HNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHNiI2IywkSSJ4R0YoJCIrYUVmVEohIiokIjVLdV9MdCEpMyNSKD4hIz0= SWhuSW5wdXR+YX5wb3NpdGl2ZX5pbnRlZ2Vyfm4sfndoZXJlfngoMCl+PX4wLjB+YW5kfngobisxKX49fjEuMEc2Ig== SSVufj1+RzYi IiIq STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI= STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi STcyLn5PdXRwdXR+dG9+dGV4dH5GaWxlRzYi STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig== SSpJbnB1dH5pc35HNiI= IiIi CUBIC SPLINE RAYLEIGH-RITZ METHOD Basis Function: A + B*X + C*X**2 + D*X**3 A B C D phi( 1 ) On (X( 1 ), X( 2 )) 0.00000000 5.00000000 -75.00000000 291.66666667 On (X( 2 ), X( 3 )) 0.33333333 -5.00000000 25.00000000 -41.66666667 phi( 2 ) On (X( 0 ), X( 1 )) 0.00000000 0.00000000 0.00000000 0.00000000 On (X( 2 ), X( 3 )) -0.20833333 8.75000000 -62.50000000 125.00000000 On (X( 3 ), X( 4 )) 1.12500000 -11.25000000 37.50000000 -41.66666667 phi( 3 ) On (X( 1 ), X( 2 )) 0.00000000 0.00000000 0.00000000 41.66666667 On (X( 2 ), X( 3 )) 0.16666667 -5.00000000 50.00000000 -125.00000000 On (X( 3 ), X( 4 )) -1.83333333 25.00000000 -100.00000000 125.00000000 On (X( 4 ), X( 5 )) 2.66666667 -20.00000000 50.00000000 -41.66666667 phi( 4 ) On (X( 2 ), X( 3 )) -0.04166667 1.25000000 -12.50000000 41.66666667 On (X( 3 ), X( 4 )) 1.29166667 -18.75000000 87.50000000 -125.00000000 On (X( 4 ), X( 5 )) -5.45833333 48.75000000 -137.50000000 125.00000000 On (X( 5 ), X( 6 )) 5.20833333 -31.25000000 62.50000000 -41.66666667 phi( 5 ) On (X( 3 ), X( 4 )) -0.33333333 5.00000000 -25.00000000 41.66666667 On (X( 4 ), X( 5 )) 4.16666667 -40.00000000 125.00000000 -125.00000000 On (X( 5 ), X( 6 )) -11.83333333 80.00000000 -175.00000000 125.00000000 On (X( 6 ), X( 7 )) 9.00000000 -45.00000000 75.00000000 -41.66666667 phi( 6 ) On (X( 4 ), X( 5 )) -1.12500000 11.25000000 -37.50000000 41.66666667 On (X( 5 ), X( 6 )) 9.54166667 -68.75000000 162.50000000 -125.00000000 On (X( 6 ), X( 7 )) -21.70833333 118.75000000 -212.50000000 125.00000000 On (X( 7 ), X( 8 )) 14.29166667 -61.25000000 87.50000000 -41.66666667 phi( 7 ) On (X( 5 ), X( 6 )) -2.66666667 20.00000000 -50.00000000 41.66666667 On (X( 6 ), X( 7 )) 18.16666667 -105.00000000 200.00000000 -125.00000000 On (X( 7 ), X( 8 )) -35.83333333 165.00000000 -250.00000000 125.00000000 On (X( 8 ), X( 9 )) 21.33333333 -80.00000000 100.00000000 -41.66666667 phi( 8 ) On (X( 6 ), X( 7 )) -5.20833333 31.25000000 -62.50000000 41.66666667 On (X( 7 ), X( 8 )) 30.79166667 -148.75000000 237.50000000 -125.00000000 On (X( 8 ), X( 9 )) -54.95833333 218.75000000 -287.50000000 125.00000000 On (X( 9 ), X( 10 )) 30.37500000 -101.25000000 112.50000000 -41.66666667 phi( 9 ) On (X( 7 ), X( 8 )) -9.00000000 45.00000000 -75.00000000 41.66666667 On (X( 8 ), X( 9 )) 48.16666667 -200.00000000 275.00000000 -125.00000000 On (X( 9 ), X( 10 )) -79.83333333 280.00000000 -325.00000000 125.00000000 On (X( 10 ), X( 11 )) 41.66666667 -125.00000000 125.00000000 -41.66666667 phi( 10 ) On (X( 8 ), X( 9 )) -14.29166667 61.25000000 -87.50000000 41.66666667 On (X( 9 ), X( 10 )) 71.04166667 -258.75000000 312.50000000 -125.00000000 On (X( 10 ), X( 11 )) -80.83333333 247.50000000 -250.00000000 83.33333333 phi( 11 ) On (X( 9 ), X( 10 )) -21.33333333 80.00000000 -100.00000000 41.66666667 On (X( 10 ), X( 11 )) 221.66666667 -730.00000000 800.00000000 -291.66666667 Coefficients: c(1), c(2), ... , c(n+1) 3.057803e-05 1.256556e+00 2.390140e+00 3.289736e+00 3.867321e+00 4.066339e+00 3.867321e+00 3.289736e+00 2.390140e+00 1.256556e+00 The approximation evaluated at the nodes: Node Value 0.00000000 0.00000000 0.10000000 0.30901648 0.20000000 0.58778555 0.30000000 0.80901689 0.40000000 0.95105659 0.50000000 0.99999997 0.60000000 0.95105659 0.70000000 0.80901689 0.80000000 0.58778555 0.90000000 0.30901648 1.00000000 0.00000000 JSFH