如何实现分段函数,然后在 MATLAB 中按特定间隔绘制它
How to implement a piecewise function and then plot it on certain intervals in MATLAB
我实际上是在尝试为三次样条插值编写代码。三次样条归结为一系列 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
。如何绘制每个立方体线段,使它们都显示在一个图上?
感谢您的帮助。
这很简单。通过为每个区间之间的三次样条定义一个匿名函数,您已经完成了一半的工作。但是,您需要确保函数中的操作是 element-wise。您目前让它在标量上运行,或者假设我们正在使用矩阵运算。不要那样做。使用 .*
代替 *
,使用 .^
代替 ^
。您需要这样做的原因是为了更容易生成样条曲线上的点,我的下一点如下。
接下来您要做的就是在相邻 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
在起点(即 x(j)
)和终点(即 x(j+1)
)之间生成 100 个点。您可以通过指定第三个参数来控制要生成的点数(例如 linspace(x(j), x(j+1), 25);
生成 25 个点)。我们使用这些 x
值并将它们代入我们的样条方程以获得我们的 y
值。然后,我们使用红线将此结果绘制在图中。完成后,我们将关键点绘制为曲线顶部的蓝色空心圆圈。
作为奖励,我 运行 您的代码具有上述绘图机制,这就是我得到的:
我实际上是在尝试为三次样条插值编写代码。三次样条归结为一系列 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
。如何绘制每个立方体线段,使它们都显示在一个图上?
感谢您的帮助。
这很简单。通过为每个区间之间的三次样条定义一个匿名函数,您已经完成了一半的工作。但是,您需要确保函数中的操作是 element-wise。您目前让它在标量上运行,或者假设我们正在使用矩阵运算。不要那样做。使用 .*
代替 *
,使用 .^
代替 ^
。您需要这样做的原因是为了更容易生成样条曲线上的点,我的下一点如下。
接下来您要做的就是在相邻 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
在起点(即 x(j)
)和终点(即 x(j+1)
)之间生成 100 个点。您可以通过指定第三个参数来控制要生成的点数(例如 linspace(x(j), x(j+1), 25);
生成 25 个点)。我们使用这些 x
值并将它们代入我们的样条方程以获得我们的 y
值。然后,我们使用红线将此结果绘制在图中。完成后,我们将关键点绘制为曲线顶部的蓝色空心圆圈。
作为奖励,我 运行 您的代码具有上述绘图机制,这就是我得到的: