用三弯矩法求解三次样条插值函数

数值分析课程设计报告书院系名称：学生姓名：专业名称：班级：时间：实验一 三次样条插值的三弯矩法一、实验目的已知数据i x ,()i i y f x =,0,,i n =及边界条件()n j x y j j 1,0),(2=,求)(x f 的三次样条插值函数)(x S .要求输出用追赶法解出的弯矩向量0[,,]n M M M =及()(),0,,,0,1,2k i S t i m k ==的值.画出)(x S y =的图形,图形中描出插值点(,)i i x y 及(,())i i t S t 分别用‘o ’和‘*’标记.二、实验原理1.用追赶法求解第二类边界条件的三弯矩方程:0010012111121111[,,]21[,,]26[,,]212[,,]n n n n n n n n n n f x x x M f x x x M M f x x x M f x x x μλμλ------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 其中1111,,j jj j j j j j j j j h h h x x h h h h μλ-+--===-++.2.得出样条函数表达式:332211111()()()()()6666j j j j j j j j j j j j j j j jx x x x M h x x M h x x S x M M y y h h h h +++++----=++-+-. 3.计算(k)(),0,,,0,1,2i S t i m k ==.三、实验结果所用数据:x=[-2.223,-1.987,-1.8465,-1.292,-1.2266,-1.1056,-0.8662,-0.6594,-0.2671,-0.0452,0.5385,1.2564,1.4398,1.5415,1.7646,1.9678,2.236];y=[0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662];d2s1= -4.5000;d2sn= -4.8967; %第二种边界条件t=[-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236]; ;(指定计算点)计算结果:-2.5-2-1.5-1-0.500.51 1.52 2.50.811.21.41.61.82四、实验分析通过实验结果我们，知道三弯矩法求出满足初始条件的三次样条函数，与其他插值函数的构造相比，三次样条插值法的计算量要小得多。

摘要求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。

这时我们可以通过数值方法求出函数积分的近似值。

在用近似值代替真实值时，遇到的问题就是近似值的代数精度是否足够。

当代数精度不足够时，很显然提高插值函数的次数是一种方法，但是考虑到数值计算的稳定性，当次数过高时，会出现龙格现象，用增大n的方法来提高数值积代数精度是不可取的。

正如我们所知道的分段线性插值，逼近程度好，但光滑性差。

分段三次Hermite插值，逼近程度好，光滑性也有所提高，但是需要增加更多的条件，不太实用。

因此，我们将介绍一种结合二者的优点的插值方法——三次样条插值。

本实验将介绍三次样条插值的三弯矩算法。

关键词：龙格现象三弯矩算法代数精度分段三次Hermite插值1、实验目的1) 通过本次实验体会并学习三次样条插值的优点。

2) 通过对三次样条插值进行编程实现，提高自己的编程能力。

3) 用实验报告的形式展现，提高自己在写论文方面的能力。

2、算法流程如果已知函数)(x f y =在节点a =x 0<x 1<⋯<x n =b,y i =f (x i ),i =0,1,2,⋯,n 处的函数值和导数值：n i x f y i i ,,2,1,0),( ==如果)(x S 满足条件：1) S (x )是一个分段的三次多项式且i i y x s =)(；2) S (x )在[a,b]具有二阶连续导数。

则称S (x )是三次样条插值函数。

S (x )的具体形式为：其中S i (x)在[x n−1,x n ]上是三次多项式S i (x )=a i x 3+b i x 2+c i x +d i 由插值条件S (x i )=y i ，i=0,1,2,…,n ，得n+1个条件。

