如何实现分段函数，然后在MATLAB中按一定间隔绘制它 [英] How to implement a piecewise function and then plot it on certain intervals in MATLAB

问题描述

我实际上正在尝试编写三次样条插值的代码.三次样条曲线可归结为一系列的n-1段，其中n是最初给出的原始坐标数，并且每个段均由某种三次函数表示.

I am actually attempting to write code for the cubic spline interpolation. Cubic spline boils down to a series of n-1 segments where n is the number of original coordinates given initially and the segments are each represented by some cubic function.

我已经弄清楚了如何获取存储在向量a,b,c,d中的每个段的所有系数和值，但是我不知道如何在不同的时间间隔上将函数绘制为分段函数.到目前为止，这是我的代码.最后一个for循环是我尝试绘制每个段的地方.

I have figured out how to get all the coefficients and values for each segment, stored in vectors a,b,c,d, but I don't know how to plot the function as a piecewise function on different intervals. Here is my code so far. The very last for loop is where I have attempted to plot each segment.

%initializations
x = [1 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6 7 8 9.2 10.5 11.3 11.6 12 12.6 13 13.3].';
y = [1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25].';

%n is the amount of coordinates
n = length(x);
%solving for a-d for all n-1 segments
a = zeros(n,1);
b = zeros(n,1);
d = zeros(n,1);

%%%%%%%%%%%%%% SOLVE FOR a's %%%%%%%%%%%%%
%Condition (b) in Definition 3.10 on pg 146
%Sj(xj) = f(xj) aka yj
for j = 1: n
    a(j) = y(j);
end

%initialize hj
h = zeros(n-1,1);

for j = 1: n-1
    h(j) = x(j+1) - x(j);
end

A = zeros(n,n);
bv = zeros(n,1); %bv = b vector

%initialize corners to 1
A(1,1) = 1;
A(n,n) = 1;

%set main diagonal
for k = 2: n-1
    A(k,k) = 2*(h(k-1) + h(k));
end

%set upper and then lower diagonals
for k = 2 : n-1
    A(k,k+1) = h(k); %h2, h3, h4...hn-1
    A(k,k-1) = h(k-1); %h1, h2, h3...hn
end

%fill up the b vector using equation in notes
%first and last spots are 0
for j = 2 : n-1
    bv(j) = 3*(((a(j+1)-a(j)) / h(j)) - ((a(j) - a(j-1)) / h(j-1)));
end

%augmented matrix
A = [A bv];



%%%%%%%%%%%% BEGIN GAUSSIAN ELIMINATION %%%%%%%%%%%%%%%
offset = 1;
%will only need n-1 iterations since "first" pivot row is unchanged
for k = 1: n-1
  %Searching from row p to row n for non-zero pivot
  for p = k : n
      if A(p,k) ~= 0;
          break;
      end
  end

  %row swapping using temp variable
  if p ~= k
      temp = A(p,:);
      A(p,:) = A(k,:);
      A(k,:) = temp;
  end


  %Eliminations to create Upper Triangular Form
  for j = k+1:n
     A(j,offset:n+1) = A(j,offset:n+1) - ((A(k, offset:n+1) * A(j,k)) / A(k,k));
  end
  offset = offset + 1;
end

c = zeros(n,1); %initializes vector of data of n rows, 1 column

%Backward Subsitution
%First, solve the nth equation
c(n) = A(n,n+1) / A(n,n);

