2
$\begingroup$

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?

$\endgroup$
13
  • 1
    $\begingroup$ Have you checked Concus, Paul. "Numerical solution of the minimal surface equation." Mathematics of Computation (1967): 340-350.? JSTOR link $\endgroup$ – nicoguaro Aug 19 '15 at 20:46
  • $\begingroup$ What makes the problem more difficult in <put your language here>? $\endgroup$ – nicoguaro Aug 20 '15 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$ – SmallElephant Aug 20 '15 at 19:41
  • $\begingroup$ Hessian of gradient of what function? Have you solved any PDE using finite differences before? $\endgroup$ – nicoguaro Aug 20 '15 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$ – SmallElephant Aug 20 '15 at 20:17
1
$\begingroup$

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

% define square domain [-1,1] x [-1,1]
n = 51;
x=linspace(-1,1,n);
y=x;
[X,Y]=meshgrid(x,x);

% build finite differences operators
dx=x(2)-x(1);
e=ones(n,1);
d0x=ones(n,1);
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",...
    "JacobPattern",lapy+lapx+grady*gradx,"display","iter");

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

surf(x,x,reshape(f,n,n),"facecolor","interp")

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;
end

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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