Matlab files from Lars

diff --git a/case1noupdate.m b/case1noupdate.m
new file mode 100644 (file)
index 0000000..5baa44c
--- /dev/null
+++ b/case1noupdate.m
@@ -0,0 +1,105 @@
+% Simple SR1 and Standard CG. 
+% Solution of min( x'Ax/2 - b'x )
+%-----------------------------------------------
+clf; 
+clear;
+close all;
+clc;
+%
+ndim = 10;   %input('Dimension of matrix:')  
+nsim = 15;   %input('Number of simulations:')
+npot = 0.5;  %input('Power of Rt*R:')             
+%
+R = randn(ndim);     
+A = (R'*R)^npot;      % A-matrix
+xsol = rand(ndim,1);  % x*
+b    = A*xsol;        % RHS
+kappa= max(eig(A))/min(eig(A)); % Condition number
+%
+% Initialization:
+x = 0 ; g = -b;
+B = eye(ndim); H = B;
+Norm2 = norm(xsol);
+NormA = sqrt(xsol'*A*xsol);
+
+%
+% Modified SR1 iterasjon
+for loop = 1:nsim
+       Bold = B;
+       Hold = H;
+       p = -H*g;
+       x = x + p;
+       y = -g;
+       g = A*x-b;
+       y = g + y;
+
+       v = y - B*p;
+       fprintf('It %2i: v''*p= %+ 2.10i ',loop,v'*p)
+       B = B + v*v'/(v'*p);
+       eigval = eig(B);
+       eigvalmin = min(eigval);
+       fprintf('eig= %+i',eigvalmin)
+       if eigvalmin < 0
+               fprintf(' No update\n')
+               B = Bold;
+       else
+               fprintf('\n')
+               w = p - H*y;
+               H = H + w*w'/(w'*y);
+       end
+
+       err2m(loop) = norm(x-xsol)/Norm2;
+       errAm(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
+end
+clf; 
+
+close all;
+
+%
+%
+% Initialization:
+x = 0 ; g = -b;
+B = eye(ndim); H = B;
+Norm2 = norm(xsol);
+NormA = sqrt(xsol'*A*xsol);
+
+%
+% Primitive SR1 iterasjon
+for loop = 1:nsim
+
+       p = -H*g;
+       x = x + p;
+       y = -g;
+       g = A*x-b;
+       y = g + y;
+       v = y - B*p;
+       B = B + v*v'/(v'*p);
+       w = p - H*y;
+       H = H + w*w'/(w'*y);
+
+       err2(loop) = norm(x-xsol)/Norm2;
+       errA(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
+end
+
+% Standard CG
+x = zeros(size(b)); g = -b ; p = -g; 
+% Iteration
+for loop = 1:nsim
+       Ap = A*p;        % Only one matrix-vector product!
+       alfa = -(p'*g)./(p'*Ap);
+       x = x + alfa*p;
+       g = g + alfa*Ap; % g = A*x-b;
+       beta = (g'*Ap)./(p'*Ap);
+       p = -g + beta*p;
+       err2CG(loop) = sqrt((x-xsol)'*(x-xsol))/Norm2;
+       errACG(loop) = sqrt((x-xsol)'*A*(x-xsol))/NormA;
+end; 
+semilogy(1:nsim, err2, 1:nsim, err2CG, ...
+1:nsim, errA, 1:nsim, errACG, 1:length(errAm), errAm, 'Linewidth', 2)        
+set(gca,'FontSize' ,14)    
+axis([0 nsim  1e-16 1])    
+legend( '2-norm SR1' , '2-norm CG','A-norm SR1 default' , 'A-norm CG','A-norm SR1 modified'  );
+xlabel('Iteration number'); ylabel('Normalized Error')
+Tittel = ['Size = ' num2str(ndim) ', Condition number = ',...
+num2str(kappa,'%8.2e\n')];
+title(Tittel);

diff --git a/constraints.m b/constraints.m
new file mode 100644 (file)
index 0000000..2ebe623
--- /dev/null
+++ b/constraints.m
@@ -0,0 +1,4 @@
+function [c,ceq] = constraints(x)
+c = [x(1)^2 + x(2)^2 - x(3)^2, -x(1)^2 - x(2)^2 - x(3)^2 + 4 , x(3) - 5, -x(1), -x(2), -x(3)];
+ceq = [];
+end

diff --git a/foppg1.m b/foppg1.m
new file mode 100644 (file)
index 0000000..d9c7ca5
--- /dev/null
+++ b/foppg1.m
@@ -0,0 +1,4 @@
+function [val] = foppg1(x)
+val = x(1)^3 - 6*x(1)^2 + 11*x(1) + x(3);
+end
+

diff --git a/oppg1.m b/oppg1.m
new file mode 100644 (file)
index 0000000..2b6db65
--- /dev/null
+++ b/oppg1.m
@@ -0,0 +1,47 @@
+clc; clear; close all;
+
+x0 = [2.1 0 2.1]';
+OPTIONS = optimset('Algorithm','active-set','TolX',1e-10,'TolCon',1e-10,'TolFun',1e-10);
+x = fmincon('foppg1',x0,[],[],[],[],[],[],'constraints',OPTIONS)
+
+% Plotting av figur (Klønete)
+x1 = 0:0.1:5;
+x2 = 0:0.1:5;
+[a,b,c] = sphere(50);
+a = 2*a; b =2*b; c = 2*c;
+a(a<0) = nan; b(b<0) = nan; c(c<0) = nan;
+for i=1:length(x1)
+    ok = 1;
+    for j=1:length(x2)
+        if sqrt(x1(i).^2 + x2(j).^2) <= 5
+            c1(i,j) = sqrt(x1(i).^2 + x2(j).^2);
+        else
+            c1(i,j) = 5;
+        end
+
+        if (4 - x1(i).^2 - x2(j).^2) >= 0
+            c2(i,j) = sqrt(4 - x1(i).^2 - x2(j).^2);
+        else
+            if ok == 1
+                c2(i,j) = 0;
+                ok = 0;
+            else
+                c2(i,j) = NaN;
+            end
+        end
+        c3(i,j) = 5;
+        f8(i,j) = 8 - x1(j).^3 + 6*x1(j).^2 - 11*x1(j);
+        f2(i,j) = 2 - x1(j).^3 + 6*x1(j).^2 - 11*x1(j);
+    end
+end
+hold on
+axis([0 5 0 5 0 5])
+surf(x2,x1,c1,'FaceColor',[1 0.5 0],'FaceAlpha',1);
+surf(a,b,c,'FaceColor','yellow')
+surf(x2,x1,c3,'FaceColor','red','FaceAlpha',1);
+mesh(x2(4:36),x1,f8(:,4:36),'EdgeColor','black','FaceColor','blue','FaceAlpha',0.5);
+mesh(x2(1:3),x1,f2(:,1:3),'EdgeColor','black','FaceColor','green','FaceAlpha',0.5);
+xlabel('x_1'); ylabel('x_2'); zlabel('x_3')
+scatter3([2,0,0,sqrt(2)],[0,sqrt(2),0,0],[2,sqrt(2),2,sqrt(2)],75,'white','filled','Marker','o','MarkerEdgeColor','black')
+legend('c_1 = 0','c_2 = 0','c_3 = 0','f = 8','f = 2','x_i')
+% Slutt plotting