在没有内置 MATLAB 函数的情况下构造具有自然边界条件的三次样条
Construction of a cubic spline with natural boundary conditions without built-in MATLAB functions
作为项目的一部分,我必须在不使用任何内置 MATLAB 函数(例如样条曲线或 csape)的情况下构建具有自然边界条件的三次样条曲线。
我尝试编写以下功能。
虽然我很确定它在计算系数 q 之前是正确的,但我无法弄清楚我最终将如何获得实际的三次多项式。当调用函数时,我现在得到的输出是 S 的 9 个不同值。
如有任何帮助或提示,我们将不胜感激。
function S=cubic_s(x,y)
n=length(x);
%construction of the tri-diagonal matrix
for j=1:n
V(j,1)=1;
V(j,2)=4;
V(j,3)=1;
end
%the first row should be (1,0,...,0) and the last (0,0,...,0,1)
V(1,2)=1; V(n,2)=1; V(2,3)=0; V(n-1,1)=0;
d=[-1 0 1];
A=spdiags(V,d,n,n);
%construction of the vector b
b=zeros(n,1);
%the first and last elements of b must equal 0
b(1)=0; b(n)=0;
%distance between two consecutive points
h=(x(n)-x(1))/(n-1);
for j=2:n-1
b(j,1)=(6/h^2)*(y(j+1)-2*y(j)+y(j-1));
end
%solving for the coefficients q
q=A\b;
%finding the polynomials with the formula for the cubic spline
for j=1:n-1
for z=x(j):0.01:x(j+1)
S(j)=(q(j)/(6*h))*(x(j+1)-z)^3+(q(j+1)/(6*h))*(z-x(j))^3+(z-x(j))* (y(j+1)/h-(q(j+1)*h)/6)+(x(j+1)-z)*(y(j)/h-(q(j)*h)/6);
end
end
你应该在每个 z-time 保存 S,见下面的图片和代码
function plot_spline
x = (0:10);
y = [1 4 3 7 1 5 2 1 6 2 3];
xx = x(1):0.01:x(2);
[XX,YY]=cubic_s(x,y);
plot(x,y,'*r', XX,YY,'-k')
function [XX,YY]=cubic_s(x,y)
n=length(x);
%construction of the tri-diagonal matrix
for j=1:n
V(j,1)=1;
V(j,2)=4;
V(j,3)=1;
end
%the first row should be (1,0,...,0) and the last (0,0,...,0,1)
V(1,2)=1; V(n,2)=1; V(2,3)=0; V(n-1,1)=0;
d=[-1 0 1];
A=spdiags(V,d,n,n);
%construction of the vector b
b=zeros(n,1);
%the first and last elements of b must equal 0
b(1)=0; b(n)=0;
%distance between two consecutive points
h=(x(n)-x(1))/(n-1);
for j=2:n-1
b(j,1)=(6/h^2)*(y(j+1)-2*y(j)+y(j-1));
end
%solving for the coefficients q
q=A\b;
%finding the polynomials with the formula for the cubic spline
enum = 1;
for j=1:n-1
for z=x(j):0.01:x(j+1)
YY(enum)=(q(j)/(6*h))*(x(j+1)-z)^3+(q(j+1)/(6*h))*(z-x(j))^3+(z-x(j))* (y(j+1)/h-(q(j+1)*h)/6)+(x(j+1)-z)*(y(j)/h-(q(j)*h)/6);
XX(enum)=z;
enum = enum+1;
end
end
作为项目的一部分,我必须在不使用任何内置 MATLAB 函数(例如样条曲线或 csape)的情况下构建具有自然边界条件的三次样条曲线。 我尝试编写以下功能。
虽然我很确定它在计算系数 q 之前是正确的,但我无法弄清楚我最终将如何获得实际的三次多项式。当调用函数时,我现在得到的输出是 S 的 9 个不同值。
如有任何帮助或提示,我们将不胜感激。
function S=cubic_s(x,y)
n=length(x);
%construction of the tri-diagonal matrix
for j=1:n
V(j,1)=1;
V(j,2)=4;
V(j,3)=1;
end
%the first row should be (1,0,...,0) and the last (0,0,...,0,1)
V(1,2)=1; V(n,2)=1; V(2,3)=0; V(n-1,1)=0;
d=[-1 0 1];
A=spdiags(V,d,n,n);
%construction of the vector b
b=zeros(n,1);
%the first and last elements of b must equal 0
b(1)=0; b(n)=0;
%distance between two consecutive points
h=(x(n)-x(1))/(n-1);
for j=2:n-1
b(j,1)=(6/h^2)*(y(j+1)-2*y(j)+y(j-1));
end
%solving for the coefficients q
q=A\b;
%finding the polynomials with the formula for the cubic spline
for j=1:n-1
for z=x(j):0.01:x(j+1)
S(j)=(q(j)/(6*h))*(x(j+1)-z)^3+(q(j+1)/(6*h))*(z-x(j))^3+(z-x(j))* (y(j+1)/h-(q(j+1)*h)/6)+(x(j+1)-z)*(y(j)/h-(q(j)*h)/6);
end
end
你应该在每个 z-time 保存 S,见下面的图片和代码
function plot_spline
x = (0:10);
y = [1 4 3 7 1 5 2 1 6 2 3];
xx = x(1):0.01:x(2);
[XX,YY]=cubic_s(x,y);
plot(x,y,'*r', XX,YY,'-k')
function [XX,YY]=cubic_s(x,y)
n=length(x);
%construction of the tri-diagonal matrix
for j=1:n
V(j,1)=1;
V(j,2)=4;
V(j,3)=1;
end
%the first row should be (1,0,...,0) and the last (0,0,...,0,1)
V(1,2)=1; V(n,2)=1; V(2,3)=0; V(n-1,1)=0;
d=[-1 0 1];
A=spdiags(V,d,n,n);
%construction of the vector b
b=zeros(n,1);
%the first and last elements of b must equal 0
b(1)=0; b(n)=0;
%distance between two consecutive points
h=(x(n)-x(1))/(n-1);
for j=2:n-1
b(j,1)=(6/h^2)*(y(j+1)-2*y(j)+y(j-1));
end
%solving for the coefficients q
q=A\b;
%finding the polynomials with the formula for the cubic spline
enum = 1;
for j=1:n-1
for z=x(j):0.01:x(j+1)
YY(enum)=(q(j)/(6*h))*(x(j+1)-z)^3+(q(j+1)/(6*h))*(z-x(j))^3+(z-x(j))* (y(j+1)/h-(q(j+1)*h)/6)+(x(j+1)-z)*(y(j)/h-(q(j)*h)/6);
XX(enum)=z;
enum = enum+1;
end
end