6

我正在尝试编写三次样条插值程序。我已经编写了程序，但是图表没有正确显示。样条曲线使用自然边界条件（开始/结束节点的第二个导数为 0）。代码在 Matlab 中，如下所示，

clear all
%Function to Interpolate
k = 10;                    %Number of Support Nodes-1
xs(1) = -1;
for j = 1:k
    xs(j+1) = -1 +2*j/k;   %Support Nodes(Equidistant)
end;
fs = 1./(25.*xs.^2+1);     %Support Ordinates
x = [-0.99:2/(2*k):0.99];  %Places to Evaluate Function
fx = 1./(25.*x.^2+1);      %Function Evaluated at x

%Cubic Spline Code(Coefficients to Calculate 2nd Derivatives)

f(1) = 2*(xs(3)-xs(1));
g(1) = xs(3)-xs(2);
r(1) = (6/(xs(3)-xs(2)))*(fs(3)-fs(2)) + (6/(xs(2)-xs(1)))*(fs(1)-fs(2));
e(1) = 0;

for i = 2:k-2
    e(i) = xs(i+1)-xs(i);
    f(i) = 2*(xs(i+2)-xs(i));
    g(i) = xs(i+2)-xs(i+1);
    r(i) = (6/(xs(i+2)-xs(i+1)))*(fs(i+2)-fs(i+1)) + ...
           (6/(xs(i+1)-xs(i)))*(fs(i)-fs(i+1));
end
e(k-1) = xs(k)-xs(k-1);
f(k-1) = 2*(xs(k+1)-xs(k-1));
r(k-1) = (6/(xs(k+1)-xs(k)))*(fs(k+1)-fs(k)) + ...
         (6/(xs(k)-xs(k-1)))*(fs(k-1)-fs(k));

%Tridiagonal System

i = 1;
A = zeros(k-1,k-1);
while i < size(A)+1;
    A(i,i) = f(i);
    if i < size(A);
        A(i,i+1) = g(i);
        A(i+1,i) = e(i);
    end
    i = i+1;
end

for i = 2:k-1                             %Decomposition
    e(i) = e(i)/f(i-1);
    f(i) = f(i)-e(i)*g(i-1);
end

for i = 2:k-1                             %Forward Substitution 
    r(i) = r(i)-e(i)*r(i-1);
end

xn(k-1)= r(k-1)/f(k-1);
for i = k-2:-1:1                          %Back Substitution
    xn(i) = (r(i)-g(i)*xn(i+1))/f(i);
end

%Interpolation

if (max(xs) <= max(x))
    error('Outside Range'); 
end

if (min(xs) >= min(x))
    error('Outside Range'); 
end


P = zeros(size(length(x),length(x)));
i = 1;
for Counter = 1:length(x)
    for j = 1:k-1
        a(j) = x(Counter)- xs(j);
    end
    i = find(a == min(a(a>=0)));
    if i == 1
        c1 = 0;
        c2 = xn(1)/6/(xs(2)-xs(1));
        c3 = fs(1)/(xs(2)-xs(1));
        c4 = fs(2)/(xs(2)-xs(1))-xn(1)*(xs(2)-xs(1))/6;
        t1 = c1*(xs(2)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(1))^3;
        t3 = c3*(xs(2)-x(Counter));
        t4 = c4*(x(Counter)-xs(1));
        P(Counter) = t1 +t2 +t3 +t4;
    else
        if i < k-1
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = xn(i+1)/6/(xs(i+1)-xs(i-1+1));
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1))-xn(i+1)*(xs(i+1)-xs(i-1+1))/6;
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        else
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = 0;
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1));
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        end    
    end
end

P = P';
P(length(x)) = NaN;

plot(x,P,x,fx)

当我运行代码时，插值函数不是对称的，并且不能正确收敛。任何人都可以就我的代码中的问题提供任何建议吗？谢谢。

4

2 回答 2

9

我很久以前在 Mathematica 中写了一个三次样条包。这是我将该包翻译成 Matlab 的过程。注意我已经有 7 年没有看过三次样条了，所以我是根据我自己的文档来做的。你应该检查我所说的一切。

