fullmegalopp

%Παρακάτω ορίζεται συνάρτηση που δημιουργεί φυσική spline με δεδομένα
%τα σημεία (x,y) σε μορφή πινάκων-γραμμή. Τα σημεία χ πρέπει να ισαπέχουν.
%Η συνάρτηση επιστρέφει τις τιμές s της spline στα σημεία t (σε μορφή πινάκων)
function [t,s]=NaturalSpline(x,y)

n=size(x,2);      %όπου n ο αριθμός των σημείων Χ
h=abs(x(1)-x(2)); %υπολογίζεται το ισοδιάστημα h ως η διαφορά δύο διαδοχικών σημείων

a1=ones(n,1)-[zeros(n-2,1);1;1]; %κατασκευάζονται οι 3 στήλες που θα αποτελέσουν
a2=4*ones(n,1)-3*[1;zeros(n-2,1);1]; %τις διαγώνιες του Α. Ο Α για κάθε n δίνεται
a3=ones(n,1)-[1;1;zeros(n-2,1)]; %από τη θεωρία των φυσικών spline.

A=spdiags([a1 a2 a3],[-1,0,1],n,n); %κατασκευάζεται ο τριδιαγώνιος πίνακας του συστήματος
%για τους συντελεστές d(i) της spline χρησιμοποιώντας την εντολή spdiags για να
%αξιοποιήσουμε την τριδιαγώνια αραιή μορφή του Α.

for ii=1:n-2

  b1(ii,1)=(6/(h^2))*(y(ii)-2*y(ii+1)+y(ii+2));
endfor
b=[0; b1; 0]; %ορίζεται ο κάθετος πίνακας b του δεξιού μέλους του συστήματος, βάσει θεωρίας

d=A\b; %επιλύεται το σύστημα

%Τέλος, χρησιμοποιούνται οι συντελεστές d(i) που υπολογίστηκαν (πίνακας d) για τον
%υπολογισμό της spline στα σημεία t (20 σημεία ανά υποδιάστημα)
t=linspace(x(1),x(n),20*(n-1));
for ii=1:n-1
  for jj=1:20
    s(20*(ii-1)+jj)=(d(ii)/(6*h))*(x(ii+1)-t(20*(ii-1)+jj))^3+(d(ii+1)/(6*h))*(t(20*(ii-1)+jj)-x(ii))^3+((y(ii+1)/h)-((d(ii+1)*h)/6))*(t(20*(ii-1)+jj)-x(ii))+((y(ii)/h)-((d(ii)*h)/6))*(x(ii+1)-t(20*(ii-1)+jj));
  endfor

endfor


endfunction