function []=blatt9() n=50; f = @(U) inhalt(U); grad_f = @(U) gradinhalt(U); z0=zeros(n^2,1); for i=1:n for j=1:n z0((i-1)*n+j,1) = sin( (j-1)/(n-1)*(i-1)/(n-1)*pi ); end end B0 =10*eye(n^2); gamma = 0.01; eta=0.5; [z, fz] = Globalisiertes_BFGS( f, grad_f, z0, B0, gamma, eta); h=1/(n-1); x=0:h:1; y=0:h:1; [X,Y]=meshgrid(x,y); l=length(fz); Z=reshape(z(:,l),n,n); surf(X,Y,Z); end function [z, fz] = Globalisiertes_BFGS( f, grad_f, z0, B0, gamma, eta) i= 1; z(:,i) = z0; fz(i) = f(z(:,i)); grafz(:,i) = grad_f(z(:,i)); B = B0; while norm( grafz(:,i) ) > 10^(-6) d = - B* grafz(:,i); smin=1; while f( z(:,i) + smin*d ) - fz(i) > gamma*smin*(grafz(:,i))'*d smin=1/2*smin; end smax=2*smin; while f( z(:,i) + smax*d ) - fz(i) <= gamma*smax*(grafz(:,i))'*d smax=2*smax; end while (grad_f( z(:,i) + smin*d ))'*d < eta*( grafz(:,i))'*d s0 = (smin+smax)/2; if f( z(:,i) + s0*d ) - fz(i) <= gamma*s0*(grafz(:,i))'*d smin = s0; else smax = s0; end end s = smin; z(:,i+1) = z(:,i) + s*d; fz(i+1) = f(z(:,i+1)); grafz(:,i+1) = grad_f(z(:,i+1)); if norm( grafz(:,i+1)) > 10^(-6) d = z(:,i+1) - z(:,i); y = grafz(:,i+1) - grafz(:,i); B = B + ( (d - B*y)*d'+ d*(d -B*y)')/(d'*y) -( ( (d -B*y)'*y )/(d'*y)^2)*d*d'; else end i = i+1; end end function gra = gradinhalt(U) n=sqrt(length(U)); gra=zeros(n^2,1); for i=1:(n-2) for j=1:(n-2) gra(i*n + j+1, 1 ) = 1/(2*(n-1))*( ( U(i*n + j +1) - U((i-1)*n + j +1))/(sqrt(( U((i-1)*n+j+1)-U((i-1)*n + j))^2 + (U((i-1)*n+j+1) - U(i*n+j+1))^2+1/(n-1)^2)) + (U(i*n+j+1)-U(i*n+j))/sqrt( (U(i*n+j+1)-U(i*n+j))^2+(U(i*n+j)-U((i-1)*n+j))^2 + 1/(n-1)^2)+ (2*U(i*n+j+1) - U(i*n+j+2)-U((i-1)*n+j+1))/sqrt((U(i*n+j+2)-U(i*n+j+1))^2+(U(i*n+j+1)-U((i-1)*n+j+1))^2+1/(n-1)^2) + (2*U(i*n+j+1)-U(i*n+j)-U((i+1)*n+j+1))/sqrt((U(i*n+j+1)-U(i*n+j))^2+(U(i*n+j+1)-U((i+1)*n+j+1))^2+1/(n-1)^2)+ (U(i*n+j+1)-U(i*n+j+2))/sqrt((U(i*n+j+2)-U(i*n+j+1))^2+(U(i*n+j+2)-U((i+1)*n+j+2))^2+1/(n-1)^2) + (U(i*n+j+1)-U((i+1)*n+j+1))/sqrt((U((i+1)*n+j+2)-U((i+1)*n+j+1))^2+(U((i+1)*n+j+1)-U(i*n+j+1))^2+1/(n-1)^2) ); end end end function [A] = inhalt(U) n=sqrt(length(U)); for i=1:n for j=1:n x( (i-1)*n+j ) = (j-1)/(n-1); y( (i-1)*n+j ) = (i-1)/(n-1); z( (i-1)*n+j ) = U( (i-1)*n+j ); end end for i=1:(n-1) for j=1:(n-1) t1( 2*(n-1)*(i-1)+j ) = (i-1)*n+j ; t2( 2*(n-1)*(i-1)+j ) = (i-1)*n+j+1; t3( 2*(n-1)*(i-1)+j ) = i*n+j+1; t1( 2*(n-1)*(i-1)+(n-1)+j ) = (i-1)*n+j; t2( 2*(n-1)*(i-1)+(n-1)+j ) = i*n+j; t3( 2*(n-1)*(i-1)+(n-1)+j ) = i*n+j+1; end end A=0; for i=1:2*(n-1)^2 A = A + 1/2*norm( [ ( z(t2(i)) - z(t1(i)))*(y(t3(i))-y(t1(i))) - (y(t2(i)) - y(t1(i)))*(z(t3(i))-z(t1(i))); (x(t2(i))-x(t1(i)))*(z(t3(i))-z(t1(i))) - (z(t2(i))-z(t1(i)))*(x(t3(i))-x(t1(i)));(x(t2(i))-x(t1(i)))*(y(t3(i))-y(t1(i))) - (y(t2(i))-y(t1(i)))*(x(t3(i))-x(t1(i)))] ); end end