sympy.nsolve, mpmath.findroots how to pass args to the function?

[Scott]

I am currently using scipy.optimize.fsolve to solve a 15 and 420
equation nonlinear system but I would like to try methods other than
Powell's hybrid. The answer to the previous become a parameter for
the next step. The nonlinear system is polynomial and was created in
sympy.

sympy.nsolve, mpmath.findroots both provide a variety of solvers but
the interfaces seem to be less friendly regarding passing parameters
to the residual function and analytical jacobian.

Advice on how to get from :
>>> def f(x1, x2):
... return x1**2 + x2, 5*x1**2 - 3*x1 + 2*x2 - 3
...
>>> findroot(f, (0, 0))
matrix(
[['-0.618033988749895'],
['-0.381966011250105']])

to something like the following would be appreciated.
>>> def f(x1, x2,a,b):
... return a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 - 3
...
>>> findroot(f, parameters=(a,b),(0, 0))
matrix(
[['-0.618033988749895'],
['-0.381966011250105']])
V/R

Scott

[Ondrej Certik]

Either by modifying our solvers, I think we should provide this
"parameters" argument. If you could send us a patch fixing it, it'd be
awesome.

In the meantime, just use global variables, e.g. something like:

params = [ <whatever you need here> ]

def f(...):
global params
<use params>

global params
params = (a, b)
findroot(f, (0, 0))

Ondrej

[Vinzent Steinberg]

On Jan 31, 3:47 am, Scott <scotta_2...@yahoo.com> wrote:
> I am currently using scipy.optimize.fsolve to solve a 15 and 420
> equation nonlinear system but I would like to try methods other than
> Powell's hybrid.  The answer to the previous become a parameter for
> the next step.  The nonlinear system is polynomial and was created in
> sympy.

As far as I know scipy has more methods to solve nonlinear systems
than mpmath, which implements only Newton's method for
multidimensional functions. This is of course nice to get roots of any
precision.

> sympy.nsolve, mpmath.findroots both provide a variety of solvers but
> the interfaces seem to be less friendly regarding passing parameters
> to the residual function and analytical jacobian.

You can pass your own norm function if you want. Why would you want to
define your own residual? nsolve() needs indeed to be improved to
accept user specified jacobians or to calculate them numerically, at
the moment it's always calculated analytically which might be
inefficient for certain real-life problems.

>    Advice on how to get from :
>         >>> def f(x1, x2):
>         ...     return x1**2 + x2, 5*x1**2 - 3*x1 + 2*x2 - 3
>         ...
>         >>> findroot(f, (0, 0))
>         matrix(
>         [['-0.618033988749895'],
>          ['-0.381966011250105']])
>
> to something like the following would be appreciated.
>         >>> def f(x1, x2,a,b):
>         ...     return a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 - 3
>         ...
>         >>> findroot(f, parameters=(a,b),(0, 0))
>         matrix(
>         [['-0.618033988749895'],
>          ['-0.381966011250105']])

I think the cleanest way would be to define a new function, like

f2 = lambda x1, x2: f(x1, x2, a, b)

Not sure if I understood your problem correctly however.

Vinzent

[Scott]

Ondrej,

I will use your global parameter suggestion in the short term and take
a look at patching the code to allow an 'args' argument like in
scipy.optimize.fsolve .

Vinzent

By residual I meant the system of nonlinear equations I am trying to
find the roots of. Such as [a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 -
3]

>>> def f(x1, x2,a,b):
... return a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 - 3
...
>>> findroot(f, parameters=(a,b),(0, 0))
matrix(
[['-0.618033988749895'],
['-0.381966011250105']])

In scipy.optimize.fsolve I did not see any flags to use a non-default
system solver.
sympy and mpmath also have the benefit of be easily built with pure
python.

I am still having issues packaging my system in a format compatible
nsolve (no evalf) or findroot (cannot create mpf from array).

I am starting with :

import numpy
import sympy

sympy.var('x1,x2,a,b')

