(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; paletteColors = 128; currentKernel; ] :[font = title; inactive; Cclosed; preserveAspect; startGroup] Approximation of Functions-part iv :[font = subtitle; inactive; preserveAspect] An interpretation based upon Theodore J. Rivilin's "An Introduction to the Approximation of Functions" :[font = subsubtitle; inactive; preserveAspect; endGroup] M.A. Lachance :[font = section; inactive; Cclosed; preserveAspect; startGroup] Polynomial interpolation :[font = text; inactive; Cclosed; preserveAspect; startGroup] Linear system approach :[font = input; preserveAspect; startGroup] left = -1; right = 1; n=4; x = Table[left + Random[](right-left),{n+1}]; y = Table[2 Random[]-1,{n+1}]; pow[0,t_] := 1; pow[i_,t_] := t^i; mat = Table[Table[x[[i1]]^j1,{j1,0,n}],{i1,Length[x]}]; dat = Table[{y[[i1]]},{i1,n+1}]; unk = Table[{a[i1]},{i1,0,n}]; {MatrixForm[mat],MatrixForm[unk],MatrixForm[dat]} :[font = input; preserveAspect; endGroup] unk = Flatten[Inverse[mat].dat]; p[t_] := Table[t^j1,{j1,0,n}].unk Plot[p[t],{t,left,right}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}] :[font = text; inactive; preserveAspect; endGroup] Exercise: Redefine x and y so that the data to be fitted have uniformly spaced points in [-1,1], and y=Abs[x-.5]. Construct the third degree polynomial interpolate of this data. How does it compare to either the best uniform approximate of the same degree? to the best least square approximate of the same degree? Exercise: Repeat the above but with x having an arcsine distribution in [-1,1], and y=Abs[x-.5]. Construct the third degree polynomial interpolate of this data. How does it compare to either the best uniform approximate of the same degree? to the best least square approximate of the same degree? ;[s] 18:0,1;10,0;19,2;21,0;25,2;27,0;62,3;71,0;129,3;134,0;318,1;328,0;354,2;356,0;366,3;373,0;431,3;436,0;621,-1; 4:9,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0;3,13,10,Courier,1,12,0,0,0;4,13,9,Times,1,12,0,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Vandermonde matrices and determinants :[font = text; inactive; preserveAspect] From a practical standpoint, the approach above is not a very good one for constructing an interpolating polynomial. From a theoretical perspective, however, the coeffcient matrix above is called a Vandermonde matrix, and its determinant has a neat closed form: ;[s] 5:0,0;237,1;255,0;265,1;276,0;314,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; startGroup] x = Table[left + Random[](right-left),{n+1}]; mat = Table[Table[x[[i1]]^j1,{j1,0,n}],{i1,Length[x]}]; Det[mat] Product[Product[(x[[i1]]-x[[j1]]),{i1,1,j1-1}],{j1,n+1}] :[font = output; output; inactive; preserveAspect] -0.0003712141801348916 ;[o] -0.000371214 :[font = output; output; inactive; preserveAspect; endGroup] -0.0003712141801348915 ;[o] -0.000371214 :[font = input; preserveAspect; startGroup] Clear[x] MatrixForm[mat = Table[Table[x[i1]^j1,{j1,0,n}],{i1,0,n}]] Simplify[Det[mat]] :[font = output; output; inactive; preserveAspect] MatrixForm[{{1, x[0], x[0]^2, x[0]^3, x[0]^4}, {1, x[1], x[1]^2, x[1]^3, x[1]^4}, {1, x[2], x[2]^2, x[2]^3, x[2]^4}, {1, x[3], x[3]^2, x[3]^3, x[3]^4}, {1, x[4], x[4]^2, x[4]^3, x[4]^4}}] ;[o] 2 3 4 1 x[0] x[0] x[0] x[0] 2 3 4 1 x[1] x[1] x[1] x[1] 2 3 4 1 x[2] x[2] x[2] x[2] 2 3 4 1 x[3] x[3] x[3] x[3] 2 3 4 1 x[4] x[4] x[4] x[4] :[font = output; output; inactive; preserveAspect; endGroup; endGroup] (-x[0] + x[1])*(x[0] - x[2])*(x[1] - x[2])*(x[0] - x[3])*(x[1] - x[3])* (-x[2] + x[3])*(x[0] - x[4])*(x[1] - x[4])*(x[2] - x[4])*(x[3] - x[4]) ;[o] (-x[0] + x[1]) (x[0] - x[2]) (x[1] - x[2]) (x[0] - x[3]) (x[1] - x[3]) (-x[2] + x[3]) (x[0] - x[4]) (x[1] - x[4]) (x[2] - x[4]) (x[3] - x[4]) :[font = text; inactive; Cclosed; preserveAspect; startGroup] Newton polynomials and divided differences :[font = text; inactive; preserveAspect; startGroup] Newton basis functions: :[font = input; preserveAspect; endGroup] Clear[pow]; pow[0,t_] := 1; pow[i_,t_] := t^i; left = -1; right = 1; n=4; f[t_] := Exp[t/4.] x={1,2,3,4,5}; y=Table[f[x[[i1]]],{i1,Length[x]}]; NewtonBase[0,t_] := 1 NewtonBase[i_,t_] := (t-x[[i]]) NewtonBase[i-1,t] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Divided difference table: :[font = input; preserveAspect; endGroup] DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; Print["The x values: ",x] Print[" "]; Print["The divided difference table: ",MatrixForm[DDtable]]; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Interpolating polynomials: :[font = input; preserveAspect; startGroup] Clear[p]; p[0,t]=y[[1]]; p[i_,t_]:= Sum[DDtable[[i1,1]] NewtonBase[i1-1,t],{i1,i+1}] MatrixForm[Table[p[i1,t],{i1,0,n}]] MatrixForm[Table[Expand[p[i1,t]],{i1,0,n}]] :[font = text; inactive; preserveAspect; endGroup; endGroup] Exercise: Change f[t] above to t^3. What is p[3,t]? p[4,t]? What is going on? Explain. Exercise: Express p[3,t] in the Horner form. ;[s] 6:0,1;10,0;92,1;102,0;124,2;136,0;138,-1; 3:3,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0;1,13,9,Times,1,12,0,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Random examples: :[font = input; preserveAspect; startGroup] x=Table[2 Random[]-1,{i1,n+1}]; y=Table[pow[3,x[[i1]]],{i1,Length[x]}]; DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; Plot[{p[0,t],p[1,t],p[2,t],p[3,t]},{t,Min[x],Max[x]}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}, Axes->False, PlotStyle->{RGBColor[.5,0,0], RGBColor[.5,.5,0], RGBColor[.5,.5,.5], RGBColor[.5,0,1]}] :[font = text; inactive; preserveAspect; endGroup; endGroup; endGroup] Exercise: Choose x so that the points come from the interval [-eps,eps], where eps is, say, .1, and y=f[x] with f[t]=Exp[t]. Compare p[0,t], p[1,t], p[2,t], and p[3,t] with the firsts four Taylor polynomials. t[0,t] = f[0] t[1,t] = f[0]+f'[0] x t[2,t] = f[0]+f'[0] x+1/2*f''[0] x^2 t[3,t] = f[0]+f'[0] x+1/2*f''[0] x^2 +1/6*f'''[0]x^3 If eps shrinks to .01, how is the comparison affected? ;[s] 6:0,1;10,0;17,2;18,0;100,2;107,0;394,-1; 3:3,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0;2,13,10,Courier,1,12,0,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Runge's function :[font = text; inactive; Cclosed; preserveAspect; startGroup] Construct Newton interpolate on odd number of uniformly spaced points :[font = input; preserveAspect; endGroup] n=16; left = -5; right = 5; Runge[x_] := (1+x^2)^-1; x=Table[left+i1/n(right-left),{i1,0,n}]; y=Table[Runge[x[[i1]]],{i1,Length[x]}]; DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; Clear[p]; p[0,t]=y[[1]]; p[i_,t_]:= Sum[DDtable[[i1,1]] NewtonBase[i1-1,t],{i1,i+1}] Plot[{Runge[t],p[n,t]},{t,left,right}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}, Axes->False] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Construct Newton interpolate on odd number of scaled Chebyshev zeros :[font = input; preserveAspect; startGroup] n=16; left = -5; right = 5; Runge[x_] := (1+x^2)^-1; x=5 Table[Cos[i1/n N[Pi]],{i1,0,n}]; y=Table[Runge[x[[i1]]],{i1,Length[x]}]; DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; Clear[p]; p[0,t]=y[[1]]; p[i_,t_]:= Sum[DDtable[[i1,1]] NewtonBase[i1-1,t],{i1,i+1}] Plot[{Runge[t],p[n,t]},{t,left,right}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}, Axes->False] :[font = text; inactive; preserveAspect; endGroup; endGroup; endGroup] Exercise:Let d denote the smallest distance between the n zeros of the Chebyshev polynomials. Change the x values above so that there is a random perturbation of each, not to exceed d/4. How does the polynomial fit of Runge's function respond to this perturbation? ;[s] 6:0,1;8,0;13,2;14,0;182,2;185,0;267,-1; 3:3,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0;2,13,9,Times,1,12,0,0,0; :[font = text; inactive; preserveAspect; startGroup] Error in interpolation :[font = text; inactive; Cclosed; preserveAspect; startGroup] Construct approximate as above :[font = input; preserveAspect] Clear[pow]; pow[0,t_] := 1; pow[i_,t_] := t^i; left = -1; right = 1; n=4; f[t_] := Sin[2.t] x=Sort[Table[Random[],{n+1}]](n+1) y=Table[f[x[[i1]]],{i1,Length[x]}]; pic= Plot[f[z],{z,x[[1]],x[[n+1]]}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],f[x[[i1]]]}],{i1,n+1}]}, PlotLabel->"function and points to be interpolated"] NewtonBase[0,t_] := 1 NewtonBase[i_,t_] := (t-x[[i]]) NewtonBase[i-1,t] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Divided difference table: :[font = input; preserveAspect; endGroup] DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Interpolating polynomials: :[font = input; preserveAspect; endGroup] Clear[p]; p[0,t]=y[[1]]; p[i_,t_]:= Sum[DDtable[[i1,1]] NewtonBase[i1-1,t],{i1,i+1}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Display result :[font = input; preserveAspect; endGroup; endGroup] MatrixForm[Table[Expand[p[i1,t]],{i1,0,n}]] :[font = text; inactive; preserveAspect; startGroup] Pick an approxiamte to analyze by choosing the degree m <=n ;[s] 3:0,0;57,1;66,0;71,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; startGroup] Clear[phi,t]; m=4; w[t_]:= Product[(t-x[[i1]]),{i1,m+1}] w[t] Plot[{w[z],f[z]-p[m,z]},{z,x[[1]],x[[4]]}, PlotStyle->{Dashing[{.02,.02}], Dashing[{1,0}]}] ;[s] 3:0,0;14,1;19,0;161,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Pick any old point xx in [x[[1]],x[[m+1]]] ;[s] 3:0,0;31,1;33,0;55,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] x xx = x[[1]]+ Random[](x[[m+1]]-x[[1]]) lambda = (f[xx]-p[m,xx])/w[xx]; (* define phi[t] which has m+2 roots *) phi[t_] := Expand[f[t]-p[m,t]-lambda w[t]]; phi[t] Plot[phi[z],{z,x[[1]],x[[m+1]]}, Prolog->{PointSize[.02], Append[ Table[Point[{x[[i1]],0}],{i1,m+1}],Point[{xx,0}]]}, PlotLabel->"phi[t] = f[t]-p[m,t]-lambda*w[t]"] ;[s] 3:0,0;101,1;111,0;343,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Which means phi's (m+1)st derivative has one root, xi ;[s] 4:0,0;52,1;56,0;62,1;65,-1; 2:2,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; startGroup] Clear[Dphi,Df,z]; Dphi[z_]:= Expand[D[phi[t],{t,m+1}]/.{t->z}] Df[z_]:= Expand[D[f[t]-p[m,t],{t,m+1}]/.{t->z}] Dphi[z] Df[z]-lambda Factorial[m+1] xi = FindRoot[Dphi[z]==0,{z,(x[[1]]+x[[m+1]])/2.}][[1,2]] :[font = output; output; inactive; preserveAspect] 21.22701818512149 + 32.*Cos[2.*z] ;[o] 21.227 + 32. Cos[2. z] :[font = output; output; inactive; preserveAspect] 21.22701818512149 + 32.*Cos[2.*z] ;[o] 21.227 + 32. Cos[2. z] :[font = output; output; inactive; preserveAspect; endGroup] 1.148037710232872 ;[o] 1.14804 :[font = input; wordwrap; preserveAspect; endGroup] Plot[Dphi[z],{z,x[[1]],x[[m+1]]}, Prolog->{PointSize[.02],Point[{xi,0}]}, PlotLabel->"(m+1)st derivative of phi[t]"] Plot[{Df[z], lambda Factorial[m+1], -lambda Factorial[m+1]}, {z,x[[1]],x[[m+1]]}, Prolog->{PointSize[.02],Point[{xi,Df[xi]}]}, PlotStyle->{Dashing[{.02,.02}], Dashing[{1,0}],Dashing[{1,0}]}, PlotLabel->"(m+1)st derivative of f[xi]=+/- lambda (m+1)!"] :[font = text; inactive; Cclosed; preserveAspect; startGroup] phi[t] and its (m+1)st derivative, AND the envelope ;[s] 3:0,0;54,1;62,0;63,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] Clear[funcNorm]; funcNorm[f_, p_, a_, b_] := If[p == Infinity, Max[Table[Abs[f[xxx]],{xxx, a,b,(b-a).001}]], NIntegrate[Abs[f[xxx]]^p, {xxx, a, b}]^(1/p)] M = funcNorm[Df,Infinity,x[[1]],x[[m+1]]] Plot[{f[z]-p[m,z], M/Factorial[m+1] w[z], -M/Factorial[m+1] w[z]}, {z,x[[1]],x[[m+1]]}, PlotStyle->{Dashing[{.02,.02}], Dashing[{1,0}],Dashing[{1,0}]}, PlotLabel->"+/- M/(n+1)! w[x]"] :[font = text; inactive; preserveAspect; endGroup; endGroup] and just what is the best envelope function? :[font = text; inactive; preserveAspect; endGroup; endGroup] Exercise: Choose as your polynomial q[t], the best unifrom quadratic approximation to the function f[t]=t^4. i) Define w[t] for your error function f[t]-q[t]. ii) Determine the value of xi so that f'''[xi]=(f[.5]-q[.5])/w[.5] *3! iii) Determine the max M of |f'''[x]| on [-1,1]. iv) Plot the error function f[t]-q[t] superimposed on the two curves M*w[t] and -M*w[t] ;[s] 2:0,1;8,0;374,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = section; inactive; Cclosed; preserveAspect; startGroup] C2 Spline interpolation :[font = text; inactive; Cclosed; preserveAspect; startGroup] Polynominals tend to oscillate: :[font = text; inactive; Cclosed; preserveAspect; startGroup] Construct Newton interpolate on odd number of scaled Chebyshev zeros :[font = input; preserveAspect] NewtonBase[0,t_] := 1 NewtonBase[i_,t_] := (t-x[[i]]) NewtonBase[i-1,t] :[font = input; preserveAspect; endGroup] n=8; left = -5; right = 5; Runge[x_] := (1+x^2)^-1; x=5 Table[Cos[i1/n N[Pi]],{i1,0,n}]; y=Table[Runge[x[[i1]]],{i1,Length[x]}]; DDtable = Table[Table[0,{i1,n+2-j1}],{j1,n+1}]; DDtable[[1]] = y; For[i=2,i<=n+1,i++, For[j=1,j<=n-i+2,j++, DDtable[[i,j]]=(DDtable[[i-1,j+1]]-DDtable[[i-1,j]])/ (x[[j+i-1]]-x[[j]]); ]; ]; Clear[p]; p[0,t]=y[[1]]; p[i_,t_]:= Sum[DDtable[[i1,1]] NewtonBase[i1-1,t],{i1,i+1}] Plot[{Runge[t],p[n,t]},{t,left,right}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}, Axes->False] :[font = text; inactive; preserveAspect; startGroup] Curvature as a measurement of bending energy :[font = text; inactive; Cclosed; preserveAspect; startGroup] Define osculating circle so that its Taylor expansion agrees with f to define curvature ;[s] 5:0,0;24,1;41,0;95,1;104,0;105,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect] Off[General::"spell1", Plot::"plnr"]; Clear[f, x, h, k, r, tan, nor, y]; Series[f[x], {x, 0, 3}] tan[x_] := f[0] + f'[0] x; nor[x_] := f[0] - x/f'[0]; tan[x] nor[x] h = -f'[0]/Sqrt[1 + f'[0]^2] r; k = f[0]+r/Sqrt[1 + f'[0]^2]; y[x_] := -Sqrt[r^2 - (x - h)^2] + k; temp = Solve[Coefficient[Series[y[x], {x, 0, 3}] - Series[f[x], {x, 0, 3}], x^2] == 0, r]; Print["r = ",r = Simplify[temp[[2,1]][[2]]]] :[font = input; preserveAspect; endGroup] f[x_] := (x + .25)^2; tan[x_] := f[0] + f'[0] x; nor[x_] := f[0] - x/f'[0]; h = -f'[0]/Sqrt[1 + f'[0]^2] r; k = f[0]+r/Sqrt[1 + f'[0]^2]; y[x_] := -Sqrt[r^2 - (x - h)^2] + k; Show[ Plot[{0, f[x], tan[x], nor[x]}, {x, -1.25, .75}, Prolog -> {PointSize[.03], Point[{h, k}], Point[{0, f[0]}]}, PlotRange -> {{-1.25, .75}, {-.5, 1.75}}, AspectRatio -> Automatic, Axes -> False], Graphics[Circle[{h,k},r]], Graphics[Text["(x-h)^2+(y-k)^2=r^2", {.25, 2r+.1}]], Graphics[Text["y=f(x)", {.85, 1.6r}]], Graphics[Text["(h,k)", {-.5, r}]], Graphics[Text["r", {-.2, r/2}]], Graphics[Text["y=f[0]+f'[0] x", {-.45, -5r/8}]], Graphics[Text["y=f[0]-x/f'[0]", {.6, -r/4}]], PlotRange -> {{-1.5, 1}, {-.5, 1.75}}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Define "energy" of a curve :[font = input; preserveAspect; endGroup; endGroup; endGroup] Clear[kf,f]; f[t_] := Expand[p[n,t]] kf[t_] := f''[t]/(1+f'[t]^2)^(3/2); Plot[{Runge[t],f[t],kf[t]}, {t,-5,5}, PlotStyle->{ Dashing[{1,0}], Dashing[{1,0}], Dashing[{.025,.05}]}] Plot[{kf[t],kf[t]^2}, {t,-5,5}, PlotStyle->{ Dashing[{1,0}], Dashing[{1,0}], Dashing[{.025,.05}]}] Print["Curvature energy = ",NIntegrate[kf[t]^2,{t,-5,5}]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Construction of C2 splines (the theory): :[font = text; inactive; Cclosed; preserveAspect; startGroup] The development of the interior constraints... :[font = text; inactive; Cclosed; preserveAspect; startGroup] Consider two-segment, piecewise linear splines p''[t] and q''[t] (b - t) za (-a + t) zb (c - t) zb (-b + t) zc p''[t] = ---------- + ----------- and q''[t]= ---------- + ----------- , integrate twice, adding -a + b -a + b -b + c -b + c linear terms to each... za, zb, and zc are unspecified as yet. :[font = input; preserveAspect; endGroup] Off[General::spell]; Clear[p,q,t,a,b,c,d,c1,c2,d1,d2]; p[t_] := (b-t)^3/(b-a)/6 za +(t-a)^3/(b-a)/6 zb + (b-t)/(b-a) c1 + (t-a)/(b-a) c2 q[t_] := (c-t)^3/(c-b)/6 zb +(t-b)^3/(c-b)/6 zc + (c-t)/(c-b) d1 + (t-b)/(c-b) d2 {p''[a],p''[b],q''[b],q''[c]} ;[s] 3:0,0;4,1;18,0;253,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Solve for {c1,c2,d1,d2} by forcing p[a]=ya, p[b]=yb=q[b], and q[c]=yc :[font = input; preserveAspect; endGroup] Clear[c1,c2,d1,d2]; c1=Apart[Solve[p[a]==ya,c1][[1,1]][[2]]] c2=Apart[Solve[p[b]==yb,c2][[1,1]][[2]]] d1=Apart[Solve[q[b]==yb,d1][[1,1]][[2]]] d2=Apart[Solve[q[c]==yc,d2][[1,1]][[2]]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Replace (b-a) with h[[k]], (c-b) with h[[k+1]], {a,b,c} with {x[[k-1]],x[[k]],x[[k+1]]}, {ya,yb,yc} with {y[[k-1]],y[[k]],y[[k+1]]} {za,zb,zc} with {z[[k-1]],z[[k]],z[[k+1]]} to get row entries of matrix :[font = input; preserveAspect; endGroup] skt=p[t]/.{(b-a)->hk}; Print["The k-th spline at t:"] Print[skt = skt/.{b->xkp1,a->xk,za->zk,zb->zkp1,ya->yk,yb->ykp1}]; Print["- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"] Print[" "]; Print["The derivative of 1-st spline at x[[1]]:"] Print[" "]; Lend=p[t]/.{(b-a)->h1}; Lend=D[Lend,t]/.{t->x1,a->x1,b->x2,za->z1,zb->z2} Lend=Lend/.{x2-x1->h1,ya->y1,yb->y2} Print[Simplify[Expand[Lend]]] Print["- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"] Print[" "]; Print["The derivative of k-th spline at x[[k+1]]:"] Print[" "]; LHS=D[skt,t]/.{t->xkp1}; Print[LHS=Simplify[LHS/.{(xkp1-xk)->hk,(xkp2-xkp1)->hkp1}]] Print["- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"] Print[" "]; Print["The derivative of (k+1)-st spline at x[[k+1]]:"] Print[" "]; RHS = q[t]/.{(c-b)->hkp1}; RHS=D[RHS,t]/.{t->xkp1}; RHS=RHS/.{b->xkp1,c->xkp2,zc->zkp2,zb->zkp1,yc->ykp2,yb->ykp1}; Print[RHS=Simplify[RHS/.{(xkp1-xk)->hk,(xkp2-xkp1)->hkp1}]]; Print["- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"] Print[" "]; Print["The derivative of last spline at x[[n]]:"] Print[" "]; Rend=q[t]/.{(c-b)->hnm1}; Rend=D[Rend,t]/.{t->xn,c->xn,b->xnm1,zc->zn,zb->znm1} Rend=Rend/.{xn-xnm1->hnm1,yc->yn,yb->ynm1} Print[Simplify[Expand[Rend]]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Equating the two sides we get an equation involving zk,zkp1,zkp2 :[font = input; preserveAspect; endGroup; endGroup] temp = Simplify[LHS-RHS]; ck = Coefficient[temp,zk] ckp1 = Coefficient[temp,zkp1] ckp2 = Coefficient[temp,zkp2] leftOver=Simplify[temp-ck zk - ckp1 zkp1 - ckp2 zkp2] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Set up and solve tridiagonal system for Rivlin's zero slope end condition ;[s] 2:0,0;49,1;74,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect] n=5; Clear[x,y,h,z] x = Table[i1,{i1,n}]; y = Table[Random[],{n}]; h = Table[x[[i1+1]]-x[[i1]],{i1,n-1}]; mat = Table[0,{n},{n}]; dat = Table[{0},{n}]; For[k=1,k<=n-2,k++, mat[[k+1,k]] = h[[k]]/6; mat[[k+1,k+1]] = (h[[k]]+h[[k+1]])/3; mat[[k+1,k+2]] = h[[k+1]]/6; dat[[k+1,1]] = -(y[[k]]/h[[k]]) + y[[k+1]]/h[[k]] + y[[k+1]]/h[[k+1]] - y[[k+2]]/h[[k+1]]; dat[[k+1,1]] = -dat[[k+1,1]]; ]; {MatrixForm[mat],MatrixForm[Array[z,n]],MatrixForm[dat]} :[font = text; inactive; Cclosed; preserveAspect; startGroup] Derivation of Rivlin's end-conditions & solve system :[font = text; inactive; Cclosed; preserveAspect; startGroup] Now the left end for s'[x1] :[font = input; preserveAspect; endGroup] temp = Simplify[Lend]; c1 = Coefficient[temp,z1] c2 = Coefficient[temp,z2] leftOver=Simplify[temp-c1 z1 - c2 z2 ] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Now the right end for s'[xn] :[font = input; preserveAspect; endGroup; endGroup] temp = Simplify[Rend]; c1 = Coefficient[Numerator[temp],znm1]/Denominator[temp] c2 = Coefficient[Numerator[temp],zn]/Denominator[temp] leftOver=Simplify[temp-c1 znm1 - c2 zn ] :[font = input; preserveAspect] mat[[1,1]]=h[[1]]/3; mat[[1,2]]=h[[1]]/6; mat[[n,n]]=h[[1]]/3; mat[[n,n-1]]=1/6; dat[[1]]=(y[[2]]-y[[1]])/h[[1]]; dat[[n]]=-(y[[n]]-y[[n-1]])/h[[n-1]]; {MatrixForm[mat],MatrixForm[Array[z,n]],MatrixForm[dat]} :[font = input; preserveAspect] aug = Transpose[Append[Transpose[mat],Flatten[dat]]]; For[k=1,k<=n-1,k++, G = IdentityMatrix[n]; G[[k+1,k]] = -aug[[k+1,k]]/aug[[k,k]]; aug=G.aug; ] MatrixForm[aug] z=Table[0,{n}]; z[[n]]=aug[[n,n+1]]/aug[[n,n]]; For[k=1,k<=n-1,k++, z[[n-k]] = (aug[[n-k,n+1]]-z[[n-k+1]]aug[[n-k,n-k+1]])/ aug[[n-k,n-k]]; ] z = Transpose[{z}]; mat.z-dat :[font = text; inactive; Cclosed; preserveAspect; startGroup] Define the Spline :[font = input; preserveAspect] Clear[s,ds,dds]; z=Flatten[z]; s[k_,t_] := ((-t + x[[k+1]])^3*z[[k]])/(6*h[[k]]) + ((-t + x[[k+1]])*(y[[k]] - (h[[k]]^2*z[[k]])/6))/h[[k]] + ((t - x[[k]])^3*z[[k+1]])/(6*h[[k]]) + ((t - x[[k]])*(y[[k+1]] - (h[[k]]^2*z[[k+1]])/6))/h[[k]]; ds[k_,t_] := D[s[k,tt],tt]/.{tt->t}; dds[k_,t_] := D[ds[k,tt],tt]/.{tt->t}; :[font = input; preserveAspect] Clear[S,dS,ddS]; S[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],s[4,x[[5]]]] dS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],ds[4,x[[5]]]] ddS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],dds[4,x[[5]]]] :[font = input; preserveAspect; endGroup; endGroup] Plot[S[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}] Plot[dS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],ds[n-1,x[[n]]]}], Table[Point[{x[[i1]],ds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] Plot[ddS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],dds[n-1,x[[n]]]}], Table[Point[{x[[i1]],dds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Set up and solve tridiagonal system for natural splines ;[s] 2:0,0;40,1;56,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect] n=5; Clear[x,y,h,z] x = Table[i1,{i1,n}]; y = Table[Random[],{n}]; h = Table[x[[i1+1]]-x[[i1]],{i1,n-1}]; mat = Table[0,{n},{n}]; dat = Table[{0},{n}]; For[k=1,k<=n-2,k++, mat[[k+1,k]] = h[[k]]/6; mat[[k+1,k+1]] = (h[[k]]+h[[k+1]])/3; mat[[k+1,k+2]] = h[[k+1]]/6; dat[[k+1,1]] = -(y[[k]]/h[[k]]) + y[[k+1]]/h[[k]] + y[[k+1]]/h[[k+1]] - y[[k+2]]/h[[k+1]]; dat[[k+1,1]] = -dat[[k+1,1]]; ]; {MatrixForm[mat],MatrixForm[Array[z,n]],MatrixForm[dat]} :[font = input; preserveAspect] mat[[1,1]]=1; mat[[n,n]]=1; dat[[1]]=0; dat[[n]]=0; {MatrixForm[mat],MatrixForm[Array[z,n]],MatrixForm[dat]} :[font = input; preserveAspect] aug = Transpose[Append[Transpose[mat],Flatten[dat]]]; For[k=1,k<=n-1,k++, G = IdentityMatrix[n]; G[[k+1,k]] = -aug[[k+1,k]]/aug[[k,k]]; aug=G.aug; ] MatrixForm[aug] z=Table[0,{n}]; z[[n]]=aug[[n,n+1]]/aug[[n,n]]; For[k=1,k<=n-1,k++, z[[n-k]] = (aug[[n-k,n+1]]-z[[n-k+1]]aug[[n-k,n-k+1]])/ aug[[n-k,n-k]]; ] z = Transpose[{z}]; mat.z-dat :[font = text; inactive; Cclosed; preserveAspect; startGroup] Define the Spline :[font = input; preserveAspect] Clear[s,ds,dds]; z=Flatten[z]; s[k_,t_] := ((-t + x[[k+1]])^3*z[[k]])/(6*h[[k]]) + ((-t + x[[k+1]])*(y[[k]] - (h[[k]]^2*z[[k]])/6))/h[[k]] + ((t - x[[k]])^3*z[[k+1]])/(6*h[[k]]) + ((t - x[[k]])*(y[[k+1]] - (h[[k]]^2*z[[k+1]])/6))/h[[k]]; ds[k_,t_] := D[s[k,tt],tt]/.{tt->t}; dds[k_,t_] := D[ds[k,tt],tt]/.{tt->t}; :[font = input; preserveAspect] Clear[S,dS,ddS]; S[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],s[4,x[[5]]]+ds[4,x[[5]]](t-x[[5]])] dS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],ds[4,x[[5]]]] ddS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],0] :[font = input; preserveAspect; endGroup; endGroup; endGroup] Plot[S[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}] Plot[dS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],ds[n-1,x[[n]]]}], Table[Point[{x[[i1]],ds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] Plot[ddS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],dds[n-1,x[[n]]]}], Table[Point[{x[[i1]],dds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] :[font = text; inactive; preserveAspect; startGroup] Set up and solve tridiagonal system for Rivlin's fundamental splines ;[s] 3:0,0;49,1;68,0;69,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect] n=5; Clear[x,y,h,z] x = Table[i1,{i1,n}]; y = Table[0,{n}]; y[[1]]=1; h = Table[x[[i1+1]]-x[[i1]],{i1,n-1}]; :[font = text; inactive; Cclosed; preserveAspect; startGroup] All the grunt work... :[font = input; preserveAspect; endGroup] mat = Table[0,{n},{n}]; dat = Table[{0},{n}]; For[k=1,k<=n-2,k++, mat[[k+1,k]] = h[[k]]/6; mat[[k+1,k+1]] = (h[[k]]+h[[k+1]])/3; mat[[k+1,k+2]] = h[[k+1]]/6; dat[[k+1,1]] = -(y[[k]]/h[[k]]) + y[[k+1]]/h[[k]] + y[[k+1]]/h[[k+1]] - y[[k+2]]/h[[k+1]]; dat[[k+1,1]] = -dat[[k+1,1]]; ]; mat[[1,1]]=h[[1]]/3; mat[[1,2]]=h[[1]]/6; mat[[n,n]]=h[[1]]/3; mat[[n,n-1]]=1/6; dat[[1]]=(y[[2]]-y[[1]])/h[[1]]; dat[[n]]=-(y[[n]]-y[[n-1]])/h[[n-1]]; aug = Transpose[Append[Transpose[mat],Flatten[dat]]]; For[k=1,k<=n-1,k++, G = IdentityMatrix[n]; G[[k+1,k]] = -aug[[k+1,k]]/aug[[k,k]]; aug=G.aug; ] z=Table[0,{n}]; z[[n]]=aug[[n,n+1]]/aug[[n,n]]; For[k=1,k<=n-1,k++, z[[n-k]] = (aug[[n-k,n+1]]-z[[n-k+1]]aug[[n-k,n-k+1]])/ aug[[n-k,n-k]]; ] z = Transpose[{z}]; Clear[s,ds,dds]; z=Flatten[z]; s[k_,t_] := ((-t + x[[k+1]])^3*z[[k]])/(6*h[[k]]) + ((-t + x[[k+1]])*(y[[k]] - (h[[k]]^2*z[[k]])/6))/h[[k]] + ((t - x[[k]])^3*z[[k+1]])/(6*h[[k]]) + ((t - x[[k]])*(y[[k+1]] - (h[[k]]^2*z[[k+1]])/6))/h[[k]]; ds[k_,t_] := D[s[k,tt],tt]/.{tt->t}; dds[k_,t_] := D[ds[k,tt],tt]/.{tt->t}; Clear[S,dS,ddS]; S[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],s[4,x[[5]]]] dS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],ds[4,x[[5]]]] ddS[t_] := Which[t=x[[1]] && t=x[[2]] && t=x[[3]] && t=x[[4]] && t=x[[5]],dds[4,x[[5]]]] :[font = text; inactive; preserveAspect; startGroup] Plot the fundamental spline :[font = input; preserveAspect; endGroup; endGroup] Plot[S[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,Length[x]}]}] Plot[dS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],ds[n-1,x[[n]]]}], Table[Point[{x[[i1]],ds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] Plot[ddS[t], {t,x[[1]]-1,x[[n]]+1}, Prolog->{PointSize[.02],Point[{x[[n]],dds[n-1,x[[n]]]}], Table[Point[{x[[i1]],dds[i1,x[[i1]]]}],{i1,Length[x]-1}]}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Compare quartic polynomial to 3-segment natural cubic spline :[font = text; inactive; Cclosed; preserveAspect; startGroup] Data to interpolate :[font = input; preserveAspect] {x0,x1,x2,x3,x4} = {0,1,2,3,4}; {y0,y1,y2,y3,y4} = {0,0,0,1,0}; dat = Transpose[{{x0,x1,x2,x3,x4}, {y0,y1,y2,y3,y4}}]; w[x_] := (x-x0)(x-x1)(x-x2)(x-x3)(x-x4) :[font = input; preserveAspect; endGroup] Show[Plot[{w[x], -w[x], 0}, {x, x0 - .2, x4 + .2}, PlotStyle -> {Dashing[{1, 0}], Dashing[{.025, .025}]}, Axes -> False], Table[ Graphics[Line[ {{x0 + i1 (x4 - x0)/40, -Abs[w[x0 + i1 (x4 - x0)/40]]}, {x0 + i1 (x4 - x0)/40, Abs[w[x0 + i1 (x4 - x0)/40]]}}]], {i1, 40}]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Simple interpolating quartic ;[s] 3:0,0;31,1;38,0;39,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] Clear[q,k,q0,q1,q2,q3,q4]; q[x_] := q4 x^4 + q3 x^3 + q2 x^2 + q1 x + q0; temp = Solve[{q[x0] == y0, q[x1] == y1, q[x2] == y2, q[x3] == y3, q[x4] == y4},{q0,q1,q2,q3,q4}]; coefs = {q0,q1,q2,q3,q4} /. Flatten[temp]; {q0,q1,q2,q3,q4}=coefs; pl1 = Plot[{q[t]},{t,x0-.1,x4+.11}, PlotStyle->{Dashing[{1,0}], Dashing[{.02,.02}]}, Prolog-> {PointSize[.02],Table[Point[dat[[i1]]],{i1,5}]}] Print["q[t] = ",q[t]] Print[" "]; Print["Approx bending energy: ", Integrate[q''[t]^2, {t, x0, x4}]," = ", ApproxQuarticEnergy=NIntegrate[q''[t]^2, {t, x0, x4}]] Print[" "]; Print["Bending energy: ", QuarticEnergy=NIntegrate[q''[t]^2/Sqrt[1+q'[t]^2]^3, {t, x0, x4}]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Natural spline interpolate ;[s] 3:0,0;10,1;24,0;37,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Define piecewise polynomial :[font = input; preserveAspect; endGroup] Clear[f,p1,p2,p3, a0,a1,a2,a3, b0,b1,b2,b3, c0,c1,c2,c3, d0,d1,d2,d3]; p1[x_] := a3 x^3 + a2 x^2 + a1 x + a0; p2[x_] := b3 x^3 + b2 x^2 + b1 x + b0; p3[x_] := c3 x^3 + c2 x^2 + c1 x + c0; p4[x_] := d3 x^3 + d2 x^2 + d1 x + d0; f[x_] := Which[x0<=x && x{PointSize[.02], Point[{x0,y0}], Point[{x1,y1}], Point[{x2,y2}], Point[{x3,y3}], Point[{x4,y4}]}] Plot[{f[t],p1[t],p2[t],p3[t],p4[t]},{t,x0-1,x4+1}, Prolog->{PointSize[.02], Point[{x0,y0}], Point[{x1,y1}], Point[{x2,y2}], Point[{x3,y3}], Point[{x4,y4}]}, PlotStyle->{Dashing[{1,0}], Dashing[{.01,.01}], Dashing[{.025,.025}], Dashing[{.005,.005}], Dashing[{.01,.01}]}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Now how does the "energy" compare? :[font = input; preserveAspect; endGroup; endGroup; endGroup] kf[x_] := Which[x0<=x && x{Dashing[{1,0}], Dashing[{.01,.01}]}] Print["Approx bending energy: ", Integrate[p1''[t]^2, {t, x0, x1}] + Integrate[p2''[t]^2, {t, x1, x2}] + Integrate[p3''[t]^2, {t, x2, x3}] + Integrate[p4''[t]^2, {t, x3, x4}]," = ", N[Integrate[p1''[t]^2, {t, x0, x1}] + Integrate[p2''[t]^2, {t, x1, x2}] + Integrate[p3''[t]^2, {t, x2, x3}] + Integrate[p4''[t]^2, {t, x3, x4}]], " < ApproxQuarticEnergy = ",ApproxQuarticEnergy] Print[" "]; Print["Bending energy: ", NIntegrate[kf[t]^2, {t, x0, x4}], " < QuarticEnergy = ",QuarticEnergy] ;[s] 3:0,1;582,0;1111,1;1141,-1; 2:1,12,10,Courier,1,12,0,0,0;2,13,10,MC,1,12,0,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] "not a Knot" spline ;[s] 3:0,0;1,1;11,0;21,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect] Clear[g,p1,p2,p3, a0,a1,a2,a3, b0,b1,b2,b3, c0,c1,c2,c3, d0,d1,d2,d3]; p1[x_] := a3 x^3 + a2 x^2 + a1 x + a0; p2[x_] := b3 x^3 + b2 x^2 + b1 x + b0; p3[x_] := c3 x^3 + c2 x^2 + c1 x + c0; p4[x_] := d3 x^3 + d2 x^2 + d1 x + d0; g[x_] := Which[ x{PointSize[.02], Point[{x0,y0}], Point[{x1,y1}], Point[{x2,y2}], Point[{x3,y3}], Point[{x4,y4}]}] Plot[{g[t],p1[t],p2[t],p3[t],p4[t]},{t,x0-1,x4+1}, Prolog->{PointSize[.02], Point[{x0,y0}], Point[{x1,y1}], Point[{x2,y2}], Point[{x3,y3}], Point[{x4,y4}]}, PlotStyle->{RGBColor[1,1,1], RGBColor[.5,.5,0], RGBColor[0,1,0], RGBColor[0,.5,.5], RGBColor[0,0,1]}] ;[s] 11:0,0;880,1;888,0;914,1;922,0;950,1;958,0;984,1;992,0;1020,1;1028,0;1052,-1; 2:6,12,10,Courier,1,12,0,0,0;5,12,10,Courier,0,12,0,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Compare curvature of natural spline with that of not-a-Knot spline ;[s] 5:0,0;21,1;35,0;49,1;66,0;67,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; startGroup] Show[pl2, pl1] :[font = postscript; PostScript; formatAsPostScript; output; inactive; preserveAspect; pictureLeft = 34; pictureWidth = 282; pictureHeight = 174] %! %%Creator: Mathematica %%AspectRatio: .61803 MathPictureStart %% Graphics /Courier findfont 10 scalefont setfont % Scaling calculations 0.18254 0.15873 0.357521 0.213301 [ [(-1)] .02381 .35752 0 2 Msboxa [(1)] .34127 .35752 0 2 Msboxa [(2)] .5 .35752 0 2 Msboxa [(3)] .65873 .35752 0 2 Msboxa [(4)] .81746 .35752 0 2 Msboxa [(5)] .97619 .35752 0 2 Msboxa [(-1.5)] .17004 .03757 1 0 Msboxa [(-1)] .17004 .14422 1 0 Msboxa [(-0.5)] .17004 .25087 1 0 Msboxa [(0.5)] .17004 .46417 1 0 Msboxa [(1)] .17004 .57082 1 0 Msboxa [ -0.001 -0.001 0 0 ] [ 1.001 .61903 0 0 ] ] MathScale % Start of Graphics 1 setlinecap 1 setlinejoin newpath [ ] 0 setdash 0 g p p .002 w .02381 .35752 m .02381 .36377 L s P [(-1)] .02381 .35752 0 2 Mshowa p .002 w .34127 .35752 m .34127 .36377 L s P [(1)] .34127 .35752 0 2 Mshowa p .002 w .5 .35752 m .5 .36377 L s P [(2)] .5 .35752 0 2 Mshowa p .002 w .65873 .35752 m .65873 .36377 L s P [(3)] .65873 .35752 0 2 Mshowa p .002 w .81746 .35752 m .81746 .36377 L s P [(4)] .81746 .35752 0 2 Mshowa p .002 w .97619 .35752 m .97619 .36377 L s P [(5)] .97619 .35752 0 2 Mshowa p .001 w .05556 .35752 m .05556 .36127 L s P p .001 w .0873 .35752 m .0873 .36127 L s P p .001 w .11905 .35752 m .11905 .36127 L s P p .001 w .15079 .35752 m .15079 .36127 L s P p .001 w .21429 .35752 m .21429 .36127 L s P p .001 w .24603 .35752 m .24603 .36127 L s P p .001 w .27778 .35752 m .27778 .36127 L s P p .001 w .30952 .35752 m .30952 .36127 L s P p .001 w .37302 .35752 m .37302 .36127 L s P p .001 w .40476 .35752 m .40476 .36127 L s P p .001 w .43651 .35752 m .43651 .36127 L s P p .001 w .46825 .35752 m .46825 .36127 L s P p .001 w .53175 .35752 m .53175 .36127 L s P p .001 w .56349 .35752 m .56349 .36127 L s P p .001 w .59524 .35752 m .59524 .36127 L s P p .001 w .62698 .35752 m .62698 .36127 L s P p .001 w .69048 .35752 m .69048 .36127 L s P p .001 w .72222 .35752 m .72222 .36127 L s P p .001 w .75397 .35752 m .75397 .36127 L s P p .001 w .78571 .35752 m .78571 .36127 L s P p .001 w .84921 .35752 m .84921 .36127 L s P p .001 w .88095 .35752 m .88095 .36127 L s P p .001 w .9127 .35752 m .9127 .36127 L s P p .001 w .94444 .35752 m .94444 .36127 L s P p .002 w 0 .35752 m 1 .35752 L s P p .002 w .18254 .03757 m .18879 .03757 L s P [(-1.5)] .17004 .03757 1 0 Mshowa p .002 w .18254 .14422 m .18879 .14422 L s P [(-1)] .17004 .14422 1 0 Mshowa p .002 w .18254 .25087 m .18879 .25087 L s P [(-0.5)] .17004 .25087 1 0 Mshowa p .002 w .18254 .46417 m .18879 .46417 L s P [(0.5)] .17004 .46417 1 0 Mshowa p .002 w .18254 .57082 m .18879 .57082 L s P [(1)] .17004 .57082 1 0 Mshowa p .001 w .18254 .0589 m .18629 .0589 L s P p .001 w .18254 .08023 m .18629 .08023 L s P p .001 w .18254 .10156 m .18629 .10156 L s P p .001 w .18254 .12289 m .18629 .12289 L s P p .001 w .18254 .16555 m .18629 .16555 L s P p .001 w .18254 .18688 m .18629 .18688 L s P p .001 w .18254 .20821 m .18629 .20821 L s P p .001 w .18254 .22954 m .18629 .22954 L s P p .001 w .18254 .2722 m .18629 .2722 L s P p .001 w .18254 .29353 m .18629 .29353 L s P p .001 w .18254 .31486 m .18629 .31486 L s P p .001 w .18254 .33619 m .18629 .33619 L s P p .001 w .18254 .37885 m .18629 .37885 L s P p .001 w .18254 .40018 m .18629 .40018 L s P p .001 w .18254 .42151 m .18629 .42151 L s P p .001 w .18254 .44284 m .18629 .44284 L s P p .001 w .18254 .4855 m .18629 .4855 L s P p .001 w .18254 .50683 m .18629 .50683 L s P p .001 w .18254 .52816 m .18629 .52816 L s P p .001 w .18254 .54949 m .18629 .54949 L s P p .001 w .18254 .01624 m .18629 .01624 L s P p .001 w .18254 .59215 m .18629 .59215 L s P p .001 w .18254 .61348 m .18629 .61348 L s P p .002 w .18254 0 m .18254 .61803 L s P P 0 0 m 1 0 L 1 .61803 L 0 .61803 L closepath clip newpath .02 w .18254 .35752 Mdot .34127 .35752 Mdot .5 .35752 Mdot .65873 .57082 Mdot .81746 .35752 Mdot p p p .004 w .02381 .33467 m .06349 .34038 L .10317 .34609 L .14286 .35181 L .18254 .35752 L .20238 .36033 L .22222 .36288 L .23214 .36397 L .24206 .36489 L .24702 .36527 L .25198 .36561 L .25694 .36588 L .25942 .36599 L .2619 .36609 L .26438 .36617 L .26687 .36624 L .26811 .36626 L .26935 .36628 L .27059 .3663 L .27183 .36631 L .27307 .36632 L .27431 .36632 L .27555 .36631 L .27679 .36631 L .27803 .36629 L .27927 .36628 L .28175 .36623 L .28423 .36615 L .28671 .36606 L .29167 .36581 L .29663 .36546 L .30159 .36502 L .31151 .36383 L .32143 .36221 L .33135 .36012 L .34127 .35752 L .38095 .34359 L .40079 .33677 L .41071 .33397 L .41567 .3328 L .42063 .33181 L .4256 .33102 L .42808 .33071 L .43056 .33045 L .4318 .33035 L .43304 .33026 L .43428 .33019 L .43552 .33013 L .43676 .33009 L .438 .33007 L Mistroke .43924 .33006 L .44048 .33007 L .44172 .3301 L .44296 .33014 L .4442 .33021 L .44544 .33029 L .44792 .33052 L .4504 .33082 L .45536 .33168 L .46032 .33288 L .46528 .33445 L .47024 .33641 L .48016 .34159 L .49008 .34857 L .5 .35752 L .51984 .38137 L .53968 .4112 L .57937 .47846 L .59921 .5107 L .61905 .53857 L .62897 .55005 L .63889 .55947 L .64385 .5633 L .64881 .5665 L .65377 .56902 L .65625 .57001 L .65873 .57082 L .65997 .57116 L .66121 .57144 L .66245 .57168 L .66369 .57188 L .66493 .57202 L .66617 .57212 L .66741 .57218 L .66865 .57219 L .66989 .57216 L .67113 .57208 L .67237 .57196 L .67361 .5718 L .67609 .57134 L .67857 .57072 L .68353 .56897 L .68849 .56659 L .69841 .55999 L .70833 .55112 L .71825 .54016 L .7381 .51274 L .77778 .4412 L .81746 .35752 L .85714 .27182 L Mistroke .89683 .18612 L .93651 .10042 L .97619 .01472 L Mfstroke P P p p [ 1 0 ] 0 setdash p .004 w .16667 .32385 m .18059 .35395 L .19451 .37624 L .20147 .38478 L .20843 .39173 L .21539 .39721 L .21887 .39943 L .22235 .40133 L .22583 .40292 L .22932 .40421 L .2328 .40522 L .23454 .40562 L .23628 .40595 L .23802 .40622 L .23889 .40633 L .23976 .40642 L .24063 .4065 L .2415 .40657 L .24237 .40662 L .24324 .40665 L .24411 .40667 L .24498 .40668 L .24585 .40667 L .24672 .40665 L .24759 .40661 L .24846 .40656 L .24933 .4065 L .2502 .40642 L .25194 .40623 L .25368 .40598 L .25716 .40535 L .26064 .40452 L .26412 .40353 L .27108 .40105 L .27804 .39801 L .30589 .38179 L .33373 .36266 L .36157 .34467 L .3755 .33717 L .38942 .33112 L .39638 .32873 L .40334 .3268 L .4103 .32536 L .41378 .32484 L .41552 .32463 L .41726 .32445 L .419 .3243 L .41987 .32424 L .42074 .32419 L .42161 .32414 L Mistroke .42248 .32411 L .42335 .32408 L .42422 .32406 L .42509 .32406 L .42596 .32406 L .42683 .32406 L .4277 .32408 L .42857 .32411 L .42944 .32414 L .43031 .32419 L .43118 .32424 L .43466 .32454 L .4364 .32474 L .43814 .32499 L .44163 .32558 L .44511 .32631 L .45207 .32823 L .45903 .33074 L .47295 .33756 L .48687 .34675 L .50079 .35824 L .52864 .38756 L .55648 .42395 L .58433 .46504 L .61217 .50761 L .64001 .54768 L .65394 .5653 L .66786 .58042 L .68178 .59233 L .68874 .59683 L .69222 .59867 L .6957 .60023 L .69918 .60149 L .70092 .602 L .70266 .60243 L .7044 .60277 L .70527 .60292 L .70614 .60304 L .70701 .60314 L .70788 .60322 L .70875 .60327 L .70962 .60331 L .71049 .60332 L .71136 .60331 L .71223 .60328 L .7131 .60322 L .71397 .60314 L .71484 .60304 L .71658 .60276 L .71832 .60239 L Mistroke .72006 .60192 L .72354 .60068 L .72703 .59902 L .73051 .59693 L .73747 .59141 L .74443 .58398 L .75139 .57453 L .76531 .549 L .77923 .51374 L .79315 .46762 L .80708 .40947 L .83492 .25205 L Mfstroke P P P P % End of Graphics MathPictureEnd :[font = output; output; inactive; preserveAspect; endGroup; endGroup; endGroup] Graphics["<<>>"] ;[o] -Graphics- :[font = section; inactive; Cclosed; preserveAspect; startGroup] B-Splines :[font = text; inactive; Cclosed; preserveAspect; startGroup] Divided difference approach to B-splines using one-sided ramp functions :[font = text; inactive; Cclosed; preserveAspect; startGroup] 1 (u3-u1)*[u1,u2,u3](x-t) + :[font = input; preserveAspect; endGroup] Off[Plot::plnr]; Clear[ramp,u,x]; ramp[i_,j_,x_] := If[i>=j,0,(x-i)] u = {0,1,2}; pieces = {}; For[i=2,i<=3,i++, Print[" "]; d0 = Expand[{ramp[u[[1]],u[[i]],x], ramp[u[[2]],u[[i]],x], ramp[u[[3]],u[[i]],x]}]; d1 = Expand[{(d0[[2]]-d0[[1]])/(u[[2]]-u[[1]]), (d0[[3]]-d0[[2]])/(u[[3]]-u[[2]])}]; d2 = Expand[{(d1[[2]]-d1[[1]])/(u[[3]]-u[[1]])}]; pieces = Append[pieces,d2[[1]]]; ] pieces =Expand[ (u[[3]]-u[[1]])*pieces]; Print["pieces = ",pieces]; Plot[Which[ x{AbsoluteThickness[2]}, PlotRange->{0,1}] ;[s] 3:0,0;4,1;14,0;755,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] 2 - (u4-u1)*[u1,u2,u3,u4](x-t) + :[font = input; preserveAspect; endGroup] Clear[ramp,u]; ramp[i_,j_,x_] := If[i>=j,0,(x-i)^2] u = {0,1,2,3}; pieces = {}; For[i=2,i<=4,i++, Print[" "]; d0 = Expand[{ramp[u[[1]],u[[i]],x], ramp[u[[2]],u[[i]],x], ramp[u[[3]],u[[i]],x], ramp[u[[4]],u[[i]],x]}]; d1 = Expand[{(d0[[2]]-d0[[1]])/(u[[2]]-u[[1]]), (d0[[3]]-d0[[2]])/(u[[3]]-u[[2]]), (d0[[4]]-d0[[3]])/(u[[4]]-u[[3]])}]; d2 = Expand[{(d1[[2]]-d1[[1]])/(u[[3]]-u[[1]]), (d1[[3]]-d1[[2]])/(u[[4]]-u[[2]])}]; d3 = Expand[{(d2[[2]]-d2[[1]])/(u[[4]]-u[[1]])}]; pieces = Append[pieces,d3[[1]]]; ] pieces = Expand[-(u[[4]]-u[[1]])*pieces]; Print["pieces = ",pieces]; Plot[Which[ x{AbsoluteThickness[2]}, PlotRange->{-1,1}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] 3 (u5-u1)*[u1,u2,u3,u4,u5](x-t) + :[font = input; preserveAspect; endGroup; endGroup] Clear[ramp,u]; ramp[i_,j_,x_] := If[i>=j,0,(x-i)^3] u = {0,1,2,3,4}; pieces = {}; For[i=2,i<=5,i++, Print[" "]; d0 = Expand[{ramp[u[[1]],u[[i]],x], ramp[u[[2]],u[[i]],x], ramp[u[[3]],u[[i]],x], ramp[u[[4]],u[[i]],x], ramp[u[[5]],u[[i]],x]}]; d1 = Expand[{(d0[[2]]-d0[[1]])/(u[[2]]-u[[1]]), (d0[[3]]-d0[[2]])/(u[[3]]-u[[2]]), (d0[[4]]-d0[[3]])/(u[[4]]-u[[3]]), (d0[[5]]-d0[[4]])/(u[[5]]-u[[4]])}]; d2 = Expand[{(d1[[2]]-d1[[1]])/(u[[3]]-u[[1]]), (d1[[3]]-d1[[2]])/(u[[4]]-u[[2]]), (d1[[4]]-d1[[3]])/(u[[5]]-u[[3]])}]; d3 = Expand[{(d2[[2]]-d2[[1]])/(u[[4]]-u[[1]]), (d2[[3]]-d2[[2]])/(u[[5]]-u[[2]])}]; d4 = Expand[{(d3[[2]]-d3[[1]])/(u[[5]]-u[[1]])}]; pieces = Append[pieces,d4[[1]]]; ] pieces = Expand[(u[[5]]-u[[1]])*pieces]; Print["pieces = ",pieces]; Plot[Which[ x{AbsoluteThickness[2]}, PlotRange->{0,1}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] deBoor's BSPLVB ;[s] 2:0,1;15,0;16,-1; 2:1,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = text; inactive; preserveAspect; startGroup] Define the general B-spline bsp[i,k,u,x],where k is the order; i=1,...,dim is the index of the splines; u is the knot set; x is where the spline is to be evaluated; dim =length(u)-k. Note that the components of the B-splines are right-continuous.Thus all are defined in intervals[a,b) except for the last one which is in[a,b]. ;[s] 13:0,0;28,1;40,0;53,1;54,0;75,1;86,0;122,1;123,0;147,1;148,0;195,1;211,0;357,-1; 2:7,13,9,Times,0,12,0,0,0;6,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] Off[General::spell1] Clear[bsp] bsp[index_, order_, u_, x_] := Block[ {j, biatx, deltar, deltal, jp1, i, il, iu, im, term, saved, left}, j = 1; biatx = Table[0, {order}]; deltar = Table[0, {order}]; deltal = Table[0, {order}]; dim = Length[u] - order; biatx[[1]] = 1.; (*Determine which interval x is in with binary search*) il = 0; iu = dim + order; While[iu - il > 1, im = Floor[(il + iu)/2]; If[u[[dim + order]] > u[[1]] && x >= u[[im]], il = im;, iu = im; ]; ]; left = il; (*If x is the last point, do the right thing.That is, remember that splines are defined in half - open intervals, except for the last one which is closed.*) If[u[[dim + order]] == x, For[jj = 0, jj < order, jj++, If[u[[dim + order - jj - 1]] < x && x <= u[[dim + order - jj]], left = dim + order - jj - 1]; ]; ]; (*... but only report the relevant value*) If[(left - index) < 0 || (left - index) >= order || x > u[[dim + order]] || x < u[[1]], Return[0.], (*Basically do BSPLVB*) While[j < order, jp1 = j + 1; deltar[[j]] = u[[left + j]] - x; deltal[[j]] = x - u[[left + 1 - j]]; saved = 0.; For[i = 1, i <=j, i++, term = biatx[[i]]/(deltar[[i]] + deltal[[jp1 - i]]); biatx[[i]] = saved + deltar[[i]]*term; saved = deltal[[jp1 - i]]*term; ]; biatx[[jp1]] = saved; j = jp1; ]; Return[biatx[[order - (left - index)]]]; ]; ] ;[s] 9:0,0;4,1;19,0;306,2;366,0;587,2;762,0;1016,2;1063,0;1729,-1; 3:5,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,65535,0,0;3,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot the uniform cubic B-splines :[font = input; preserveAspect; endGroup] order = 4; u = {-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7}; dim = Length[u] - order; Plot[{bsp[1, order, u, s], bsp[2, order, u, s], bsp[3, order, u, s], bsp[4, order, u, s], bsp[5, order, u, s], bsp[6, order, u, s], bsp[7, order, u, s]}, {s, u[[order]], u[[dim + 1]]}, Prolog -> {PointSize[.02], Table[Point[{u[[i1]], 0}], {i1, Length[u]}]}, PlotRange -> {{First[u], Last[u]}, {0, 1}}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot the psuedo-uniform cubic B-splines :[font = input; preserveAspect; endGroup] order = 4; u = {0, 0, 0, 0, 1, 2, 3, 4, 4, 4, 4}; dim = Length[u] - order; Plot[{bsp[1, order, u, s], bsp[2, order, u, s], bsp[3, order, u, s], bsp[4, order, u, s], bsp[5, order, u, s], bsp[6, order, u, s], bsp[7, order, u, s]}, {s, u[[order]], u[[dim + 1]]}, Prolog -> {PointSize[.02], Table[If[u[[i1]] == First[u] || u[[i1]] == Last[u], Point[{u[[i1]], Which[ i1 == 1 || i1 == 9, -.3, i1 == 2 || i1 == 10, -.2, i1 == 3 || i1 == 11, -.1, i1 == 4 || i1 == 8, 0.]}], Point[{u[[i1]], 0}]], {i1, Length[u]}]}, PlotRange -> {{First[u], Last[u]}, {-.4, 1}}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot the nonuniform cubic B-splines :[font = input; preserveAspect; endGroup] order = 4; u = Sort[{0, 0, 0, 0, 4Random[], 4Random[], 4Random[], 4, 4, 4, 4}]; dim = Length[u] - order; Plot[{bsp[1, order, u, s], bsp[2, order, u, s], bsp[3, order, u, s], bsp[4, order, u, s], bsp[5, order, u, s], bsp[6, order, u, s], bsp[7, order, u, s]}, {s, u[[1]], u[[dim + order]]}, Prolog -> {PointSize[.02], Table[If[u[[i1]] == First[u] || u[[i1]] == Last[u], Point[{u[[i1]], Which[ i1 == 1 || i1 == 9, -.3, i1 == 2 || i1 == 10, -.2, i1 == 3 || i1 == 11, -.1, i1 == 4 || i1 == 8, 0.]}], Point[{u[[i1]], 0}]], {i1, Length[u]}]}, PlotRange -> {{First[u], Last[u]}, {-.4, 1}}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot the cubic Bezier polynomials: :[font = input; preserveAspect; endGroup; endGroup] order = 4; u = {0,0,0,0,1,1,1,1}; dim = Length[u]-order; Plot[{bsp[1,order,u,s], bsp[2,order,u,s], bsp[3,order,u,s], bsp[4,order,u,s]}, {s,u[[order]],u[[dim+1]]}, Prolog->{PointSize[.02], Table[If[u[[i1]]==First[u] || u[[i1]]==Last[u], Point[{u[[i1]], Which[ i1==1 || i1==5,-.3, i1==2 || i1==6,-.2, i1==3 || i1==7,-.1, i1==4 || i1==8,0.]}], Point[{u[[i1]],0}]],{i1,Length[u]}]}, PlotRange->{{First[u],Last[u]},{-.4,1}}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Splines and derivatives :[font = text; inactive; Cclosed; dontPreserveAspect; startGroup] Define the general B-spline bsp[i,k,u,x], and derivatives dbsp[i,k,u,x], ddbsp[i,k,u,x], dddbsp[i,k,u,x], where k is the order; i=1,...,dim is the index of the splines; u is the knot set; x is where the spline is to be evaluated; dim = length(u)-k. Note that the components of the B-splines are right-continuous. Thus all are defined in intervals [a,b) except for the last one which is in [a,b]. ;[s] 19:0,0;28,1;40,0;58,1;104,0;118,1;119,0;139,1;140,0;185,1;186,0;209,1;210,0;256,1;259,0;269,1;270,0;272,1;273,0;426,-1; 2:10,13,9,Times,0,12,0,0,0;9,13,9,Times,1,12,0,0,0; :[font = input; initialization; dontPreserveAspect; startGroup] *) bsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, biatx = Table[0,{d1,k}], deltar = Table[0,{d1,k}], deltal = Table[0,{d1,k}], j = 1, jp1 = 0, i1 = 0, term = 0, saved = 0, left = 0, vacuous = 0, outside = 1 }, biatx[[1]] = 1.; (* If interval of support is empty, then vacuous *) If[knots[[index+order]]-knots[[index]]<10^-6, vacuous = 1 ]; (* Determine which subinterval xx is in *) For[i1=index,i10 || vacuous>0, (* Then ... *) 0, (* Else ... Basically do BSPLVB *) While[j=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; For[i1=1,i1<=2,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]; b = 1/delta[[2]]; result = a*bsp[index+1,order-1,u,xx]+ b*bsp[index,order-1,u,xx]; result = (order-1)*result ] :[font = text; inactive; dontPreserveAspect] Define the second derivative of general B-spline ddbsp[i,k,u,x]. ;[s] 3:0,0;49,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect] ddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,5}], a = 0,b = 0,c = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; For[i1=1,i1<=5,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = 1/delta[[1]]/delta[[3]]; b = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]; c = 1/delta[[2]]/delta[[5]]; result = a*bsp[index+2,order-2,u,xx]+ b*bsp[index+1,order-2,u,xx]+ c*bsp[index,order-2,u,xx]; result = (order-1)*(order-2)*result ]; :[font = text; inactive; dontPreserveAspect] Define the third derivative of general B-spline dddbsp[i,k,u,x]. ;[s] 3:0,0;48,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect; endGroup] dddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,9}], a = 0,b = 0,c = 0,d = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; If[index>=-2 && index<=dim, delta[[6]] = u[[index+order]]-u[[index+3]]]; If[index>=-1 && index<=dim+1, delta[[7]] = u[[index+order-1]]-u[[index+2]]]; If[index>=0 && index<=dim+2, delta[[8]] = u[[index+order-2]]-u[[index+1]]]; If[index>=1 && index<=dim+3, delta[[9]] = u[[index+order-3]]-u[[index]]]; For[i1=1,i1<=9,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]/delta[[3]]/delta[[6]]; b = (1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[7]]+ 1/delta[[1]]/delta[[3]]/delta[[7]]; c = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[8]]- 1/delta[[2]]/delta[[5]]/delta[[8]]; d = 1/delta[[2]]/delta[[5]]/delta[[9]]; result = a*bsp[index+3,order-3,u,xx]+ b*bsp[index+2,order-3,u,xx]+ c*bsp[index+1,order-3,u,xx]+ d*bsp[index,order-3,u,xx]; result = (order-1)*(order-2)*(order-3)*result ]; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot family of splines and respective derivatives... :[font = input; preserveAspect] u = {0,0,0,0,1,2,3,3,3,3}; :[font = input; preserveAspect] Plot[{bsp[1,4,u,s], bsp[2,4,u,s], bsp[3,4,u,s], bsp[4,4,u,s], bsp[5,4,u,s], bsp[6,4,u,s]}, {s,0,2.9999}] :[font = input; preserveAspect] Plot[{dbsp[1,4,u,s], dbsp[2,4,u,s], dbsp[3,4,u,s], dbsp[4,4,u,s], dbsp[5,4,u,s], dbsp[6,4,u,s]}, {s,0,2.9999}] :[font = input; preserveAspect] Plot[{ddbsp[1,4,u,s], ddbsp[2,4,u,s], ddbsp[3,4,u,s], ddbsp[4,4,u,s], ddbsp[5,4,u,s], ddbsp[6,4,u,s]}, {s,0,2.9999}] :[font = input; preserveAspect; endGroup] Plot[{dddbsp[1,4,u,s], dddbsp[2,4,u,s], dddbsp[3,4,u,s], dddbsp[4,4,u,s], dddbsp[5,4,u,s], dddbsp[6,4,u,s]}, {s,0,2.9999}] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Plot individual spline and respective derivatives :[font = input; preserveAspect] u = {0,0,0,0,1,2,3,3,3,3}; kk=4; :[font = text; inactive; Cclosed; preserveAspect; startGroup] The grunt work... ;[s] 3:0,0;10,1;15,0;33,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,10,Courier,1,12,0,0,0; :[font = input; preserveAspect; endGroup] max=1.1 Max[Abs[Table[dbsp[kk,4,u,(u[[i1]]+u[[i1+1]])/2],{i1,3,6}]]]; pic1=Plot[{bsp[kk,4,u,s], dbsp[kk,4,u,s]}, {s,0,2.9999}, PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]}, PlotRange->max{-1,1}]; max=1.1 Max[Abs[Table[ddbsp[kk,4,u,(u[[i1]]+u[[i1+1]])/2],{i1,3,6}]]]; pic2=Plot[{dbsp[kk,4,u,s], ddbsp[kk,4,u,s]}, {s,0,2.9999}, PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]}, PlotRange->max{-1,1}]; max=1.1 Max[Abs[Table[dddbsp[kk,4,u,(u[[i1]]+u[[i1+1]])/2],{i1,3,6}]]]; pic3=Plot[{ddbsp[kk,4,u,s], dddbsp[kk,4,u,s]}, {s,0,2.9999}, PlotStyle->{Dashing[{1,0}],Dashing[{.02,.02}]}, PlotRange->max{-1,1}]; :[font = input; preserveAspect; endGroup; endGroup] Show[GraphicsArray[{pic1,pic2,pic3}]] :[font = text; inactive; Cclosed; preserveAspect; startGroup] Now, to fit some curve data... :[font = text; inactive; Cclosed; dontPreserveAspect; startGroup] Define the general B-spline bsp[i,k,u,x], where k is the order; i=1,...,dim is the index of the splines; u is the knot set; x is where the spline is to be evaluated; dim = length(u)-k. Note that the components of the B-splines are right-continuous. Thus all are defined in intervals [a,b) except for the last one which is in [a,b]. ;[s] 17:0,0;28,1;40,0;54,1;55,0;75,1;76,0;121,1;122,0;145,1;146,0;192,1;195,0;205,1;206,0;208,1;209,0;362,-1; 2:9,13,9,Times,0,12,0,0,0;8,13,9,Times,1,12,0,0,0; :[font = input; initialization; dontPreserveAspect] *) bsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, biatx = Table[0,{d1,k}], deltar = Table[0,{d1,k}], deltal = Table[0,{d1,k}], j = 1, jp1 = 0, i1 = 0, term = 0, saved = 0, left = 0, vacuous = 0, outside = 1 }, biatx[[1]] = 1.; (* If interval of support is empty, then vacuous *) If[knots[[index+order]]-knots[[index]]<10^-6, vacuous = 1 ]; (* Determine which subinterval xx is in *) For[i1=index,i10 || vacuous>0, (* Then ... *) 0, (* Else ... Basically do BSPLVB *) While[j=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; For[i1=1,i1<=2,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]; b = 1/delta[[2]]; result = a*bsp[index+1,order-1,u,xx]+ b*bsp[index,order-1,u,xx]; result = (order-1)*result ] :[font = text; inactive; dontPreserveAspect] Define the second derivative of general B-spline ddbsp[i,k,u,x]. ;[s] 3:0,0;49,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect] ddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,5}], a = 0,b = 0,c = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; For[i1=1,i1<=5,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = 1/delta[[1]]/delta[[3]]; b = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]; c = 1/delta[[2]]/delta[[5]]; result = a*bsp[index+2,order-2,u,xx]+ b*bsp[index+1,order-2,u,xx]+ c*bsp[index,order-2,u,xx]; result = (order-1)*(order-2)*result ]; :[font = text; inactive; dontPreserveAspect] Define the third derivative of general B-spline dddbsp[i,k,u,x]. ;[s] 3:0,0;48,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect] dddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,9}], a = 0,b = 0,c = 0,d = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; If[index>=-2 && index<=dim, delta[[6]] = u[[index+order]]-u[[index+3]]]; If[index>=-1 && index<=dim+1, delta[[7]] = u[[index+order-1]]-u[[index+2]]]; If[index>=0 && index<=dim+2, delta[[8]] = u[[index+order-2]]-u[[index+1]]]; If[index>=1 && index<=dim+3, delta[[9]] = u[[index+order-3]]-u[[index]]]; For[i1=1,i1<=9,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]/delta[[3]]/delta[[6]]; b = (1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[7]]+ 1/delta[[1]]/delta[[3]]/delta[[7]]; c = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[8]]- 1/delta[[2]]/delta[[5]]/delta[[8]]; d = 1/delta[[2]]/delta[[5]]/delta[[9]]; result = a*bsp[index+3,order-3,u,xx]+ b*bsp[index+2,order-3,u,xx]+ c*bsp[index+1,order-3,u,xx]+ d*bsp[index,order-3,u,xx]; result = (order-1)*(order-2)*(order-3)*result ]; :[font = input; preserveAspect; endGroup] rangePlus[pts_] := Block[ {m,ptsT,min,max}, m=Length[pts[[1]]]; ptsT=Transpose[pts]; min = Table[Min[ptsT[[i1]]],{i1,m}]; max = Table[Max[ptsT[[i1]]],{i1,m}]; Return[Transpose[{min,max}]];] :[font = text; inactive; Cclosed; preserveAspect; startGroup] In the following example we fit conventional ordered pairs (xi,yi), i=1,...,dim ;[s] 3:0,0;32,1;44,0;80,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] u={-3,-2,-1,0,1,2,3,4,5,6,7}; k=4; n=Length[u]; dim=n-k; x=Table[u[[k]]+(i1-1)/(dim-1)(u[[n-k+1]]-u[[k]]),{i1,dim}]; y=Table[Random[],{dim}]; mat = Table[Table[bsp[j1,k,u,x[[i1]]],{j1,dim}],{i1,dim}]; dat = Transpose[{y}]; coef=Inverse[mat].dat; curve[x_] := Sum[coef[[i1,1]] bsp[i1,k,u,x],{i1,dim}]; Plot[curve[x],{x,u[[k]],u[[n-k+1]]}, Prolog->{PointSize[.02], Table[Point[{x[[i1]],y[[i1]]}],{i1,dim}]} ] ;[s] 2:0,1;35,0;414,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] In the following example we fit ordered pairs (xi,yi), i=1,...,dim, parametrically. ;[s] 3:0,0;68,1;82,0;84,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup] u={-3,-2,-1,0,1,2,3,4}; u={0,0,0,0,1,1,1,1}; k=4; n=Length[u]; dim=n-k; t=Table[u[[k]]+(i1-1)/(dim-1)(u[[n-k+1]]-u[[k]]),{i1,dim}]; xy=Table[{Random[],Random[]},{dim}]; mat = Table[Table[bsp[j1,k,u,t[[i1]]],{j1,dim}],{i1,dim}]; dat = xy; coef=Inverse[mat].dat Clear[curve]; curve[t_] := Sum[coef[[i1]] bsp[i1,k,u,t],{i1,dim}]; Off[ParametricPlot::ppcom]; ParametricPlot[curve[t],{t,u[[k]],u[[n-k+1]]}, Prolog->{PointSize[.02], Table[Point[xy[[i1]]],{i1,dim}]} ] Show[%, Graphics[Line[ Table[coef[[i1]],{i1,dim}]]], Table[Graphics[Text["O",coef[[i1]]]],{i1,dim}], PlotRange->rangePlus[coef]] ;[s] 4:0,1;50,0;334,2;355,0;604,-1; 3:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;1,12,10,Courier,0,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] In the following example we fit ordered triples (xi,yi,zi), i=1,...,dim, parametrically. ;[s] 5:0,0;40,1;47,0;73,1;87,0;89,-1; 2:3,13,9,Times,0,12,0,0,0;2,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup; endGroup] u={-3,-2,-1,0,1,2,3,4}; u={0,0,0,0,1,1,1,1}; k=4; n=Length[u]; dim=n-k; t=Table[u[[k]]+(i1-1)/(dim-1)(u[[n-k+1]]-u[[k]]),{i1,dim}]; xyz=Table[{Random[],Random[],Random[]},{dim}]; mat = Table[Table[bsp[j1,k,u,t[[i1]]],{j1,dim}],{i1,dim}]; dat = xyz; coef=Inverse[mat].dat Clear[curve]; curve[t_] := Sum[coef[[i1]] bsp[i1,k,u,t],{i1,dim}]; Off[ParametricPlot3D::ppcom]; ParametricPlot3D[curve[t],{t,u[[k]],u[[n-k+1]]}] Show[%, Graphics3D[{PointSize[.02], Table[Point[xyz[[i1]]],{i1,dim}]}], Graphics3D[Line[ Table[coef[[i1]],{i1,dim}]]], Table[Graphics3D[Text["O",coef[[i1]]]],{i1,dim}], PlotRange->rangePlus[coef]] ;[s] 4:0,1;50,0;345,2;368,0;626,-1; 3:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;1,12,10,Courier,0,12,65535,0,0; :[font = text; inactive; Cclosed; preserveAspect; startGroup] Now, to fit some surface data... :[font = text; inactive; Cclosed; dontPreserveAspect; startGroup] Define the general B-spline bsp[i,k,u,x], where k is the order; i=1,...,dim is the index of the splines; u is the knot set; x is where the spline is to be evaluated; dim = length(u)-k. Note that the components of the B-splines are right-continuous. Thus all are defined in intervals [a,b) except for the last one which is in [a,b]. ;[s] 17:0,0;28,1;40,0;54,1;55,0;75,1;76,0;121,1;122,0;145,1;146,0;192,1;195,0;205,1;206,0;208,1;209,0;362,-1; 2:9,13,9,Times,0,12,0,0,0;8,13,9,Times,1,12,0,0,0; :[font = input; initialization; dontPreserveAspect] *) bsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, biatx = Table[0,{d1,k}], deltar = Table[0,{d1,k}], deltal = Table[0,{d1,k}], j = 1, jp1 = 0, i1 = 0, term = 0, saved = 0, left = 0, vacuous = 0, outside = 1 }, biatx[[1]] = 1.; (* If interval of support is empty, then vacuous *) If[knots[[index+order]]-knots[[index]]<10^-6, vacuous = 1 ]; (* Determine which subinterval xx is in *) For[i1=index,i10 || vacuous>0, (* Then ... *) 0, (* Else ... Basically do BSPLVB *) While[j=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; For[i1=1,i1<=2,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]; b = 1/delta[[2]]; result = a*bsp[index+1,order-1,u,xx]+ b*bsp[index,order-1,u,xx]; result = (order-1)*result ] :[font = text; inactive; dontPreserveAspect] Define the second derivative of general B-spline ddbsp[i,k,u,x]. ;[s] 3:0,0;49,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect] ddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,5}], a = 0,b = 0,c = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; For[i1=1,i1<=5,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = 1/delta[[1]]/delta[[3]]; b = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]; c = 1/delta[[2]]/delta[[5]]; result = a*bsp[index+2,order-2,u,xx]+ b*bsp[index+1,order-2,u,xx]+ c*bsp[index,order-2,u,xx]; result = (order-1)*(order-2)*result ]; :[font = text; inactive; dontPreserveAspect] Define the third derivative of general B-spline dddbsp[i,k,u,x]. ;[s] 3:0,0;48,1;63,0;65,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; dontPreserveAspect] dddbsp[i_,k_,u_,x_] := Block[ (* Declare and initialize local variables *) { index = i, order = k, knots = u, xx = x, dim = Length[u]-k, delta = Table[0,{i1,1,9}], a = 0,b = 0,c = 0,d = 0 }, Clear[i1]; If[index>=0 && index<=dim, delta[[1]] = u[[index+order]]-u[[index+1]]]; If[index>=1 && index<=dim+1, delta[[2]] = u[[index+order-1]]-u[[index]]]; If[index>=-1 && index<=dim, delta[[3]] = u[[index+order]]-u[[index+2]]]; If[index>=0 && index<=dim+1, delta[[4]] = u[[index+order-1]]-u[[index+1]]]; If[index>=1 && index<=dim+2, delta[[5]] = u[[index+order-2]]-u[[index]]]; If[index>=-2 && index<=dim, delta[[6]] = u[[index+order]]-u[[index+3]]]; If[index>=-1 && index<=dim+1, delta[[7]] = u[[index+order-1]]-u[[index+2]]]; If[index>=0 && index<=dim+2, delta[[8]] = u[[index+order-2]]-u[[index+1]]]; If[index>=1 && index<=dim+3, delta[[9]] = u[[index+order-3]]-u[[index]]]; For[i1=1,i1<=9,i1++, If[delta[[i1]]==0, delta[[i1]] = 1]; ]; a = -1/delta[[1]]/delta[[3]]/delta[[6]]; b = (1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[7]]+ 1/delta[[1]]/delta[[3]]/delta[[7]]; c = -(1/delta[[1]]+1/delta[[2]])/delta[[4]]/delta[[8]]- 1/delta[[2]]/delta[[5]]/delta[[8]]; d = 1/delta[[2]]/delta[[5]]/delta[[9]]; result = a*bsp[index+3,order-3,u,xx]+ b*bsp[index+2,order-3,u,xx]+ c*bsp[index+1,order-3,u,xx]+ d*bsp[index,order-3,u,xx]; result = (order-1)*(order-2)*(order-3)*result ]; :[font = input; preserveAspect; endGroup] rangePlus[pts_] := Block[ {m,ptsT,min,max}, m=Length[pts[[1]]]; ptsT=Transpose[pts]; min = Table[Min[ptsT[[i1]]],{i1,m}]; max = Table[Max[ptsT[[i1]]],{i1,m}]; Return[Transpose[{min,max}]];] :[font = text; inactive; Cclosed; preserveAspect; startGroup] In the following example we fit ordered triples (xi,yi,f[xi,yi]), i=1,...,dim, parametrically. ;[s] 3:0,0;79,1;93,0;95,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = input; preserveAspect; endGroup; endGroup; endGroup] u={0,0,0,0,1,2,2,2,2}; v=u; k=4; f[x_,y_] := N[Cos[x^2+y^2]] nu=Length[u]; nv=Length[v]; dimu=nu-k; dimv=nv-k; dim = dimu*dimv; x=Table[u[[k]]+(i1-1)/(dimu-1)(u[[nu-k+1]]-u[[k]]),{i1,dimu}] y=Table[v[[k]]+(j1-1)/(dimv-1)(v[[nv-k+1]]-v[[k]]),{j1,dimv}] xy=Table[{x[[Mod[l1-1,dimu]+1]],y[[Floor[(l1-1)/dimv]+1]]},{l1,dim}]; z=Table[{f[x[[Mod[l1-1,dimu]+1]],y[[Floor[(l1-1)/dimv]+1]]]},{l1,dim}]; Bsp[l_,k_,u_,v_,s_,t_]:= bsp[Mod[l-1,dimu]+1,k,u,s] * bsp[Floor[(l-1)/dimv]+1,k,v,t] mat = Table[Table[ Bsp[j1,k,u,v,x[[Mod[i1-1,dimu]+1]],y[[Floor[(i1-1)/dimv]+1]]], {j1,dim}],{i1,dim}]; dat = z; coef=Inverse[mat].dat Clear[surface]; surface[s_,t_] := Sum[coef[[l1,1]] Bsp[l1,k,u,v,s,t],{l1,dim}]; Plot3D[surface[s,t],{s,u[[k]],u[[nu-k+1]]},{t,v[[k]],v[[nv-k+1]]} ] Plot3D[f[s,t],{s,u[[k]],u[[nu-k+1]]},{t,v[[k]],v[[nv-k+1]]} ] ;[s] 2:0,1;60,0;844,-1; 2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; ^*)