Minimal surface solution in Python

The question states that we need to solve a nonlinear partial differential equation. However Wikipedia states that 'They are difficult to study: there are almost no general techniques that work for all such equations, and usually each individual equation has to be studied as a separate problem.' However, you didn't give the equation! And does Matlab sometimes use genetic algorithms to arrive at its surfaces? That is, does it use a rule of thumb to make a best guess and then tries out small variations in the component squares until no smaller surface can be found. Implementing that kind of solution would be laborious but not conceptually difficult (assuming you like that sort of thing). Also remember that the calculus of continuous functions is just a special case of the calculus of all the linear approximations of functions (the increment is set to zero instead of some finite value). This was made clear to me by reading the books of J L Bell on smooth infinitesimal analysis - just use that algebra with finite increments and leave the resulting factors in the derivations instead of 'neglecting' them.

Obviously Matlab and SciPy understand TPS in different ways. The Matlab implementation looks correct.SciPy treats TPS the same way as others RBFs, so you could implement it correctly in Python yourself - it would be enough to form a matrix of the related linear equation system and solve it to receive TPS' coefficients.

You can use FEniCS:

from fenics import (    UnitSquareMesh,    FunctionSpace,    Expression,    interpolate,    assemble,    sqrt,    inner,    grad,    dx,    TrialFunction,    TestFunction,    Function,    solve,    DirichletBC,    DomainBoundary,    MPI,    XDMFFile,)# Create mesh and define function spacemesh = UnitSquareMesh(100, 100)V = FunctionSpace(mesh, "Lagrange", 2)# initial guess (its boundary values specify the Dirichlet boundary conditions)# (larger coefficient in front of the sin term makes the problem "more nonlinear")u0 = Expression("a*sin(2.5*pi*x[1])*x[0]", a=0.2, degree=5)u = interpolate(u0, V)print(    "initial surface area: {}".format(assemble(sqrt(1 + inner(grad(u), grad(u))) * dx)))# Define the linearized weak formulation for the Newton iterationdu = TrialFunction(V)v = TestFunction(V)q = (1 + inner(grad(u), grad(u))) ** (-0.5)a = (    q * inner(grad(du), grad(v)) * dx    - q ** 3 * inner(grad(u), grad(du)) * inner(grad(u), grad(v)) * dx)L = -q * inner(grad(u), grad(v)) * dxdu = Function(V)# Newton iterationtol = 1.0e-5maxiter = 30for iter in range(maxiter):    # compute the Newton increment by solving the linearized problem;    # note that the increment has *homogeneous* Dirichlet boundary conditions    solve(a == L, du, DirichletBC(V, 0.0, DomainBoundary()))    u.vector()[:] += du.vector()  # update the solution    eps = sqrt(        abs(assemble(inner(grad(du), grad(du)) * dx))    )  # check increment size as convergence test    area = assemble(sqrt(1 + inner(grad(u), grad(u))) * dx)    print(        f"iteration{iter + 1:3d}  H1 seminorm of delta: {eps:10.2e}  area: {area:13.5e}"    )    if eps < tol:        breakif eps > tol:    print("no convergence after {} Newton iterations".format(iter + 1))else:    print("convergence after {} Newton iterations".format(iter + 1))with XDMFFile(MPI.comm_world, "out.xdmf") as xdmf_file:    xdmf_file.write(u)

(Modified from http://www-users.math.umn.edu/~arnold/8445/programs/minimalsurf-newton.py.)