f1=list(a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 - 3)
x=list(x1,x2)
x0=numpy.array([.5,.5])
para=list(a,b)

a=3
b=2
x01=copy(x0)
nsolve(f1,x,x01)

How to best relate f1,x and para list is the question?

[Scott]

Making my function, F4, play interface with findroot required that I
be more verbose than desired.

(F4 runs under scipy fsolve)
def F4(dof2):
global dof1,m,G,a,dt
args=tuple(dof1)[:-cons]+tuple(dof2)+(m,G,a,dt)
return F3(*args)

dof1 is a list of symbols ,list([x12, z12, x22, z22, x32, z32, px12,
pz12, px22, pz22, px32, pz32, lam1, lam2, lam3])
F3 is the is the nonlinear system lambdified in terms of all
parameters and degrees of freedom.

If F4 is rewritten as F5 then sympy.mpmath.findroot(F5,list(dof1))
accepts the inputs and returns a proper solution.

def F5(x12, z12, x22, z22, x32, z32, px12, pz12, px22, pz22, px32,
pz32, lam1, lam2, lam3):
global dof1,m,G,a,dt
args=tuple(dof1)[:-cons]+(x12, z12, x22, z22, x32, z32, px12,
pz12, px22, pz22, px32, pz32, lam1, lam2, lam3)+(m,G,a,dt)
return F3(*args)

sympy.mpmath.findroot(F5,list(dof1))

Is there a clever way to get from F4 to F5 other than cutting and
pasting?

[smichr]

Along the lambda lines is this suggestion:

>>> from sympy import *
>>> def myfun(x, params):
... a, b = params
... return x**2 + Rational(a, b)
...
>>> params=[1,2]
>>> f = lambda x: myfun(x, p)
>>> f(1)
3/2
>>> params=[3,5]
>>> f(1)
8/5

On Feb 1, 8:30 am, Scott <scotta_2...@yahoo.com> wrote:
> Is there a clever way to get from F4 to F5 other than cutting and
> pasting?

Let python emit the function for you?

>>> def emit(globs, params):
... sg = str(globs)
... sp = str(params)
... return 'def F5%s:\n\tglobal %s\n\targs=tuple(dofl)[:-cons]+%s+%s\n
\treturn F3(*args)' % (sp, sg[1:-1], sg, sp)
...
>>> f5=emit(globs, params)
>>> exec f5 # you now have F5 in your namespace

I can say more if this looks like what you want. After the exec line
you have a function named F5 ready to use, i.e. you could pass F5 to a
solver:

>>> sympy.mpmath.findroot(F5,list(dof1))
etc...

BTW, I don't see where you've defined `cons`

[Scott]

smirch

Thanks for the tip, Your suggested strin manipulations along with the
exec function seem to work for a simple case.

Cheers

Scott

[Scott]

Below is my fast and dirty approach for making an mpmath friendly
nonlinear system from one configured for scipy.optimize.fsolve.

Thanks for the tips.

import sympy as sy
import numpy

sy.var('x1,x2,a,b')
global a,b
# assume f_list is a "long" list of "functions" generated by another
sympy script
# f_list may contain more than 100 terms
f_list=list([a*x1**2 + x2, b*5*x1**2 - 3*x1 + 2*x2 - 3])
x=list([x1,x2])#degrees of freedom or unknowns
parameters=list([a,b])#knowns
x0=numpy.array([1,1])#
args=x+parameters
b=1
a=2
# make the symbolic functions more 'numerical'
F1=sy.lambdify(args,f_list)
#reshuffle F1 as an fsolve freindly function with global parameters
def F2(x):
global a,b
parameters=list([a,b])
args=tuple(x)+tuple(parameters)
return F1(*args)

#F2(x)-> F3(x1,x2) for mpmath.findroot
f3='def F3(%s):\n\treturn F2(tuple(%s))'%(str(x)[1:-1],str(x))
exec f3

print(sy.mpmath.findroot(F3,list(x0)))