基本问题是给定n数据点(x(1), y(1)) , ... , (x(n), y(n))，我们希望计算分段三次插值。插值定义为

   S(x) = {  Sk(x)   when x(k) <= x <= x(k+1)
          {  0       otherwise

这里 Sk(x) 是形式为的三次多项式

  Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3

样条曲线的属性是：

  1. 样条通过数据点Sk(x(k)) = y(k)
  2. 样条在端点处是连续的，因此在插值区间的任何地方都是连续的Sk(x(k+1)) = Sk+1(x(k+1))
  3. 样条具有连续的一阶导数Sk'(x(k+1)) = Sk+1'(x(k+1))
  4. 样条具有连续的二阶导数Sk''(x(k+1)) = Sk+1''(x(k+1))

要从一组数据点构造三次样条曲线，我们需要求解系数 、 和sk0每个sk1三次多项式。那是一个未知数。属性 1 提供约束，属性 2、3、4 各提供一个附加约束。因此，我们有约束和未知数。这留下了两个自由度。我们通过在开始和结束节点处设置二阶导数为零来固定这些自由度。sk2sk3n-14*(n-1) = 4*n - 4nn-2n + 3*(n-2) = 4*n - 64*n - 4

让和。m(k) = Sk''(x(k))_ 以下三项递推关系成立h(k) = x(k+1) - x(k)d(k) = (y(k+1) - y(k))/h(k)

  h(k-1)*m(k-1) + 2*(h(k-1) + h(k))*m(k) + h(k)*m(k+1) = 6*(d(k) - d(k-1))

m(k) 是我们希望求解的未知数。和由输入数据定义h(k)。d(k)这个三项递归关系定义了一个三对角线性系统。一旦m(k)确定 的系数Sk由下式给出

   sk0 = y(k)
   sk1 = d(k) - h(k)*(2*m(k) + m(k-1))/6
   sk2 = m(k)/2
   sk3 = m(k+1) - m(k)/(6*h(k))

好的，这就是您完全定义计算三次样条的算法所需要知道的所有数学知识。这是在Matlab中：

function [s0,s1,s2,s3]=cubic_spline(x,y)

if any(size(x) ~= size(y)) || size(x,2) ~= 1
   error('inputs x and y must be column vectors of equal length');
end

n = length(x)

h = x(2:n) - x(1:n-1);
d = (y(2:n) - y(1:n-1))./h;

lower = h(1:end-1);
main  = 2*(h(1:end-1) + h(2:end));
upper = h(2:end);

T = spdiags([lower main upper], [-1 0 1], n-2, n-2);
rhs = 6*(d(2:end)-d(1:end-1));

m = T\rhs;

% Use natural boundary conditions where second derivative
% is zero at the endpoints

m = [ 0; m; 0];

s0 = y;
s1 = d - h.*(2*m(1:end-1) + m(2:end))/6;
s2 = m/2;
s3 =(m(2:end)-m(1:end-1))./(6*h);

这是一些绘制三次样条的代码：

function plot_cubic_spline(x,s0,s1,s2,s3)

n = length(x);

inner_points = 20;

for i=1:n-1
   xx = linspace(x(i),x(i+1),inner_points);
   xi = repmat(x(i),1,inner_points);
   yy = s0(i) + s1(i)*(xx-xi) + ... 
     s2(i)*(xx-xi).^2 + s3(i)*(xx - xi).^3;
   plot(xx,yy,'b')
   plot(x(i),0,'r');
end

这是一个构造三次样条曲线并绘制在著名的龙格函数上的函数：

function cubic_driver(num_points)

runge = @(x) 1./(1+ 25*x.^2);

x = linspace(-1,1,num_points);
y = runge(x);

[s0,s1,s2,s3] = cubic_spline(x',y');

plot_points = 1000;
xx = linspace(-1,1,plot_points);
yy = runge(xx);

plot(xx,yy,'g');
hold on;
plot_cubic_spline(x,s0,s1,s2,s3);

您可以通过在 Matlab 提示符下运行以下命令来查看它的实际效果

 >> cubic_driver(5)
 >> clf 
 >> cubic_driver(10)
 >> clf
 >> cubic_driver(20)

当你有 20 个节点时，你的插值在视觉上与 Runge 函数没有区别。

对 Matlab 代码的一些评论：我不使用任何 for 或 while 循环。我能够矢量化所有操作。我用 快速形成稀疏的三对角矩阵spdiags。我使用反斜杠运算符解决它。我指望 Tim Davis 的 UMFPACK 来处理分解和正向和反向求解。

希望有帮助。该代码可作为 github 上的要点https://gist.github.com/1269709

于 2011-10-07T07:36:52.827 回答
3

样条函数中有一个错误，由 (n-2) 生成的 (n-2) 应该是对称的： 

lower = h(2:end);
main  = 2*(h(1:end-1) + h(2:end));
upper = h(1:end-1);

http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f3-3.pdf

于 2013-03-26T06:25:27.983 回答