is it possible to just pass the approximate Jacobian to the derivative call? Since it's so complicated in my case to get the exact Jacobian I would like to do sth. like
du = [(F(u + h) - F(u - h))/(2*h)]
So you can't do this by passing something to derivative (which instead computes the jacobian for you given a residual).
[Buesing, Henrik] Yes, true. I mixed this up. I want to do solve(F == 0, u, J=approximate symbolic Jacobian)
You can instead of course build an approximate, symbolic, Jacobian "by hand" and provide that to the solver: you'll then end up doing quasi- newton, rather than newton.
[Buesing, Henrik] That's what I would like to do. In my multi-physics case (see attachement) the Jacobian would have a 2x2 block structure. I'm unsure how to provide this... I tried to write sth. down for the blocks (see test.py), but this is obviously wrong. Thank you for your help! Henrik