边界条件一：S ′(x 0)=y 0′,S ′(x n )=y n ′边界条件二：S ′′(x 0)=y 0′′,S ′′(x n )=y n ′′边界条件三：假定函数y =f(x)是以b-a 为周期的周期函数，这时要求S(x)也是周期函数，即{S (x 0+0)=S(x n −0)S ′(x 0+0)=S ′(x n −0)S ′′(x 0+0)=S ′′(x n +0)()()()()⎪⎪⎩⎪⎪⎨⎧∈∈∈=-]12,121,01,[,...............][,][,n n n x x x x s x x x x s x x x x s x s针对三种边界条件的求解方法的不同，可以分为三转角算法和三弯矩算法，本实验将介绍和学习三转角算法。

三转角样条插值摘要:对于插值多项式的次数问题，人们一般认为：插值多项式的次数越高，精度越高。

比如线性插值的误差就比抛物线插值的误差要大。

但是在20世纪初龙格（Runge ）现象说明插值多项式并非次数越高精度越高。

人们利用样条插值来拟合曲线，可以利用足够的的数据点，还有公式简单，运算量节省等好处，其中最常用的样条函数是三次样条。

在课本《现代数字数学与计算》中采用的是三弯矩方法推导得出三次样条插值多项式的计算公式，本文将采用三转角法推出三次样条插值多项式的计算公式，并用其拟合几个函数。

在发现等距节点在拟合不连续函数的结果不让人满意时，使用Chebyshev 插值节点构造不均匀节点来减少其拟合的误差，最后用三转角样条插值函数拟合画出手写字母。

三转角法的理论：对于b t t t a n =<<<=...10，考虑],[2b a C s ∈的三次样条插值函数，使得n i t f f t s i i i ,...,1,0),()(===，且满足第一类边界条件00'0')()(β==t f t s ，n n n t f t s β==)()(''。

验证在区间n i t t i i ,...,2,1],,[1=-上满足)()()()()(11t v t u t q f t p f t s i i i i i i i i i ββ+++=--其中1--=i i i t t h ，)](2[)()(132--+-=i i ii i t t h ht t t p ， )](2[)()(321i i ii i t t h h t t t q ---=- ，212)()()(ii i i h t t t t t u ---=， 221)()()(ii i i h t t t t t v --=- ，而 i β可由递推公式)(3)11(31)11(21212112121111ii i i i ii i i i i ii ih f h f h h f h h h h -++++++--+-=+++βββ所确定。

数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。

实验函数:dt ex f xt ⎰∞--=2221)(πx 0.0 0.1 0.2 0.3 0.4 F(x)0.50000.53980.57930.61790.7554求f(0.13)和f(0.36)的近似值实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。

实验名称 实验4.3三次样条插值函数（P126）4.5三次样条插值函数的收敛性（P127）实验时间姓名班级学号成绩实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。

对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。

实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。

实验要求：（1）随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2）三次样条插值函数的思想最早产生于工业部门。

作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：x0 1 2 3 4 5 6 7 8 9 10 ky0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 ky 0.8 0.2 k算法描述：拉格朗日插值：错误!未找到引用源。

其中错误!未找到引用源。

是拉格朗日基函数，其表达式为：()∏≠=--=nij j j i ji x x x x x l 0)()(牛顿插值：))...()(](,...,,[....))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i ji j i j i三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-1,x i ]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131i i ii i i i i i i i i i i i i i i i i i x x x h yM h M h h y x M M h h y y h x x Mi h x x M x S -------∈-+-+---+-+-=式中Mi=)(i x S ''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n nn nn ih y y mn h M M m h y y h M M对于第二类边界条件有⎩⎨⎧=+=+-nn n n d M M d M M 2210100μλ其中nn n n nn n M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ 对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j));end;%拉格朗日基函数f=f+l*Y(i);endfprintf('%d\n',f)return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi)%X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1));else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi) % X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*( xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i) )*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2); hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29];dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1); B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)-M(i))/(6* h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一张图上4.5.1增加插值节点观察误差变化从上面三张图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线。

/* 三次样条插值计算算法*/#include "math.h "#include "stdio.h "#include "stdlib.h "/*N：已知节点数N+1R：欲求插值点数R+1x，y为给定函数f(x)的节点值{x(i)} (x(i) <x(i+1)) ,以及相应的函数值{f(i)} 0 <=i <=NP0=f(x0)的二阶导数；Pn=f(xn)的二阶导数u:存插值点{u(i)} 0 <=i <=R求得的结果s(ui)放入s[R+1] 0 <=i <=R返回0表示成功，1表示失败*/int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[]){/*声明局部变量*/double *h; /*存放步长:{hi} 0 <=i <=N-1 */double *a; /*存放系数矩阵{ai} 1 <=i <=N ；分量0没有利用*/ double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ；*/double *ba; /*存放中间结果0 <=i <=N-1*/double *m; /*存放方程组的解{m(i)} 0 <=i <=N ；*/int i,k;double p1,p2,p3,p4;/*分配空间*/if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);/*第一步：计算方程组的系数*/for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++)c[k]=1-a[k];for(k=1;k <N;k++)g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]); c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;/*第二步：用追赶法解方程组求{m(i)} */ba[0]=c[0]/2;g[0]=g[0]/2;for(i=1;i <N;i++){af[i]=2-a[i]*ba[i-1];g[i]=(g[i]-a[i]*g[i-1])/af[i];ba[i]=c[i]/af[i];}af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N]; /*P110 公式：6.32*/ for(i=N-1;i> =0;i--)m[i]=g[i]-ba[i]*m[i+1];/*第三步：求值*/for(i=0;i <=R;i++){/*判断u(i)属于哪一个子区间，即确定k */if(u[i] <x[0] || u[i]> x[N]){/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 1;}k=0;while(u[i]> x[k+1])k++;//p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); //p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 0;}void main(){int N,R;double *x,*y,*u,*s;double P0,Pn;int i;/*验证算法：*/N=7;R=6;/*分配空间*/if(!(x=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(y=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(u=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(s=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9 463;u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;P0=-0.4794;Pn=-0.9463;if(!SPL( N, R, x, y, P0, Pn, u, s)){/*打印结果*/printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\ns= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf( "\nsin= ");for(i=0;i <=R;i++)printf( "%9.5f ",sin(u[i]));}/*释放空间*/free(x);free(y);free(u);free(s);}/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/ //输入327.7 4.128 4.329 4.130 3.013.0 -4.0//输出输出三次样条插值函数：1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)//三次样条插值函数#include<iostream>#include<iomanip>using namespace std;const int MAX = 50;float x[MAX], y[MAX], h[MAX];float c[MAX], a[MAX], fxym[MAX];float f(int x1, int x2, int x3){float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);return (a - b)/(x[x3] - x[x1]);} //求差分void cal_m(int n){ //用追赶法求解出弯矩向量M……float B[MAX];B[0] = c[0] / 2;for(int i = 1; i < n; i++)B[i] = c[i] / (2 - a[i]*B[i-1]);fxym[0] = fxym[0] / 2;for(i = 1; i <= n; i++)fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);for(i = n-1; i >= 0; i--)fxym[i] = fxym[i] - B[i]*fxym[i+1];}void printout(int n);int main(){int n,i; char ch;do{cout<<"Please put in the number of the dots:";cin>>n;for(i = 0; i <= n; i++){cout<<"Please put in X"<<i<<':';cin>>x[i]; //cout<<endl;cout<<"Please put in Y"<<i<<':';cin>>y[i]; //cout<<endl;}for(i = 0; i < n; i++) //求步长h[i] = x[i+1] - x[i];cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";int t;float f0, f1;cin>>t;switch(t){case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";cin>>f0>>f1;c[0] = 1; a[n] = 1;fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];break;case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";cin>>f0>>f1;c[0] = a[n] = 0;fxym[0] = 2*f0; fxym[n] = 2*f1;break;default:cout<<"不可用\n";//待定};//switchfor(i = 1; i < n; i++)fxym[i] = 6 * f(i-1, i, i+1);for(i = 1; i < n; i++){a[i] = h[i-1] / (h[i] + h[i-1]);c[i] = 1 - a[i];}a[n] = h[n-1] / (h[n-1] + h[n]);cal_m(n);cout<<"\n输出三次样条插值函数：\n";printout(n);cout<<"Do you to have anther try ? y/n :";cin>>ch;}while(ch == 'y' || ch == 'Y');return 0;}void printout(int n){cout<<setprecision(6);for(int i = 0; i < n; i++){cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";/*cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";cout<<endl;*/float t = fxym[i]/(6*h[i]);if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";else cout<<-t<<"*(x - "<<x[i+1]<<")^3";t = fxym[i+1]/(6*h[i]);if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";cout<<"\n\t";t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";cout<<endl<<endl;}cout<<endl;}。

非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。

%%%%%%%%%%%%%%%%%%%function [newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b)% 三转角样条插值% 将插值点分两次输入，x0 y0 单独输入% 边值条件a的一阶导数 y1a 和b的一阶导数 y1bn=length(x);m=length(y);if m~=nerror('x or y 输入有误，再来');endv=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);w=2*ones(n-1,1);h0=x(1)-x0;h=zeros(n-1,1);for k=1:n-1h(k)=x(k+1)-x(k);endv(1)=h0/(h0+h(1));u(1)=1-v(1);d(1)=3*(v(1)*(y(2)-y(1))/h(1)+u(1)*((y(1)-y0))/h0);%for k=2:n-1v(k)=h(k-1)/(h(k-1)+h(k));u(k)=1-v(k);d(k)=3*(v(k)*(y(k+1)-y(k))/h(k)+u(k)*(y(k)-y(k-1))/h(k-1));endd(1)=d(1)-u(1)*y1a;d(n-1)=d(n-1)-v(n-1)*y1b;newv=v(1:n-2,:);newu=u(2:n-1,:);%%%%%%%%%%%%function intersanzhj(x,y,x0,y0,y1a,y1b)% 三转角样条插值%第一部分n=length(x);m=length(y);if m~=nerror('x or y 输入有误，再来');end%重新定义hh=zeros(n,1);h(1)=x(1)-x0;for k=2:nh(k)=x(k)-x(k-1);end% 调用三转角函数[a,b,c,d]=sanzhj(x,y,x0,y0,y1a,y1b);% 三对角方程m=chase(a,b,c,d);% 求MM=[1;m;0];% 求插值函数fprintf('三次样条（三转角）插值的函数表达式\n');syms X ;fprintf('S0--1:\n');S(1)=collect((y0/h(1).^3)*(X-x(1)).^2*(h(1)+2*(X-x0))+(y(1)/h(1).^3)*(X-x0).^2*(h(1)+2*(x(1)-X))+(M(1)/h(1).^2)*(X-x0)*(X-x(1)).^2+(M(2)/h(1).^2)*(X-x(1))*(X-x0).^2);for k=2:nfprintf('S%d--%d:\n',k-1,k);S(k)=collect((y(k-1)/h(k).^3)*(X-x(k)).^2*(h(k)+2*(X-x(k-1)))+(y(k)/h(k).^3)*(X-x(k-1)).^2*(h(k)+2*(x(k)-X))+(M(k)/h(k).^2)*(X-x(k-1))*(X-x(k)).^2+(M(k+1)/h(k).^2)*(X-x(k))*(X-x(k-1)).^2);endS=S.';disp(S);fprintf('以上为样条函数（三转角）解析式，显示为手写习惯如下：\n');pretty(S);%第二部分%是否继续运行程序myloop=input('继续运行程序输入“1”，否则输入“0”\n');if myloopwhile myloopxi=input('输入需要计算的点的值，并按回车键\n');if xi>x0|xi<x(n)fprintf('现在开始计算输入点的插值函数值……\n');elsefprintf('输入数值不在插值范围内，请重新输入\n');xi=input('输入需要计算的点的值，并按回车键……\n'); end% 确定输入的数值应该使用哪个解析式newx=[x0;x];[r,suoy]=min(abs(newx-xi));fprintf('输入点的插值函数值为:\n\n');fprintf('\t');if xi<=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input('继续计算输入“1”，终止计算输入“0”\n'); endelsereturn;end。

