1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
|
//Example 6.20
clc
clear
function [p] = cubicsplin(x,y)
// Fits point data to cubic spline fit
n = length(x);
a = y(1:n-1); // Spline Initials
M1 = zeros(3*(n-1));
M2 = zeros(3*(n-1),1);
// Point Substitutions
for i = 1:n-1
M1(i,i) = x(i+1) - x(i);
M1(i,i+n-1) = (x(i+1) - x(i))^2;
M1(i,i+2*(n-1)) = (x(i+1) - x(i))^3;
M2(i) = y(i+1) - y(i);
end
// Knot equations
for i = 1:n-2
// Derivative (S') continuity
M1(i+n-1,i) = 1;
M1(i+n-1,i+1) = -1;
M1(i+n-1,i+n-1) = 2*(x(i+1)-x(i));
M1(i+n-1,i+2*(n-1)) = 3*(x(i+1)-x(i))^2;
// S'' continuity
M1(i+2*n-3,i+n-1) = 2;
M1(i+2*n-3,i+n) = -2;
M1(i+2*n-3,i+2*(n-1)) = 6*(x(i+1)-x(i));
end
// Given BC
M1(3*n-4,n) = 1;
M1(3*n-3,2*n-2) = 1;
M1(3*n-3,3*n-3) = 3*(3-2);
var = M1\M2;
var = round(var);
b = var(1:n-1);
c = var(n:2*(n-1));
d = var(2*(n-1)+1:3*(n-1));
p = [d c b a(:)];
endfunction
x = 0:3;
y = [1 4 0 -2];
p = cubicsplin(x,y);
for i = 1:length(p(:,1))
disp(strcat(["S",string(i-1),"(x) ="]))
disp(poly(p(i,:),"X",["coeff"]))
end
|
