Beetje tekenen, periospline functie is afgewerkt

Author: Xden21

Tante_Sidonia.gif

@@ -1,6 +1,6 @@
 % Excercise 2.2
 plotval = linspace(-pi, pi, 1000);
-x = linspace(-pi, pi, 50);
+x = linspace(-pi, pi, 20);
 t = linspace(-pi, pi, 200);
 
 func = @(x) sin(x) + sin(4.*x)./2;
@@ -9,11 +9,14 @@
 f = func(x(1:end-1));
 
 yEst = periospline(x,f,t);
-close;
+close all;
 plot(plotval, yplotval);
 hold on;
 plot(x(1:end-1), f, 'ro')
 plot(t, yEst, '-go');
+title('x=20')
+xlabel('x');
+ylabel('y')
 hold off;
 
 x1 = linspace(-pi,pi ,10);
@@ -22,23 +25,29 @@
 figure;
 plot(t, abs(func(t)-yEst1));
 title('x=10');
+xlabel('x');
+ylabel('error')
 
 x2 = linspace(-pi,pi ,50);
 f2 = func(x2(1:end-1));
 yEst2 = periospline(x2,f2,t);
 figure;
 plot(t, abs(func(t)-yEst2));
 title('x=50');
+xlabel('x');
+ylabel('error')
 
-errMax = zeros(1,50);
+errMax = zeros(1,40);
 index = 1;
-for i=10:10:500
+for i=5:5:200
     x = linspace(-pi,pi ,i);
     f = func(x(1:end-1));
     yEst = periospline(x,f,t);
     errMax(index) = max(abs(func(t)-yEst));
     index = index + 1;    
 end
 figure;
-plot(10:10:500, errMax);
-title('MaxErr');
+plot(5:5:200, errMax);
+title('MaxErr');
+xlabel('aantal of abscissen');
+ylabel('max fout');

excercise22.m

@@ -0,0 +1,9 @@
+% Excercise 2.3
+[x,y] = click(); % Gives us x_1 to x_n-1 and y_1 to y_n-1
+                  %x_n = x_1 and y_n = y_1, this happens implicitly in
+                  %periospline
+                  
+t = 1:1:(length(x)+1);
+evalX = linspace(1, length(x)+1, 1000);
+estimate = periospline(t,[x;y], evalX);
+plot(estimate(1,:), estimate(2,:));

excercise23.m

@@ -0,0 +1,8 @@
+% Exercise 2.4
+
+% Excercise 2.3
+open sidonia.mat;
+t = 1:1:(length(x)+1);
+evalX = linspace(1, length(x)+1, 1000);
+estimate = periospline(t,[x;y], evalX);
+plot(estimate(1,:), estimate(2,:));

excercise24.m

@@ -4,7 +4,7 @@
     n = length(x)-1;
     [d,~] = size(f);
     N = length(t);
-    y = zeros(d,N);
+    y = zeros(N, d);
     f = f';
 
     % Setup matrix & right hand side to determine s''(x_i)
@@ -20,12 +20,12 @@
             mainDiag(i-1) = 2.*((x(end) - x(end-1)) + (x(i) - x(i-1))); % 2*(x_n - x_n-1 + x_1 - x_0)
             downDiag(i-1) = x(i) - x(i-1); % x_1 - x_0
 
-            RHS(i-1,:) = ((f(i,:) - f(i-1,:))./(x(i) - x(i-1))) - ((f(i,:) - f(end,:))./(x(end)-x(end-1))); %f_n = f_0, f_n-1 = f(end)
+            RHS(i-1,:) = ((f(i,:) - f(i-1,:))./(x(i) - x(i-1))) - ((f(i-1,:) - f(end,:))./(x(end)-x(end-1))); %f_n = f_0, f_n-1 = f(end)
 
         elseif i == n+1
             % Other diags don't have values here
              mainDiag(i-1) = 2.*(x(i) - x(i-2)); % 2*(x_i-1 - x_i-2 + x_i - x_i-1)
-             RHS(i-1,:) = ((f(1,:) - f(i-1,:))./(x(i) - x(i-1))) - ((f(i-1,:) - f(i-2,:))./(x(i-1)-x(i-2))); %f_n = f_0, f_n-1 = f(end)
+             RHS(i-1,:) = ((f(1,:) - f(end,:))./(x(i) - x(i-1))) - ((f(end,:) - f(end-1,:))./(x(i-1)-x(i-2))); %f_n = f_0, f_n-1 = f(end)
 
         else
             upDiag(i-1) = x(i) - x(i-1); %x_i - x_i-1
@@ -42,11 +42,12 @@
     % Add corner coeff to matrix
     Coef(1,end) = x(end) - x(end-1);
     Coef(end,1) = x(end) - x(end-1);
+    RHS = 6.*RHS;
 
     % Solve all at once. This gives us the s''(x_i) for each function set
     % This is a matrix where each row gives us the s''(x_i) for a function set 
     % s''(x_0) to s''(x_n-1), s''(x_n) = s''(x_0)
-    SCoef = (Coef\RHS)
+    SCoef = (Coef\RHS);
 
     % For each p_i(x) determine c_1i & c_2i and put in array
     % We use a matrix where each colomn for one function set => n x d
