This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Minimization

0 votes

Hello,

I would like to solve the following least square problem (part of Anderson acceleration fixed point iterations):
I have computed the residuals (Function(V)) R=(r1,r2, ..., rk) of a fixed point iteration
and I would like to find (a=a_1,a_2, ..., a_n) such that \sum a_i =1
and \|aR\|_(L^2) to be minimum.

I would be grateful for your help.

Best,
Fotini

asked Jan 13, 2016 by foteinikara FEniCS Novice (120 points)

1 Answer

0 votes

To find the minimum you differentiate the functional f with respect to dofs u
and solve df/du = 0.

answered Jan 15, 2016 by Kent-Andre Mardal FEniCS Expert (14,380 points)

Thank you for your answer, but I do not understand it. Can you please give me some more details because I am new in fenics and I do not know how to start?

I need to define a new function which is a linear combination of finite element functions,
that is I need to define
f(a_1,a_2, ..., a_m) = dot(a_1v_1+a_2v_2+...+a_m v_m, a_1v_1+a_2v_2+...+a_m v_m), and then to find the real values a_1,a_2, ..., a_m such that \sum a_i = 1 and f(a_1,a_2,...,a_m) is minimum.

First of all, I do not know how to define the function (the number m is given by the user and it is not fixed, so the length of (a_1,a_2, ..., a_m) is not fixed). Shall I use numpy arrays or something else? Secondly, I need to minimize with respect to real values (a_1,a_2, ..., a_m).

I hope that my problem is clearer now. Thank you again for your help!

Here is an example with minimization and constraint (but no PDE)

from dolfin import *
import math

mesh = UnitSquareMesh(1,1)

V = FunctionSpace(mesh, "DG", 0)
W = MixedFunctionSpace([V, V])

UU = Function(W)
UU.vector()[:] = 2.7

uu = TrialFunction(W)
k, l = TestFunctions(W)
U, V = split(UU)

f = (U2 + V2 - 4)kdx # minimization problem
g = (U-1)ldx # constraint 1, using Lagrange multiplier

F = f+g

J = derivative(F, UU, uu)

nl_problem = NonlinearVariationalProblem(F, UU, J=J)
solver = NonlinearVariationalSolver(nl_problem)

solver.parameters["newton_solver"]["relaxation_parameter"] = 1.0
solver.parameters["newton_solver"]["maximum_iterations"] = 40
solver.parameters["newton_solver"]["absolute_tolerance"] = 1.0e-6

solver.solve()

print "The true solution sqrt(2)=%e and the numerical solution is %e " % (math.sqrt(2), max(UU.vector().array()[:]))

...