Как реализовать кусочную функцию, а затем построить ее на определенных интервалах в MATLAB

Question

Я на самом деле пытаюсь написать код для интерполяции кубического сплайна. Кубический сплайн сводится к серии n-1 сегменты где n это число исходных координат, заданных изначально, и каждый сегмент представлен некоторой кубической функцией.

Я выяснил, как получить все коэффициенты и значения для каждого сегмента, сохраненные в векторах a,b,c,d, но я не знаю, как построить функцию как кусочную функцию на разных интервалах. Вот мой код до сих пор. Самый последний цикл for - это место, где я пытался построить каждый сегмент.

%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 правильно для каждого сегмента. Как мне построить каждый кубический сегмент так, чтобы все они отображались на одном графике?

Я ценю помощь.

Answer 1

Это довольно легко. Вы уже сделали половину работы, определив анонимную функцию для кубического сплайна между каждым интервалом. Однако вам нужно убедиться, что операции в функции поэлементны. В настоящее время он работает на скалярах или предполагает, что мы используем матричные операции. Не делай этого. использование .* вместо * а также .^ вместо ^, Причина, по которой вам нужно это сделать, состоит в том, чтобы значительно упростить создание точек на сплайне, откуда следует мой следующий пункт.

Все, что вам нужно сделать, это определить группу x точки в интервале, определенном соседним x ключевые точки и подставьте их в свою функцию, затем нарисуйте результат.... так что-то вроде этого:

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 значения с linspace которые находятся между текущими x ключевой момент и тот, что рядом с ним. По умолчанию, linspace генерирует 100 точек между начальной точкой (т.е. x(j)) и конечная точка (т.е. x(j+1)). Вы можете контролировать, сколько точек вы хотите сгенерировать, указав третий параметр (так что-то вроде linspace(x(j), x(j+1), 25); сгенерировать 25 баллов). Мы используем эти x значения и подставить их в наше уравнение сплайна, чтобы получить наш y ценности. Затем мы наносим этот результат на фигуру, используя красную линию. Как только мы закончим, мы нанесем ключевые точки в виде синих открытых кружков на вершине кривой.

В качестве бонуса я запустил ваш код с помощью вышеуказанного механизма построения графиков, и вот что я получаю: