I face to the following problem:

$$(1+u_x^2)u_{yy} - 2u_xu_yu_{xy} + (1+u_y^2)u_{xx}=0.$$ Problem needs to be discretized and assembled. Does anybody know how to proceed in Matlab?

  • 1
    $\begingroup$ Have you checked Concus, Paul. "Numerical solution of the minimal surface equation." Mathematics of Computation (1967): 340-350.? JSTOR link $\endgroup$
    – nicoguaro
    Commented Aug 19, 2015 at 20:46
  • $\begingroup$ What makes the problem more difficult in <put your language here>? $\endgroup$
    – nicoguaro
    Commented Aug 20, 2015 at 19:24
  • $\begingroup$ @nicoguaro in <Matlab> ... nothing more difficult - I guess it is same difficulty, I am just confused... after assemble, do I end up with one matrix and following system Ax = b? (I know - stupid question - I am really newbie in PDE and discretization in general)... also, afterwords - I need to obtain hessian and gradient - do u have any suggestions? $\endgroup$ Commented Aug 20, 2015 at 19:41
  • $\begingroup$ Hessian of gradient of what function? Have you solved any PDE using finite differences before? $\endgroup$
    – nicoguaro
    Commented Aug 20, 2015 at 19:42
  • 1
    $\begingroup$ seems like it is, you can write it in form as : $$A = \int \int (1 + u_x^2 + u_y^2)^{1/2} dx dy $$ and according to : en.wikipedia.org/wiki/Minimal_surface#cite_note-Meusnier1785-3 it should be $\endgroup$ Commented Aug 20, 2015 at 20:17

1 Answer 1


Here is a fast implementation using sparse matrices and sparse Jacobian estimation:

% define square domain [-1,1] x [-1,1]
n = 51;

% build finite differences operators
grad = spdiags([-e 0*e e],-1:1,n,n)/2/dx;
% use Kronecker product to build matrix of d/dx and d/dy
gradx = kron(grad,speye(n,n));
grady = kron(speye(n,n),grad);
lap = spdiags([e -2*e e],-1:1,n,n)/dx^2;
% use Kronecker product to build matrix of d/dx^2 and d/dy^2
lapx = kron(lap,speye(n,n));
lapy = kron(speye(n,n),lap);

% Dirichlet boundary condition
bdy = find(X(:)==x(1) | X(:)==x(end) | Y(:)==y(1) | Y(:)==y(end));
cnd = 0.5*sign(sin(4*(X(bdy)+Y(bdy))));

% initial value of iterate
f0 = zeros(n^2,1)

opt = optimoptions("fsolve","Algorithm","trust-region",...

f = fsolve(@(f) fun(f,gradx,grady,lapx,lapy,bdy,cnd),f0,opt);


function out=fun(f,gradx,grady,lapx,lapy,bdy,cnd)
    fx = gradx*f;
    fy = grady*f;
    % equation of mininal surface. Surface is the graph of z=f(x,y)
    % inside domain
    out = (1+fy.^2).*(lapx*f)+(1+fx.^2).*(lapy*f)-2*fx.*fy.*(gradx*(grady*f));
    % at the boundary
    out(bdy) = f(bdy)-cnd;

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.