1、三次样条插值在船舶邦戎曲线中的应用船建学院 B1301095 wj1、计算原理1、三次样条插值原理三次样条插值多项式 是一种分段函数，它的应用范围很广，本文探讨)(xSn该方法在船舶静力学曲线计算和绘制中的应用。节点 ix分成的每个小区间 上是 3 次多项式，其在011( )naxxb1,i此区间上的表达式如下： 2233 111111()()()()()6 66,.iiiiiiiii ihxhxSxxMxyMyh n，因此，只要确定了 的值，就确定了整个表达式， 的计算方法如下：i i令： 11116()6(,)iii i iiiii iiihhyydfx，则 满足如下 n-1 个方程：iM
2、112,21iiiidn，对于第一种边界条件下有0110 11 ),(62,hfxMxfnnnn如果令 那么,)(6,1,),(6,1 1010 nnhxfdhfxfd解就可以写为 nnndM101011022 2、船舶静力学中的邦戎曲线船舶邦戎曲线是由一组船舶横剖面的面积曲线组成的，其中每条曲线表示该处横剖面在不同水线以下浸入水中的面积。邦戎曲线是船体纵向积分的基础，利用它可以计算船舶在不同吃水下的排水体积和浮心位置，进而为船舶的稳性与强度计算提供基本数据。因而邦戎曲线的精确性直接影响到船舶的安全性。传统邦戎曲线的计算以型值表为基础，利用梯形积分法，把船体某一横截面上各水线之间的面积近似成梯
3、形，然后把这些小梯形的面积求和得到，但梯形法只有一阶代数精度，对稳性计算要求较高的液货船来说，似乎不够精确。本文提出先把某一横截面上各水线面之间的水线面半宽用 1 个一元三次函数近似表达，然后在该水线面之间对这个一元三次函数进行积分，得到各水线面之间的小图形面积，再对这些小图形面积进行叠加，得到某一横截面的面积。3、追赶法计算对角占优方程组在实际问题中,经常遇到以下形式的方程组 nnnn kkkk dxbacxbadxcxb 1121123212 1这种方程组的系数矩阵 A 为三对角矩阵,即 nnkk baccba11221追赶法实际上是高斯消去法的一种简化形式,它同样分消元与回代两个过程。求
4、解基本步骤如下：(1) 计算 的递推公式i, .1/cb/(),23,1iiian(2) 解 Ly=f, 11/yf 1()/(),;iiiiiyfb(3) 解 Ux=y.1,2,niiixxn将计算系数 及 的过程称为追的过程。12n nyy将计算方程组的解 的过程称为追赶的过程。11nxx4、Simpson 公式Simpson 公式是一个数值积分公式，在计算一些多项式函数（三次或三次一下）的定积分时会得出精确值。但容易验证它对于 = 通常是不准确的，)(xf4因此，Simpson 公式实际上具有三次代数精度，但是配合三次样条插值使用，三次代数精度已满足。将定积分 的积分区间 分隔为 等分，
5、各节点为badxfI)(,bannjhaxj ,10,abh得到 Simpson 公式：)(2(4)6bfffS2、Matlab 程序1、三次样条插值主程序：clc;clear;%global x y h s n; %定义全局变量%y=0 0.5 1 1.5 2 3 4;%船舶吃水深度%x=0 1.653 2.089 2.388 2.619 3.049 3.654;%对应水深下的船宽%y_0=0;y_n=0;%第一类边界条件%n=length(x);for i=1:n-1h(i)=x(i+1)-x(i); %计算分化间隔%endfor i=2:n-1u(i-1)=h(i-1)/(h(i-1)+
6、h(i); %计算三对角的下对角线1到n-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）%endu(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);%调用追
7、赶法函数，n为方程组的阶数%for i=1:n-1%计 算每个分化的三次多项式，conv%a=y(i)-M(i)*h(i)2/6;b=y(i+1)-M(i+1)*h(i)2/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%for i=1:n-1e=linspace(x(i),x(i+1),10);f=polyval(s(i,:),e);%代入e计算多项式
8、的值%plot(e,f);hold on;%根据以上多项式，做出图像，hold on%end%syms X; %定义符号变量%for i=1:n-1p(i)=vpa(poly2sym(s(i,:),X),4); %表示出每个分段的多项式表达式%end%for i=1:36x0(i)=0+i*0.1;%吃水从0到设计水线%A(i)=simpson(x0(i);%调用simpson函数%endplot(x0,A);%绘制对应 站位上的邦戎曲 线（横截面面积与吃水的关系）%2、追赶法子程序：%追赶发 %function x=zhuigan(a,b,c,f,n)bt(1)=c(1)/b(1);for
9、i=2:n-1 %求解%bt(i)=c(i)/(b(i)-a(i-1)*bt(i-1);endy(1)=f(1)/b(1);for i=2:n %LUx=f,先计算Ly=f%y(i)=(f(i)-a(i-1)*y(i-1)/(b(i)-a(i-1)*bt(i-1);endx(n)=y(n);for i=n-1:-1:1%计算Ux=y%x(i)=y(i)-bt(i)*x(i+1);endend3、Simpson法子程序%simpson法计 算数值积分（配合三次 样条曲线）%function A=simpson(x0)%返回面积A，输入任一水深x0%global x y h s n; %申明全局变
10、量%A=0;%面 积初始化 %for i=1:n-1if (x0-x(i)0 %判断x0 所在的分片位置%break;endendfor i=1:k-1 %mean_i为中间位置，从h(1) 到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_x0);%计算余下
11、的面积，2*是因为船舶面积是对称%end3、计算结果首先把型值表里的数据按照表 1 的格式输入到 Excel 中，以备计算程序调用，然后把与计算有关的数据用 C#语言实时地调用到程序中，根据相关数据计算得出船舶横剖面的三次样条插值函数，最后用辛普森积分法对该函数进行积分，得到邦戎曲线。计算流程图如下：以某 14300DWT 散货船为例（下图） ，计算第 2 站的横剖面曲线去，并计算吃水从 03.6m 的浸水面积，绘制邦戎曲线。基于三次样条插值计算得到二号站位处的横剖面曲线图，如下图所示：图：二号站位的横剖面曲线对以上得到的横剖面半宽图利用 Simpon 积分得到邦戎曲线，如下图所示：图：二号站
12、位的邦戎曲线各分划的三次函数，表示如下：0，1.653： 0.1058*X3 + 0.008074*X2 - 3.469*10(-17)*X 1.653， 2.089：0.1072*X3 + 0.0009915*X2 + 0.01171*X - 0.0064512.089， 2.388：0.5718*X3 - 2.91*X2 + 6.093*X - 4.2412.388， 2.619：12.62*X2 - 1.596*X3 - 31.0*X + 25.282.619， 3.049：0.1255*X3 - 0.906*X2 + 4.432*X - 5.6473.049， 3.654：22.72*X2 - 2.458*X3 - 67.61*X + 67.574、总结与展望利用三次样条插值计算船舶中的邦戎曲线，具有一定的专业价值。在后续的工作中，可以计算费尔索夫曲线、船型稳性力臂插值曲线以及破损船舶的浮性和稳性。将数值计算的知识应用到船舶与海洋工程专业中去，即提高对本专业的认识和 matlab 编程的能力，同时更加深刻的体会到计算方法的神奇之处！（wujun）
