### PETSc and SNES in DOLFIN: Jacobian test fails with O(1) relative error

```
Testing hand-coded Jacobian, if the ratio is
O(1.e-8), the hand-coded Jacobian is probably correct.
Run with -snes_test_display to show difference
of hand-coded and finite difference Jacobian.
Norm of matrix ratio 1.00017, difference 31232.9 (user-defined state)
Norm of function ratio 1., difference 4371.8 (user-defined state)
Norm of matrix ratio 1.00017, difference 31232.9 (constant state -1.0)
Norm of function ratio 1.00015, difference 4411. (constant state -1.0)
Norm of matrix ratio 1.00017, difference 31232.9 (constant state 1.0)
Norm of function ratio 1.00015, difference 4411. (constant state 1.0)
```

After running, it exits with an error (normal behavior for test method). I am using a Jacobian computed automatically, so this result seems strange to me -- shouldn't the Jacobian computed with DOLFIN be exact, up to numerical precision? Is this the expected output from running the SNES test?

I actually thought that the whole point of this was linearization, and I would indeed expect the result to be strongly dependent on the initial guess/solution. Perhaps I'm misunderstanding something.

You only have to call `fenics.derivative`

once; but the nonlinear solver will always at some point have to evaluate the discrete Jacobian at a given solution point to form a linear system.

For my own use case, I am always solving a time-dependent initial boundary value problem (IBVP), so I simply initialize the solution either with the initial values of my IBVP or with the solution from the previous time step. That being said, I have observed at least once case where this somehow fails and leaving the default solution initialization (which should be zeros, at least when compiling with gcc) works miraculously.

```
from dolfin import *
import math
mesh = UnitSquareMesh(10,10)
# scalar case: output is easy to understand
S = FiniteElement("Lagrange",mesh.ufl_cell(),1)
# vector case: output uses an ugly index notation and requires a bit more
# patience to understand, but still eventually makes sense
#S = VectorElement("Lagrange",mesh.ufl_cell(),1)
V = FunctionSpace(mesh,S)
u = Function(V)
u.rename("u","u")
v = Function(V)
v.rename("N_A","N_A")
a = exp(sqrt(inner(u,v)))*dx
deriv_a = derivative(a,u)
u_trial = Function(V)
u_trial.rename("N_B","N_B")
action_deriv_a = action(deriv_a,u_trial)
print action_deriv_a
```

Could you perhaps try this same test with a more simple nonlinear PDE, and post the MWE (minimum working example)?

I'm interested in where this goes, since I've had some odd issues with

`fenics.derivative`

in the past, and never figured out what was wrong. (as of now,`fenics.derivative`

works fine for my problem)Also I'm guessing this isn't practical for your case, but: Can you derive the Jacobian (Gateaux derivative) by hand and test that? I did this a few times when I was having issues, and it's quite tedious.