function minimalsurface
M = 20;
h = 1.0/M;
[x,y] = meshgrid(0:h:1,0:h:1);
V = reshape(( (x).^2+(y).^2 )',(M+1)*(M+1),1);
[SOL,iter] = Newton(V,M,h);
figure(1);
surf(reshape(SOL,M+1,M+1)');
title(['Minimal Surface Example 1: Newton iteration = ' num2str(iter)]);
V = reshape(( (x-.5).^2+(y-.5).^2 )',(M+1)*(M+1),1);
[SOL,iter] = Newton(V,M,h);
figure(2);
surf(reshape(SOL,M+1,M+1)');
title(['Minimal Surface Example 2: Newton iteration = ' num2str(iter)]);
V = reshape(( cos(pi*(x-.5)) )',(M+1)*(M+1),1);
[SOL,iter] = Newton(V,M,h);
figure(3);
surf(reshape(SOL,M+1,M+1)');
title(['Minimal Surface Example 3: Newton iteration = ' num2str(iter)]);
V = reshape(( cos(pi*(x-.5))-cos(pi*(y-.5)) )',(M+1)*(M+1),1);
[SOL,iter] = Newton(V,M,h);
figure(4);
surf(reshape(SOL,M+1,M+1)');
title(['Minimal Surface Example 4: Newton iteration = ' num2str(iter)]);
function [Vnew, iter] = Newton(V,M,h)
for iter=1:1000
FV = F(V,M,h);
DirectionVector = ComputeDirectionVector(V,FV,M,h);
Vnew = V - DirectionVector;
marchingstepsize = 1.0;
while (norm(F(Vnew,M,h))>(1-(1e-4)*marchingstepsize)*norm(FV))
marchingstepsize = .5 * marchingstepsize;
Vnew = V - marchingstepsize * DirectionVector;
end
if (norm(Vnew-V)<1e-6)
break;
end
V = Vnew;
end
function [FV] = F(V,M,h)
v = reshape(V,M+1,M+1)';
ds = Ds(v,M,h);
fv = Dxf(Dxb(v,M,h)./ds,M,h) + Dyf(Dyb(v,M,h)./ds,M,h);
fv([1,end],:) = 0;
fv(:,[1,end]) = 0;
FV = reshape(fv',(M+1)*(M+1),1);
function [DirectionVector] = ComputeDirectionVector(V,FV,M,h)
DFV = -speye((M+1)*(M+1),(M+1)*(M+1));
vv = reshape(V,M+1,M+1)';
ds = 1./(Ds(vv,M,h)).^3/h^2;
dxdy = Dxb(vv,M,h).*Dyb(vv,M,h).*ds;
dxdx = ds.*(1+Dxb(vv,M,h).^2)-dxdy;
dydy = ds.*(1+Dyb(vv,M,h).^2)-dxdy;
for m=2:M
for l=2:M
row = (l-1)*(M+1) + m;
col = row + [-M-1, -M, -1, 1, M, M+1];
value = [dxdx(l,m), dxdy(l,m+1), dydy(l,m), ...
dydy(l,m+1), dxdy(l+1,m), dxdx(l+1,m) ];
DFV(row,col) = value;
DFV(row,row) = -sum(value);
end
end
DirectionVector = DFV \ FV;
function [dv] = Dxb(vv,M,h)
dv = [zeros(M+1,1), diff(vv,1,2)/h];
function [dv] = Dxf(vv,M,h)
dv = [diff(vv,1,2)/h, zeros(M+1,1)];
function [dv] = Dyb(vv,M,h)
dv = [zeros(1,M+1); diff(vv,1,1)/h];
function [dv] = Dyf(vv,M,h)
dv = [diff(vv,1,1)/h; zeros(1,M+1)];
function [dv] = Ds(vv,M,h)
dv = sqrt(1 + (Dxb(vv,M,h)).^2 + (Dyb(vv,M,h)).^2);