三次样条程序 [英] Cubic Spline Program
问题描述
我正在尝试编写三次样条插值程序.我已经编写了程序，但是图形无法正确显示.样条线使用自然边界条件(开始/结束节点处的第二个导数为0).该代码在Matlab中，如下所示，
I'm trying to write a cubic spline interpolation program. I have written the program but, the graph is not coming out correctly. The spline uses natural boundary conditions(second dervative at start/end node are 0). The code is in Matlab and is shown below,
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)
运行代码时，插值函数不对称，并且无法正确收敛.谁能在我的代码中提供有关问题的任何建议?谢谢.
When I run the code, the interpolation function is not symmetric and, it doesn't converge correctly. Can anyone offer any suggestions about problems in my code? Thanks.
推荐答案
很久以前，我在Mathematica中编写了三次样条程序包.这是我将该软件包翻译成Matlab的内容.请注意，大约7年以来我一直没有看过三次样条曲线，因此我将其作为我自己的文档的基础.你应该检查我说的一切.
I wrote a cubic spline package in Mathematica a long time ago. Here is my translation of that package into Matlab. Note I haven't looked at cubic splines in about 7 years, so I'm basing this off my own documentation. You should check everything I say.
基本问题是我们得到了n个数据点(x(1), y(1)) , ... , (x(n), y(n))，我们希望计算分段三次插值.插值定义为
The basic problem is we are given n data points (x(1), y(1)) , ... , (x(n), y(n)) and we wish to calculate a piecewise cubic interpolant. The interpolant is defined as
S(x) = { Sk(x) when x(k) <= x <= x(k+1)
{ 0 otherwise
Sk(x)是形式为三次多项式
Here Sk(x) is a cubic polynomial of the form
Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3
样条曲线的属性是:
- 样条线穿过数据点
Sk(x(k)) = y(k) - 样条曲线在端点处是连续的，因此在插值区间
Sk(x(k+1)) = Sk+1(x(k+1))中到处都是连续的. - 样条具有连续的一阶导数
Sk'(x(k+1)) = Sk+1'(x(k+1)) - 样条具有连续的二阶导数
Sk''(x(k+1)) = Sk+1''(x(k+1))
- The spline pass through the data point
Sk(x(k)) = y(k) - The spline is continuous at the end-points and thus continuous everywhere in the interpolation interval
Sk(x(k+1)) = Sk+1(x(k+1)) - The spline has continuous first derivative
Sk'(x(k+1)) = Sk+1'(x(k+1)) - The spline has continuous second derivative
Sk''(x(k+1)) = Sk+1''(x(k+1))
要从一组数据点构建三次样条，我们需要求解系数
每个n-1三次多项式的sk0，sk1，sk2和sk3.这总共是4*(n-1) = 4*n - 4个未知数.属性1提供n约束，属性2、3、4每个提供一个附加的n-2约束.因此，我们有n + 3*(n-2) = 4*n - 6约束和4*n - 4未知数.这留下了两个自由度.我们通过在起点和终点将二阶导数设置为零来修正这些自由度.
To construct a cubic spline from a set of data point we need to solve for the coefficients
sk0, sk1, sk2 and sk3 for each of the n-1 cubic polynomials. That is a total of 4*(n-1) = 4*n - 4 unknowns. Property 1 supplies n constraints, and properties 2,3,4 each supply an additional n-2 constraints. Thus we have n + 3*(n-2) = 4*n - 6 constraints and 4*n - 4 unknowns. This leaves two degrees of freedom. We fix these degrees of freedom by setting the second derivative equal to zero at the start and end nodes.
让m(k) = Sk''(x(k))，h(k) = x(k+1) - x(k)和d(k) = (y(k+1) - y(k))/h(k).以下
三项递归关系成立
Let m(k) = Sk''(x(k)) , h(k) = x(k+1) - x(k) and d(k) = (y(k+1) - y(k))/h(k). The following
three-term recurrence relation holds
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的系数由
The m(k) are unknowns we wish to solve for. The h(k) and d(k) are defined by the input data.
This three-term recurrence relation defines a tridiagonal linear system. Once the m(k) are determined the coefficients for Sk are given by
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中:
Okay that is all the math you need to know to completely define the algorithm to compute a cubic spline. Here it is in 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);
下面是一些绘制三次样条的代码:
Here is some code to plot a cubic spline:
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
这是一个构造三次样条并在著名的Runge函数上绘图的函数:
Here is a function that constructs a cubic spline and plots in on the famous Runge function:
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提示符下运行以下命令来查看它的运行情况
You can see it in action by running the following at the Matlab prompt
>> cubic_driver(5)
>> clf
>> cubic_driver(10)
>> clf
>> cubic_driver(20)
当您有二十个节点时，插值与Runge函数在视觉上是无法区分的.
By the time you have twenty nodes your interpolant is visually indistinguishable from the Runge function.
关于Matlab代码的一些注释:我不使用任何for或while循环.我能够向量化所有操作.我很快用spdiags形成了稀疏的三对角矩阵.我使用反斜杠运算符解决它.我指望蒂姆·戴维斯(Tim Davis)的UMFPACK来处理分解以及向前和向后求解.
Some comments on the Matlab code: I don't use any for or while loops. I am able to vectorize all operations. I quickly form the sparse tridiagonal matrix with spdiags. I solve it using the backslash operator. I counting on Tim Davis's UMFPACK to handle the decomposition and forward and backward solves.
希望有帮助.该代码可在github https://gist.github.com/1269709
Hope that helps. The code is available as a gist on github https://gist.github.com/1269709
这篇关于三次样条程序的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