@@ -60,7 +61,6 @@
         end
         C2(i-1, :) = f(i-1, :)./(x(i) - x(i-1)) - (x(i) - x(i-1))./6.0.*SCoef(i-1, :);
     end
-    
     % For eval determine in which interval the point lies, and use the
     % according p_i(x) (valid in interval [x_i-1, x_i])      
 
@@ -69,23 +69,20 @@
     for i = 1:N
        x_i = 2; % The top of the interval t(i) is in
        funcNum = 1; % The number of the spline function
-       if(t(i) > x(end)) || (t(i) < x(1))
-           % Figure out how to get it back in de base region (its periodic
-           % after all)
-           continue;
-       else
-           while(x_i < n+1) && (~(t(i) < x(x_i)) && (t(i)>x(x_i-1)))
-               x_i = x_i + 1;
-               funcNum = funcNum + 1;
-           end           
-       end
+       while(x_i < n+1) && (~(t(i) < x(x_i)) && (t(i)>x(x_i-1)))
+           x_i = x_i + 1;
+           funcNum = funcNum + 1;
+       end           
        if x_i == n+1
-           s_i = 2;
+           s_i = 1;
        else
            s_i = x_i;
        end
 
-       y(:, i) = ((t(i) - x(x_i-1)).^3)./(6.0.*(x(x_i) - x(x_i-1))).*SCoef(s_i) - ((t(i) - x(x_i)).^3)./(6.0*(x(x_i) - x(x_i -1))).*SCoef(x_i-1) + C1(funcNum, :)'.*(t(i) - x(x_i-1)) + C2(funcNum, :)'.*(x(x_i) - t(i));
+       y(i, :) = ((t(i) - x(x_i-1)).^3)./(6.0.*(x(x_i) - x(x_i-1))).*SCoef(s_i, :)...
+           - ((t(i) - x(x_i)).^3)./(6.0.*(x(x_i) - x(x_i -1))).*SCoef(x_i-1, :)...
+           + C1(funcNum, :).*(t(i) - x(x_i-1)) + C2(funcNum, :).*(x(x_i) - t(i));
     end    
+   y = y';
 end
 

periospline.m

@@ -0,0 +1,23 @@
+plotval = linspace(-pi, pi, 1000);
+plotval2 = linspace(pi, 3*pi, 1000);
+x = linspace(-pi, pi, 10);
+x2 = linspace(pi, 3*pi, 10);
+t = linspace(-pi, pi, 2000);
+t2 = linspace(pi, 3*pi, 2000);
+
+func = @(x) sin(x);
+
+yplotval = func(plotval);
+yplotval2 = func(plotval2);
+f = func(x(1:end-1));
+f2 =func(x2(1:end-1)); 
+yEst = periospline(x,f,t);
+yEst2 = periospline(x2,f2,t2);
+close;
+plot(plotval, yplotval, 'b');
+hold on;
+plot(plotval2, yplotval2, 'b');
+plot(x(1:end-1), f, 'ro')
+plot(t, yEst, 'go');
+plot(t2, yEst2,'go');
+hold off;

sidonia.mat

@@ -0,0 +1,52 @@
+classdef perioSplineTest < matlab.unittest.TestCase
+    %PERIOSPLINETEST Summary of this class goes here
+    %   Detailed explanation goes here
+     methods(TestMethodSetup)
+        function addPath(testCase)
+            addpath('../');
+        end
+    end
+    
+   
+    
+    methods(Test)
+        function testCoefMatrix(Test)
+            x = [1 3 7 15 31];
+            n = length(x)-1;
+            upDiag = zeros(n-1,1);
+            mainDiag = zeros(n,1);
+            downDiag = zeros(n-1,1);
+
+            for i =2:n+1
+                if i == 2
+                    % First one's special (used the last value)
+                    upDiag(i-1) = x(i) - x(i-1); % x_1 - x_0
+                    mainDiag(i-1) = 2.*((x(end) - x(end-1)) + (x(i) - x(i-1))); % 2*(x_n - x_n-1 + x_1 - x_0)
+                    downDiag(i-1) = x(i) - x(i-1); % x_1 - x_0
+
+                elseif i == n+1
+                    % Other diags don't have values here
+                     mainDiag(i-1) = 2.*(x(i) - x(i-2)); % 2*(x_i-1 - x_i-2 + x_i - x_i-1)
+                else
+                    upDiag(i-1) = x(i) - x(i-1); %x_i - x_i-1
+                    mainDiag(i-1) = 2.*(x(i) - x(i-2)); % 2*(x_i-1 - x_i-2 + x_i - x_i-1)
+                    downDiag(i-1) = x(i) - x(i-1); % x_i - x_i-1
+
+                end
+
+            end
+            Coef = full(gallery('tridiag', downDiag, mainDiag, upDiag));
+
+            % Add corner coeff to matrix
+            Coef(1,end) = x(end) - x(end-1);
+            Coef(end,1) = x(end) - x(end-1);
+            
+            Test.verifyEqual(Coef, [36 2  0 16; 2 12 4 0; 0 04 24 8; 16 0 8 48]); 
+        end
+        
+        function testRHSMatrix(Test)
+            
+        end
+    end
+end
+