样条插值及应用深入研究

学院： 研究生学院 专业： 机械工程 组号： 39 成绩：

第三章 算法应用

一．《样条插值及应用》方法怎么用？（程序设计）？

3.1.1 一般程序设计

三弯矩插值法程序:

Clear[x,y,a,b,c,n,M] x[1]=27.7; x[2]=28; x[3]=29; x[4]=30; y[1]=4.1; y[2]=4.3; y[3]=4.1; y[4]=3.0;

B=Table[{x[i],y[i]},{i,1,4}]; y'[1]=3.0; y'[4]=-4.0; h[1]=0.3; h[2]=1; h[3]=1;

a[2]=h[1]/(h[1]+h[2]); a[3]=h[2]/(h[2]+h[3]); a[4]=1; b[1]=1; b[2]=1-a[2]; b[3]=1-a[3];

c[1]=6/h[1]((y[2]-y[1])/h[1]-y'[1]);

c[j_]:=6((y[j+1]-y[j])/h[j]-(y[j]-y[j-1])/h[j-1])/(h[j-1]+h[j]);

-25-

学院： 研究生学院 专业： 机械工程 组号： 39 成绩：

c[4]=6/h[4-1](y'[4]-(y[4]-y[4-1])/h[4-1]);

A=Table[Switch[i-j,-1,b[j-1],0,2,1,a[j+1],_,0],{i,1,4},{j,1,4}]; MatrixForm[%] CC=Table[c[j],{j,1,4}]; MatrixForm[%] LinearSolve[A,CC]; MatrixForm[%];

M[j_]:=LinearSolve[A,CC][[j]] Table[M[j],{j,1,4}]

S[j_]:=M[j+1](x-x[j])^3/(6h[j])-M[j](x-x[j+1])^3/(6h[j])+ (y[j+1]-M[j+1]h[j]^2/6)(x-x[j])/h[j]- (y[j]-M[j]h[j]^2/6)(x-x[j+1])/h[j] Table[S[j],{j,1,3}]; Expand[%]; MatrixForm[%]

g1=Plot[%[[1]],{x,27.7,28}] g2=Plot[%%[[2]],{x,28,29}] g3=Plot[%%%[[3]],{x,29,30}]

g4=ListPlot[B,Prolog->AbsolutePointSize[15]] Show[g1,g2,g3,g4,Prolog->AbsolutePointSize[15]]

3.1.2 举例验证

例子：已知函数f?x?的数值表如下：

x f?x? f'?x?

2 3 1 -26-

4 7 6 13 -1 学院： 研究生学院 专业： 机械工程 组号： 39 成绩：

试求f?x?在[2,6]上的三次样条插值函数

解：这是第一类边界条件的问题 ，n?2，hi?h?2，

?2???1?0?由公式知 a2??0?1，

?02?20??M0??c0???????1??M1???c1?

????2???M2??c2??1?h011?,?1?

h0?h1226?y1?y0'?'?c0???y0??3(f[x,x]?y010)?3(2?1)?3 ?h0?h0?c2?6?'y2?y1?'??y??3(y?f[x1,x2])??12 22??h1?h1?c1?6f[x0,x1,x2]?3 2得方程组

?2M0?M1?3??0.5M0?2M1?0.5M2?1.5 ?M?2M??122?1解得 M0?0.25,M1?0.25,M2??7.25 故所求的三次养条差值函数

(x?x0)3(x?x1)3S0(x)?M1?M0?6h06h0?(??y0M0yM?h0)(x?x1)?(1?1h0)(x?x0)h06h06

15178(x?4)3?(x?2)3?(x?4)?(x?2),x?[2,4] 48241235298107(x?6)3?(x?4)3?(x?6)?(x?4),x?[4,6] 2448312-27-

同理有

S1(x)??即：

学院： 研究生学院 专业： 机械工程 组号： 39 成绩：

5178?133?(x?4)?(x?2)?(x?4)?(x?2),x?[2,4]??4824123S(x)??

??5(x?6)3?29(x?4)3?8(x?6)?107(x?4),x?[4,6]?48312?24Maths程序如下：

Clear[x,y,a,b,c,n,M] x[i_]:=2i; y[1]=3; y[2]=7; y[3]=13;

B=Table[{x[i],y[i]},{i,1,3}]; y'[1]=1; y'[3]=-1; h[j_]:=2;

a[j_]:=h[j-1]/(h[j-1]+h[j]); a[3]=1; b[1]=1; b[j_]:=1-a[j];

c[1]=6/h[1]((y[2]-y[1])/h[1]-y'[1]);

c[j_]:=6((y[j+1]-y[j])/h[j]-(y[j]-y[j-1])/h[j-1])/(h[j-1]+h[j]); c[3]=6/h[3-1](y'[3]-(y[3]-y[3-1])/h[3-1]);

A=Table[Switch[i-j,-1,b[j-1],0,2,1,a[j+1],_,0],{i,1,3},{j,1,3}]; MatrixForm[%] CC=Table[c[j],{j,1,3}]; MatrixForm[%] LinearSolve[A,CC]; MatrixForm[%];

M[j_]:=LinearSolve[A,CC][[j]] Table[M[j],{j,1,3}]

S[j_]:=M[j+1](x-x[j])^3/(6h[j])-M[j](x-x[j+1])^3/(6h[j])+

-28-