%=============================================================================
%(1) Vorbesetzungen
clear
Ordnung=4;
aa=0.05;b=0.4; %Intervall
tau=[aa,aa,0.06,0.06,0.075,0.075,0.09,0.09,0.13,0.13,0.21,0.21,b,b]; %14 Stück
taueinfach=[aa,0.06,0.075,0.09,0.13,0.21,b]; %Interpolationsknoten einfach
t=[aa aa tau b b]; %alle 18 Splineknoten  
%B_1 ist auf [t_1,t_5[, B_2 auf [t_2,t_6[, ..., B_{14} auf [t_{14},t_{18}] definiert.
%auszurechnen B_1(tau_1),   B_2(tau_1), ... ,   B_14(tau_1)
%             B_1'(tau_2),  B_2'(tau_2), ... ,  B_14'(tau_2)  %tau_2=tau_1,...
%             .
%             .
%             .
%             B_1(tau_13),  B_2(tau_13), ... ,  B_14(tau_13)
%             B_1'(tau_14), B_2'(tau_14), ... , B_14'(tau_14)  %tau_14=tau_13 
%             Im Programm heißt diese die Matrix Bneu
%Benötigt werden die folgenden Unterprogramme:

                 %B=bsplinesfct(k,Knoten,Knoteneinfach);
                 %alle B-Splines der Ordnung k zur Knotenfolge Knoten, 
                 %werden an den Stellen Knoteneinfach ausgerechnet.

                 %s=splinesfct(Knoten,a,k,x);
                 %Zu gegebenen Koeffizienten a, wird der Spline 
                 %s = sum a_jB_j(x) der Ordnung k an den Stellen x ausgerechnet. 
                 %Die Stellen x müssen in einem Intervall 
                 %[Knoten(j),Knoten(j+1)[  liegen, sonst Fehlermeldung.
%=============================================================================
%(2a) Berechnung der Matrix der Werte
B=bsplinesfct(Ordnung,t,taueinfach);      %B hat  die Größe 14 x 7.
%(2b) Berechnung der Matrix der Ableitungen
Bkm1=bsplinesfct(Ordnung-1,t,taueinfach); %B-Splines mit einer um 1 niedrigeren 
                                          %Ordnung. Gibt Splines Bkm1 mit Nummern 
                                          %von 1 bis 13. Das Progrann bsplinesfct 
                                          %reduziert die Anzahl der Endknoten 
                                          %immer auf Ordnung, in t sind dann 
                                          %also nur noch 16 Knoten, und es gibt
                                          %daher 16 - (Ordnung - 1) = 13 Zeilen.  
Bkm1=[zeros(1,length(taueinfach));Bkm1];  %vorne wird eine Zeile mit Nullen 
                                          %davorgehängt, da die vierfachen Knoten 
                                          %B_{1,3}=0 produzieren. 
                                          %Dieser Spline wird aber vom Programm 
                                          %bsplinesfct nicht ausgerechnet.
Bstrichh(1,:)=(Ordnung-1)*(-Bkm1(2,:)/(t(5)-t(2)));
for j=2:13
  Bstrichh(j,:)=(Ordnung-1)*(Bkm1(j,:)/(t(j+3)-t(j))-Bkm1(j+1,:)/(t(j+4)-t(j+1)));
end; %for j
Bstrichh(14,:)=(Ordnung-1)*(Bkm1(14,:)/(t(17)-t(14)));
Bneu=zeros(14,14); %Bneu wird aus B und Bstrichh zusammengesetzt
Bneu(1:2:13,:)=B';
Bneu(2:2:14,:)=Bstrichh';  %Die Matrix zerfällt in 4er Blöcke
%=============================================================================
%Alternativ zur Berechnung der Ableitung: %wird aber nicht gebraucht.
%I=eye(14);
%ast=[];
%for j=1:14
%  a=I(j,:);
%  as=astrich(t,a,4);
%  ast=[ast;as]; %Koeffizientenmatrix
%end; %for 
%=============================================================================
%(3) Jetzt Definition der rechten Seite
y=taueinfach.*sin(1./taueinfach);
ystrich=sin(1./taueinfach)-cos(1./taueinfach)./taueinfach;
yneu=zeros(1,14);
yneu(1:2:13)=y;
yneu(2:2:14)=ystrich;
yneu=yneu';
%=============================================================================
%(4) L\"osung des Gleichungssystems, L\"osung ist a
a=Bneu\yneu; 
%=============================================================================
%(5) Zeichnen des Splines (geht nur Intervallweise)
for j=1:6
  x=linspace(taueinfach(j),taueinfach(j+1));
  x(length(x))=[]; %halboffenes Intervall
  s=splinesfct(t,a,Ordnung,x); 
  plot(x,s,'b--'); if j==1, hold on, end; %if
end; %for
%=============================================================================
%(6) Zeichnen der interpolierten Funktion t*sin(1/t) in [0.05,0.4]
tt=linspace(aa,b,201);
yy=tt.*sin(1./tt);
plot(tt,yy,'r',[aa b],[0 0],'b',taueinfach,zeros(size(taueinfach)),'go')
hold off
%print -dpsc intspl.ps  %postscript-Version des (farbigen) Bildes herstellen
return
%=============================================================================
% (7) Ergebnisse:
% rechte Seite        Koeffizienten   lfd. Nr.
%=============================================
%    0.0456              0.0456          1
%   -7.2487              0.0215          2
%   -0.0491             -0.0783          3
%    8.7579             -0.0053          4
%    0.0520              0.0966          5
%
%   -8.9063              0.0075          6
%   -0.0894             -0.0780          7 
%   -2.2742             -0.1197          8
%    0.1283              0.1317          9
%   -0.2513              0.1216         10
%
%   -0.2097             -0.1768         11
%   -1.2345             -0.2879         12
%    0.2394              0.0746         13
%    2.6013              0.2394         14

Matrix 14 x 14
%==================================================================================================
%    1   0    0    0    0    0        0       0        0       0        0       0        0        0
% -300 300    0    0    0    0        0       0        0       0        0       0        0        0
%    0   0  0.6  0.4    0    0        0       0        0       0        0       0        0        0
%    0   0 -120  120    0    0        0       0        0       0        0       0        0        0
%    0   0    0    0  0.5  0.5        0       0        0       0        0       0        0        0
%    0   0    0    0 -100  100        0       0        0       0        0       0        0        0
%    0   0    0    0    0    0   0.7273  0.2727        0       0        0       0        0        0
%    0   0    0    0    0    0 -54.5455 54.5455        0       0        0       0        0        0
%    0   0    0    0    0    0        0       0   0.6667  0.3333        0       0        0        0
%    0   0    0    0    0    0        0       0 -25.0000 25.0000        0       0        0        0
%    0   0    0    0    0    0        0       0        0       0   0.7037  0.2963        0        0
%    0   0    0    0    0    0        0       0        0       0 -11.1111 11.1111        0        0
%    0   0    0    0    0    0        0       0        0       0        0       0        0   1.0000
%    0   0    0    0    0    0        0       0        0       0        0       0 -15.7895  15.7895
%==================================================================================================