20 年月日实验标准题已知如下数据：且端点约束条件为f’(-1)=5，f’(3.50)=29.16，用三次样条插值的三弯矩法求f （-0.02）和f （2.56）。

一、理论依据：三次样条函数点a=X 0<X 1<X 2<...<X n-1<X n =b 将区间[a,b]分成n 个小区间，若函数S(x)满足:(1)在每个子区间[x i ，x i+1]（i=0,1,2,3,n-1)上S(x)是三次多项式 (2) S(x)∈C 2[a,b]则称S(x)是区间[a,b]上的三次样条函数。

求f(x)在[a,b]上的三次样条函数，可设S(x)=a i x 3+b i x 2+c i x+d i ,x ∈[x i ,X i+1](i=0,1,2,...n-1)其中，a i ,b i ,C i ,d;为待定常数。

二、操作原理：依据给定的x 、f （x ）数据点计算下面各项的值：()()[]()11111i 111i ,,6,,6)1,,2,1(1+--+++++=-+=-⋅⋅⋅=+=-=+=i i i i i i i ii ii i i i ii ix x x f x x f x x f h h d n i h h h h h h μλμ()()⎥⎦⎤⎢⎣⎡-'=⎥⎦⎤⎢⎣⎡'-=-n n n nx x f f h d f x x f h d ,6,61n 01010⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡----n n n d d d M M M M 1n 101n 1011n 11d 212212λμλμ三弯矩方程组可以统一写成AM=d 的形式，系数矩阵A 为三对角矩阵，可以使用追赶法求解M 。

将所求的M 带入下面的式子计算出三次样条插值函数：121111121131131i )6()6(6)(6)()(+++++++++++--+--+-+-=i i i i j i i i i i i i i i i h x x h M f h x x h M f h x x M h x x M x S 三、实验代码： 主函数：b=[-1. -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.5;-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31]; fa=5;fb=29.16; c=[-0.02;2.54]; [sx] = scyt(b,fa,fb,c) 调用函数：LU 求解三对角方程： function [z]=lux(a,b)% lux 追赶法求解三对角函数% a 系数矩阵 b 列向量，z 为求解结果 [l,u]=lu(a); [n,~]=size(b); z=zeros(n,1);z(1)=b(1); for i=2:nz(i)=b(i)-l(i,i -1)*z(i -1); endz(n)=z(n)/u(n,n);for j=(n -1):-1:1 %计算x 值 z(j)=(z(j)-u(j,j+1)*z(j+1))/u(j,j); end end三次样条函数插值函数： function [sx] = scyt(b,fa,fb,c)% scyt三次样条函数插值% b已知数据2*1 fa为下端约束条件fb为上端约束条件c为需要计算的列% sx为计算出的某点插值函数值b=b';[f,y]=size(b);for i=2:f-1b(i,y+1)=b(i,1)-b(i-1,1); %h1b(i+1,y+1)=b(i+1,1)-b(i,1); %h2b(i,y+2)=b(i,y+1)/(b(i,y+1)+b(i+1,y+1)); %计算出μb(i,y+3)=6*((b(i+1,2)-b(i,2))/b(i+1,y+1)-(b(i,2)-b(i-1,2))/b(i,y+1))/(b(i,y+1)+b(i+1,y+1)); %计算出dend[f,y]=size(b);b(1,y)=6/b(2,3)*((b(2,2)-b(1,2))/b(2,3)-fa); %补齐d0，dnb(f,y)=6/b(f,3)*(fb-(b(f,2)-b(f-1,2))/b(f,3));a=2*eye(f);a(1,2)=1;a(f,f-1)=1;s=1; %系数矩阵建立afor j=2:f-1a(j,s)=b(j,4);a(j,s+2)=1-b(j,4);s=s+1;endd=b(:,y); %非齐次项常数项dm=lux(a,d); %求解待定常数m弯矩值[t,~]=size(c);st=0;for g=1:t %判断计算数字所用方程式并计算for h=1:fif c(g,1)>b(h,1)st=st+1;endendif st==fst=st-1;endif st==0st=1;endk=st+1;x=c(g,1);st=0;x1=b(k-1,1);x2=b(k,1);h=b(k,3);m1=m(k-1,1);m2=m(k,1);f1=b(k-1,2);f2=b(k,2);sx(g,1)=(m1*(x2-x)^3+m2*(x-x1)^3)/(6*h)+(f1-m1*h^2/6)*(x2-x)/h+(f2-m2*h^2/6)*( x-x1)/h;endend四、求解结果：sx =[-3.1058;4.7834]所以求f（-0.02）=-3.1058、f（2.56）=4.7834用Romberg 算法求以下积分，允许误差eps=0.00001，()dx x x x x 24.131sin 753+⎰。

三次样条插值0 引⾔三次样条插值以构造简单，使⽤⽅便，拟合准确，具有“保凸”的重要性质等特点成为了常⽤的插值⽅法。

⼀般三次样条插值解算过程中通过追赶法求解三弯矩阵，但使⽤计算机求解时会表现出解的精度不⾼的问题，导致其计算结果⽆法应⽤到⼯程实践之中。

因此需要找出⼀种提⾼解精度的⽅法。

1 基本概念三次样条函数的定义：在区间内对于给定的函数值，其中，如果函数满⾜条件：（1）在每个⼦区间，上都是不⾼于三次的多项式；（2）、、在上都连续；（3），。

则称为函数关于节点的三次样条函数。

想要求解三次样条插值函数，只需在每个⼦区间上确定⼀个三次多项式共有4个系数，确定它们需要 4n 个条件，因此要完全确定共需 4n 个条件。

由所满⾜的条件（1）、（2）、（3），可确定个条件，仍然缺少两个条件。

这两个条件通常由实际问题对三次样条插值函数在端点的状态要求给出，也称之为边界条件，常见的边界条件有：1）夹持边界条件（Clamped Spline）：给定两端点的⼀阶导数值，即，；2）⾃然边界条件（Natural Spline）：使两端点的⼆阶导数值为零，即；3）⾮扭结边界条件（Not-A-Knot Spline）：强制第⼀个插值点的三阶导数值等于第⼆个点的三阶导数值，最后第⼀个点的三阶导数值等于最后第⼆个点的三阶导数值，即，。

2 计算⽅法设三次样条函数，（0），，，由三次样条函数定义（1）（2）（3）可得：，（1）如下构造式（1）矩阵：（2）由式（1）可知：，，，，（3）1）在夹持边界条件时，，，，；，，，；2）在⾃然边界条件时，，，，；，，，；3）在⾮扭结边界条件时，，，，；，，，；由n个未知数的⾮齐次⽅程组有惟⼀解的充分必要条件是，可知矩阵⽅程（2）在以上三种情况下都有惟⼀解。

对矩阵⽅程（2）采⽤⾼斯列主元消去法即可求解得出。

最后，代⼊式（0）可以得出：，，，，3 应⽤算例有点集，在⾮扭结边界条件下进⾏插值。

同时使⽤Matlab R2010a和⽂章所述⽅法进⾏插值计算，对⽐计算结果。

用弯矩分量求三次样条插值函数的弯矩刘远鹏;刘光好;金少华【摘要】根据参考弯矩图解法原理,将相邻节点函数差和节点函数的导数变换为侧移分量和转角分量.利用各区间转角分量、侧移分量与函数的二阶导数形成弯矩-位移方程,以及连续梁支座相邻两跨转角方程,结合连续梁弯矩图的性质分析样条插值函数,做参考弯矩图;由此得到三次样条插值函数的支座弯矩.通过算例说明未知量参考弯矩图、荷载参考弯矩图、转角分量图和侧移分量图的使用方法.【期刊名称】《河北工业大学学报》【年(卷),期】2016(045)004【总页数】6页(P93-98)【关键词】特征点;转角分量;侧移分量;参考弯矩图;三次样条插值函数【作者】刘远鹏;刘光好;金少华【作者单位】河北工业大学工程训练中心,天津300401;河北工业大学土木工程学院,天津300401;河北工业大学理学院,天津300401【正文语种】中文【中图分类】TU311.4三次样条插值函数有广泛应用[1-3]，且有各种表达方法[4-5]，计算方法也多种．文献 [4]则采用构造三角阵的逆阵方法简化节点等距问题．通常可以将三次样条插值函数比拟为连续梁的挠曲线；其二阶导数视为连续梁的弯矩．因此三次样条函数的二阶导数分布就相当于无荷载作用有支座沉陷的抗弯刚度EI=1（本文采用无量纲）的连续梁弯矩图．用力学计算这种连续梁方法很多，可以利用广义图乘法及其弯矩分量[4-7]分析连续梁的原理．利用连续梁弯矩图和单跨梁弯矩图（分别是结点位移弯矩图和杆端位移弯矩图．在不致混淆时的也简称为弯矩图）的性质；求弯矩图从而获得结点弯矩，从而产生插值函数．整体计算过程简单，力学意义清楚．其计算量也得到优化．首先分析节点区间内的插值函数（这相当于单跨梁的挠曲线）．如下常见表达式：根据结构力学常用方式，对三次样条函数的连续梁称弦转角为：设连续梁第k跨梁的跨度l，据EI=1，则线刚度，侧移分量：左右两端的转角分量[5-6]统称这3个量为梁弯矩的位移分量，即弯矩分量．称梁内1/3点和2/3点分别为左特征点和右特征点[5-6]；设为L点和R点．利用式 (1)计算特征点的弯矩，与位移分量比较：这2个方程是梁的弯矩-位移方程，是梁的普遍形式的简化[5]．因连续梁相邻两跨的杆端转角分量表达式不同，但是杆端连接有相同转角，则：这个方程为连续梁的结点转角方程[7]．根据样条函数的二阶导数，相邻两跨的特征点弯矩Mk1+1/3、Mk1+2/3与杆端的弯矩Mk1、Mk有如下的关系：由此可知，弯矩分量有2种，它们是特征点的弯矩分解．若梁的弦转角和杆端转角为正方向，那么弯矩分量可用图1方法表示．连续梁弯矩图的性质：1）区间内2个特征点的侧移分量数值相同，分别在杆的上下两侧，与弦转角方向一致．如图1a）．2）转角分量的方向分别与其所表示的转角方向一致．如图1b）．3）在每个区间内，弯矩图是直线，如图1c）．4）特征点的弯矩及其分量满足梁的弯矩-位移方程和连续梁结点转角方程．连续梁相邻两跨的转角分，如图1d）．只要区间内插值函数是三次多项式；无论采用那种表达式等，上述性质都成立．而分析样条函数过程中，对结点转角的估计相当于对特征点的转角分量的估计．2.1 参考弯矩图分解计算过程中连续梁弯矩图是用已知量和未知量共同表示的，则为区别称其为参考弯矩图（用M表示）．那么据线性关系，将参考弯矩图M分解为[6]：其中：为荷载参考弯矩图；为未知量参考弯矩图．荷载参考弯矩图是利用连续梁的支座位移和一个已知条件（一端的已知弯矩或转角）计算形成；且表达式不含未知量．而未知量参考弯矩图根据参考弯矩图的未知量计算形成的弯矩图；不含已知条件（位移条件或杆端弯矩条件）．其各弯矩表达式是x的一次齐次式．2个弯矩图分别有各自的转角分量和侧移分量，分别用符号和表示．荷载参考弯矩图和未知量参考弯矩图都是杆内无荷载的连续弯矩图，则在每个弯矩图内分别具有连续梁弯矩图的上述4条性质．对于具体的计算过程来说，的值由已知的边界条件和0到n的xk，yk的值确定；这样对于来说yk的值全为0，则也全为0，带入式 (6)、式 (7)后得，，根据公式后得．例如对于第1种边界条件的参考弯矩图；根据其未知量，有荷载参考弯矩图的和；而对应的未知量参考弯矩图有．根据三次样条插值函数连续梁弯矩图的上述性质，通过作参考弯矩图方法求支座弯矩．将支座弯矩代入公式 (1)就生成区间内三次插值样条函数．2.2 弯矩分量图与算例已知函数值分布如表1所示[8]．边界条件：，求三次样条插值函数：所以跨度：．根据支座位移计算侧移分量：将它们标在连续梁上形成侧移分量图，如图2a）．利用式 (4)和式 (5)计算已知转角的转角分量：．将它们标在连续梁特征点上，如图2b）．同时设左端弯矩为未知量x．提取图2a）和图2b）的位移分量和转角分量，根据式(6)计算第1跨的L点弯矩：作荷载参考弯矩图．将1跨 L点弯矩标在连续梁上，如图3a）．因最左端弯矩等于0，根据式 (9)、式 10，作第1跨参考弯矩图直线，产生R特征点和右端的弯矩，如图3a）．根据式 (7)转角分量将其标在转角分量图上，如图2b）．利用结点转角方程计算第2跨左端转角分量，如图2b）．同理依次分析第2跨和第3跨的杆端弯矩和特征点的弯矩；依次产生各跨参考弯矩图直线．同时得到转角分量的分布．如图2b）（称其为转角分量图）．因此荷载参考弯矩图如图3a）所示．根据未知量假设作未知量参考弯矩图．这个图的A端转角等于0，且无侧移；则各个特征点弯矩等于转角分量．利用A端弯矩等于未知量x，A端转角分量等于0作第1跨弯矩图，如图3b）．利用特征点未知量参考弯矩与杆长成反比关系，则依次作第2跨、第3跨和第4跨的未知量参考弯矩图，如图3b）．荷载参考弯矩图与未知量参考弯矩图相加形成参考弯矩图．因第4跨的右转角分量和侧移分量是已知量，其和是R特征点弯矩；与参考弯矩图的特征点弯矩比较表达式建立方程：因此x=-2.028 6，计算各杆端的弯矩值，用式 (11)计算各杆端的支座弯矩：最后结果与其它方法一致．具体计算过程中，由于各个分量均有系数2，可以先将系数2提出，最后得到的数值乘以2以获得各段点弯矩．这样可以进一步简化计算过程．如果左端边界条件与此例相同，那么所需未知量以及建立的参考弯矩图的形式可不变；而依据右端边界条件建立不同的方程．对于第2种边界条件和第3种边界条件也可以用类似的方法求解．2.3 算例的MATLAB程序S为边界条件，com_la为侧移分量，com_l为左端转角分量，com_r为右端转角分量，MF为荷载弯矩（l为左端，r为右端），Mx为未知量弯矩（l为左端，r为右端），M为每个结点的弯矩．本文利用连续梁参考弯矩图解法形成分析计算三次插值样条函数的连续梁结点弯矩方法．其中利用参考弯矩图的两个层次分解：荷载未知量图分解，以及弯矩位移分量图分解．利用了梁的特征点的弯矩的侧移分量、转角分量，以及弯矩位移方程和连续梁结点转角方程．通过未知量参考弯矩图、荷载参考弯矩图、转角分量图和侧移分量图分析产生未知量方程，进一步得到三次插值样条函数的连续梁结点弯矩．计算过程仅需一个未知量，具有明显的并行计算特点．【相关文献】[1]陈盛，赵东标，陆永华，等．基于二维轮廓点云的螺纹中径计算 [J]．光学精密工程，2015，23（6）：1791-1799．[2]常景娜，高慧斌，乔冠宇．一种基于传感器的测量目标运动参数方法研究[J]．传感器与微系统，2015，34（11）：19-22．[3]梁栋，刘良栋，何英姿．月球精确软着陆最优标称轨迹在轨制导方法 [J]．中国空间科学技术，2011（6）：27-35．[4]Moon BS．Anexplicitsolution for thecubic spline interpolation for functionofasinglevariable[J]．Applied Mathematicsand Computation，2001，117（2-3）：251-255．[5]刘光好，赵光鉴．刚架参考弯矩图解法的侧移问题 [C]//崔京浩．第十一届全国结构工程学术会议论文集（VolⅠ）．长沙，2002：395-398．[6]刘光好，张扬，王俊杰，等．连续梁支座反力、弯矩、位移及三支座反力方程 [J]．公路交通科技（应用技术版），2007，31（3）：106-108．[7]Liu Guanghao，Yang Shaobin，LiuYuanpeng．Three-columniationendmomentequationof framewithsupportsinkageandangle[C]//2011INternationA l Conferenceon Remote Sensing，Environmentand Transportation Engineering（Vol5）．Nanjing,2011：4670-4682．[8]朱晓临．数值分析 [M]．合肥：中国科学技术大学出版社，2010：148-150．。