=di
=d?
a2x1
+b?
x?
+c?
X3
*
akXk」+bkXk+CkXk书
=dk
an」Xn_2+bn二Xn」+Cn」Xn
=dn」
a“Xn」+SXn
=dn
这种方程组的系数矩阵A为三对角矩阵，即
b1c1
a2
+
C2
A=
ak
bk
Ck
+
an_Lbn_lCn_i
anbnJ
追赶法实际上是高斯消去法的一种简化形式，它同样分消元与回代两个过
程。
求解基本步骤如下：
（1）计算的递推公式
★G/b—iy/g—a「）,i=2,3,|1（,n—1.
（2）解Ly=f
yi=fi/bi,%=（£—ayj/g—aj=2,3,|l（,n;
⑶解Ux=y
Xn=yn,Xi=%-jxi,i=n-1,n-2」H2,1.
将计算系数mTH-；_i及y」y2山-yn的过程称为追的过
程。
将计算方程组的解Xn—•Xn」Tl（'■Xi的过程称为追赶的过程。
4、Simpson公式
Simpson公式是一个数值积分公式，在计算一些多项式函数（三次或三次一
下）的定积分时会得出精确值。
但容易验证它对于f（x）=x4通常是不准确的，因
此，Simpson公式实际上具有三次代数精度，但是配合三次样条插值使用，三次代数精度已满足。
b
将定积分I二f（x）dx的积分区间［a,b］分隔为n等分，各节点为
a
ab、
）f（b）]
2
得至USimpson公式:
b—a
S[f（a）4f（
6
、Matlab程序
1、三次样条插值主程序：
clc;clear;
%%
globalxyhsn;%定义全局变量%
y=[00.511.5234];%船舶吃水深度%
x=[01.6532.0892.3882.6193.0493.654];%对应水深下的船宽%
y_0=0;y_n=0;%第一类边界条件%
n=length（x）;
fori=1:
n-1
h（i）=x（i+1）-x（i）;%计算分化间隔%
end
fori=2:
n-1
u（i-1）=h（i-1）/（h（i-1）+h（i））;%计算三对角的下对角线1Un-2%
v（i）=h（i）/（h（i-1）+h（i））;%计算三对角的上对角线2到n-1%d（i）=6*（（y（i+1）-y（i））/h（i）-（y（i）-y（i-1））/h（i-1））/（h（i-1）+h（i））;
%计算d向量（1到n-1）%
end
u（n-1）=1;v
（1）=1;%补充定义三对角元素%
d
（1）=6*（（y
（2）-y
（1））/h
（1）-y_0）/h
（1）;
d（n）=6*（y_n-（y（n）-y（n-1））/h（n-1））/h（n-1）;%补充定义d向量%
w（1:
n）=2;%定义三对角的中对角线1到n%
M=zhuigan（u,w,v,d,n）;%调用追赶法函数，n为方程组的阶数%
%%
fori=1:
n-1
%计算每个分化的三次多项式，conv%
a=y（i）-M（i）*h（i）A2/6;
b=y（i+1）-M（i+1）*h（i）A2/6;
s（i,:
）=（M（i）*conv（conv（[-1,x（i+1）],[-1,x（i+1）]）,[-1,x（i+1）]）+M（i+
1）*conv（conv（[1,-x（i）],[1,-x（i）]）,[1,-x（i）]））/（6*h（i））+（a*[0,0,-1,x（i+1）]+b*[0,0,1,-x（i）]）/h（i）;
end
%%
fori=1:
n-1
e=linspace（x（i）,x（i+1）,10）;
f=polyval（s（i,:
）,e）;%代入e计算多项式的值%
plot（e,f）;
holdon;
%根据以上多项式，做出图像，holdon%
end
%%
symsX;%定义符号变量%
fori=1:
n-1
p（i）=vpa（poly2sym（s（i,:
）,X）,4）;%表示出每个分段的多项式表达式
end
%%
fori=1:
36
x0（i）=0+i*0.1;%吃水从0到设计水线%
A（i）=simpson（x0（i））;%调用simpson函数%
end
plot（x0,A）;
%绘制对应站位上的邦戎曲线（横截面面积与吃水的关系）%
2、追赶法子程序：
%追赶发%
functionx=zhuigan（a,b,c,f,n）
bt
（1）=c
（1）/b
（1）;
fori=2:
n-1%求解3%
bt（i）=c（i）/（b（i）-a（i-1）*bt（i-1））;
end
y
（1）=f
（1）/b
（1）;
fori=2:
n%LUx=f,先计算Ly=f%
y（i）=（f（i）-a（i-1）*y（i-1））/（b（i）-a（i-1）*bt（i-1））;
end
x（n）=y（n）;
fori=n-1:
-1:
1%计算Ux=y%
x（i）=y（i）-bt（i）*x（i+1）;
end
end
3、Simpsor法子程序
%simpson法计算数值积分（配合三次样条曲线）%
functionA=simpson（x0）%返回面积A,输入任一水深x0%
globalxyhsn;%申明全局变量%
A=0;%面积初始化%
%%
fori=1:
n-1
if（x0-x（i））>0&&（x0-x（i+1））<0
k=i;
%判断x0所在的分片位置%
break;
end
end
fori=1:
k-1%mean_i为中间位置，从h⑴到h（n-1）叠加计算面积%
mean_i=polyval（s（i,:
）,x（i）/2+x（i+1）/2）;
A=A+2*h（i）/6*（y（i）+y（i+1）+4*mean_i）;
end
%%
y_x0=polyval（s（k,:
）,x0）;
step=x0-x（k）;
mean_k=polyval（s（k,:
）,x（k）/2+x0/2）;
A=A+2*step/6*（y（k）+4*mean_k+y_xO）;
%计算余下的面积，2*是因为船舶面积是对称%
end
三、计算结果
首先把型值表里的数据按照表1的格式输入到Excel中，以备计算程序调用，然后把与计算有关的数据用C#语言实时地调用到程序中，根据相关数据计算得出
船舶横剖面的三次样条插值函数，最后用辛普森积分法对该函数进行积分，得到
邦戎曲线。
表1某船部分型值表
曲号/水緩号
HI.
500UI.
100041”
15001*（.
皿加1
3000U1
4D00WI.
L5
94
109S
LW
1603
1&KI
174S
2U25
2.0
310
1653
2倔9
2斑
2619
3049
3654
2.5
1052
3171
4U50
4加
5218
6LS2
7029
计算流程图如下:
的浸水面积，绘制邦戎曲线。
图：
二号站位的横剖面曲线
对以上得到的横剖面半宽图利用Simpon积分得到邦戎曲线，如下图所示:
图：
二号站位的邦戎曲线
各分划的三次函数，表示如下：
[0，1.653]:
0.1058*XA3+0.008074*XA2-3.469*10A（-17）*X
[1.653,2.089]:
0.1072*XA3+0.0009915*乂八2+0.01171*X-0.006451
[2.089,2.388]：
0.5718*XA3-2.91*XA2+6.093*X-4.241
[2.388，2.619]：
12.62*XA2-1.596*XA3-31.0*X+25.28
[2.619，3.049]：
0.1255*XA3-0.906*XA2+4.432*X-5.647
[3.049，3.654]：
22.72*XA2-2.458*XA3-67.61*X+67.57
四、总结与展望
利用三次样条插值计算船舶中的邦戎曲线，具有一定的专业价值。
在后续的
工作中，可以计算费尔索夫曲线、船型稳性力臂插值曲线以及破损船舶的浮性和稳性。
将数值计算的知识应用到船舶与海洋工程专业中去，即提高对本专业的认
识和matlab编程的能力，同时更加深刻的体会到计算方法的神奇之处！
（wujun）
升级会员 