%%%%%%%%%%%%%%%%% SOLVE FOR C's %%%%%%%%%%%%%%%%%%
%now solve the n-1 : 1 equations (the rest of them going backwards
for j = n-1:-1:1 %-1 means decrement
    c(j) = A(j,n+1);
    for k = j+1:n
        c(j) = c(j) - A(j,k)*c(k); 
    end
    c(j) = c(j)/A(j,j);
end

%%%%%%%%%%%%% SOLVE FOR B's and D's %%%%%%%%%%%%%%%%%%%%
for j = n-1 : -1 : 1
    b(j) = ((a(j+1)-a(j)) / h(j)) - (h(j)*(2*c(j) + c(j+1)) / 3);
    d(j) = (c(j+1) - c(j)) / 3*h(j);
end

%series of equation segments

for j = 1 : n-1
    f = @(x) a(j) + b(j)*(x-x(j)) + c(j)*(x-x(j))^2 + d(j)*(x-x(j))^3;
end
plot(x,y,'o');

让我们假设我已经为每个段正确计算了向量a,b,c,d.如何绘制每个立方段，以便它们全部显示在单个图上?

Let's assume that I have calculated vectors a,b,c,d correctly for each segment. How do I plot each cubic segment such that they all appear graphed on a single plot?

感谢您的帮助.

推荐答案

那很容易.通过为每个间隔之间的三次样条定义一个匿名函数，您已经完成了一半的工作.但是，您需要确保函数中的操作是 element-wise .您目前可以让它在标量上运行，或者假设我们正在使用矩阵运算.不要那样做使用.*代替*，并使用.^代替^.之所以需要这样做，是为了使样条曲线上的点生成变得容易得多，我的下一个点紧随其后.

That's pretty easy. You've already done half of the work by defining an anonymous function that is for the cubic spline in between each interval. However, you need to make sure that the operations in the function are element-wise. You currently have it operating on scalars, or assuming that we are using matrix operations. Don't do that. Use .* instead of * and .^ instead of ^. The reason why you need to do this is to make generating the points on the spline a lot easier, where my next point follows.

接下来要做的就是在相邻的x关键点定义的时间间隔内定义一堆x点，并将它们替换为函数，然后绘制结果..就像这样:

All you have to do next is define a bunch of x points within the interval defined by the neighbouring x key points and substitute them into your function, then plot the result.... so something like this:

figure;
hold on;
for j = 1 : n-1
    f = @(x) a(j) + b(j).*(x-x(j)) + c(j).*(x-x(j)).^2 + d(j)*(x-x(j)).^3; %// Change function to element-wise operations - be careful
    x0 = linspace(x(j), x(j+1)); %// Define set of points
    y0 = f(x0); %// Find output points
    plot(x0, y0, 'r'); %// Plot the line in between the key points
end
plot(x, y, 'bo');

我们生成一个新图形，然后使用hold on，这样当我们多次调用plot时，我们会将结果附加到同一图形上.接下来，对于每组三次样条系数，定义一个样条函数，然后使用x值. html"rel =" nofollow noreferrer> linspace ，它们位于当前x关键点与旁边的关键点之间.默认情况下，linspace在起点(即x(j))和终点(即x(j+1))之间生成100个点.您可以通过指定第三个参数来控制要生成的点数(例如linspace(x(j), x(j+1), 25);之类的以生成25个点).我们使用这些x值，并将它们代入样条方程式以获得我们的y值.然后，我们使用红线将此结果绘制在图形上.完成后，我们将关键点绘制为曲线上方的蓝色空心圆.

We spawn a new figure, then use hold on so that when we call plot multiple times, we append the results to the same figure. Next, for each set of cubic spline coefficients we have, define a spline function, then generate a bunch of x values with linspace that are between the current x key point and the one beside it. By default, linspace generates 100 points between a start point (i.e. x(j)) and end point (i.e. x(j+1)). You can control how many points you want to generate by specifying a third parameter (so something like linspace(x(j), x(j+1), 25); to generate 25 points). We use these x values and substitute them into our spline equation to get our y values. We then plot this result on the figure using a red line. Once we're done, we plot the key points as blue open circles on top of the curve.

作为奖励，我使用上述绘图机制运行了您的代码，这就是我得到的:

As a bonus, I ran your code with the above plotting mechanism, and this is what I get:

这篇关于如何实现分段函数，然后在MATLAB中按一定间隔绘制它的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！