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
|
function cubicspline(X,Y)
n = length(X);
//c(1) = 0;
//a = zeros(n-1,n-1);
for(i = 2:n-1)
c(i-1) = (6/(X(i+1)-X(i-1)))*((Y(i+1)-Y(i))/(X(i+1)-X(i)) - (Y(i)-Y(i-1))/ (X(i)-X(i-1)));
end
for(i = 2:n-1)
a(i-1,i-1) = (X(i)-X(i-1))/(X(i+1)-X(i-1));
a(i-1,i) = 2;
a(i-1,i+1) = (X(i+1)-X(i))/(X(i+1)-X(i-1));
end
[m,n] = size(a);
b = a(:,2:n-1); //For the case of natural splines, double derivative is zero at the first and last of data points. So removing the first and last columns since these are the only non-zero terms in the respective columns.
//[r2 c2] = size(b);
//disp(c)
//disp(b)
x = inv(b)*c;
//disp(x)
x1(1,1) = 0;
x1(n,1) = 0;
x1(2:n-1,1) = x(:,1);
printf('The values of second derivatives at the data points are :\n')
disp(x1)
x = poly(0,'x');
for(i = 2:n)
f(i-1) = [Y(i-1)*(X(i)-x)/(X(i)-X(i-1)) + Y(i)*(x-X(i-1))/(X(i)-X(i-1))] + x1(i-1)/6*[(X(i)-x)^3/(X(i)-X(i-1))-(X(i)-X(i-1))*(X(i)-x)] + x1(i)/6*[(x-X(i-1))^3/(X(i)-X(i-1))-(X(i)-X(i-1))*(x-X(i-1))];
end
printf('\nThe expressions for the cubic splines are \n')
disp(f)
for(i = 1:n-1)
p = X(i):0.01:X(i+1);
q = horner(f(i),p);
plot(p,q)
end
plot(X,Y,'ks')
xlabel('x')
ylabel('y')
endfunction
|
