[SciPy-User] Request help on fsolve

Ralf Gommers ralf.gommers@gmail....
Mon Oct 1 12:42:12 CDT 2012


On Mon, Oct 1, 2012 at 12:41 AM, The Helmbolds <helmrp@yahoo.com> wrote:

> Please help me out here. I’m trying to rewrite the docstring for the
> `fsolve.py` routine
> located on my machine in: C:/users/owner/scipy/scipy/optimize/minpack.py
>
> The specific issue I’m having difficulty with is understanding the outputs
> described in fsolve’s docstring as:
>    'fjac': the orthogonal matrix, q, produced by the QR factorization of
> the final approximate Jacobian matrix, stored column wise
>    'r': upper triangular matrix produced by QR factorization of same matrix
>
> These are described in SciPy’s minpack/hybrd.f file as:
>    ‘fjac’ is an output n by n array which contains the orthogonal matrix q
> produced by the qr factorization of the final approximate jacobian.
>    ‘r’ is an output array of length lr which contains the upper triangular
> matrix produced by the qr factorization of the final approximate jacobian,
> stored rowwise.
>
> For ease in writing, in what follows let’s use the symbols ‘Jend’ for the
> final approximate Jacobian matrix, and use ‘Q’ and ‘R’ for its QR
> decomposition matrices. Now consider the problem of finding the solution to
> the following three nonlinear equations in three unknowns (u, v, w), which
> we will refer to as ‘E’:
>     2 * a * u + b * v + d - w * v = 0
>     b * u + 2 * c * v + e - w * u = 0
>     u * v - f = 0
> where (a, b, c, d, e, f ) = (2, 3, 7, 8, 9, 2). For inputs to fsolve, we
> identify (u, v, w) = (x[0], x[1], x[2]).
>
> Now fsolve gives the solution array:
>  [uend vend wend] = [  1.79838825   1.11210691  16.66195357].
> With these values, the above three equations E are satisfied to an
> accuracy of about 9 significant figures.
>
> The Jacobian matrix for the three LHS functions in E is:
>    J = np.matrix([[2*a, b-w, -v], [b-w, 2*c, -u], [v, u, 0.]])
> Note that it’s symmetric, and if we compute its value using the above
> fsolve’s ‘end’ solution values we get:
>  Jend = [[  4.          19.66195357   1.11210691],
>             [ 19.66195357  14.           1.79838825],
>             [  1.11210691   1.79838825   0.        ]]
> Using SciPy’s linalg package, this Jend has the QR decomposition:
>  Qend =  [[-0.28013447 -0.91516674 -0.28981807]
>             [ 0.95679602 -0.24168763 -0.16164302]
>             [ 0.07788487 -0.32257856  0.94333293]]
>  Rend =  [[-14.278857    17.08226116  -1.40915124]
>             [ -0.           9.69946027   1.45241144]
>             [ -0.           0.           0.61300558]]
> and Qend * Rend = Jend to within about 15 significant figures.
> However, fsolve gives the QR decomposition:
>   qretm =  [[-0.64093238  0.75748326  0.1241966 ]
>             [-0.62403598 -0.60841098  0.4903215 ]
>             [-0.44697291 -0.23675978 -0.8626471 ]]
>   rret =  [ -7.77806716  30.02199802  -0.819055   -10.74878184
> 2.00090268  1.02706198]
> and converting rret to a NumPy matrix gives:
>     rretm =  [[ -7.77806716  30.02199802  -0.819055  ]
>             [  0.         -10.74878184   2.00090268]
>             [  0.           0.           1.02706198]]
> Now qret and rretm bear no obvious relation to Qend and Rend. Although
> qretm is orthogonal to about 16 significant figures, we find the product:
>  qretm * rretm =  [[  4.98521509 -27.38409295   2.16816676]
>             [  4.85379376 -12.19513008  -0.2026608 ]
>             [  3.47658529 -10.87414051  -0.99362993]]
> which bears no obvious relationship to Jend.
>
> The hybrdj.f routine in minpack refers to a permutation matrix, p, such
> that we should have in our notation:
>  p*Jend = qretm*rretm,
> but fsolve apparently does not return the matrix p, and I don’t see any
> permutation of Jend that would equal qretm*rretm.
>
> If we reinterpret rret as meaning the matrix:
>  rretaltm =  [[ -7.77806716  30.02199802 -10.74878184]
>             [  0.          -0.819055     2.00090268]
>             [  0.           0.           1.02706198]]
> then we get the product:
>  qretm * rretaltm =  [[  4.98521509 -19.86249109   8.53245022]
>             [  4.85379376 -18.2364849    5.99384603]
>             [  3.47658529 -13.22510045   3.44468895]]
> which again bears no obvious relationship to Jend. Using the transpose of
> qretm in the above product is no help.
>
> So please help me out here. What are the fjac and r values that fsolve
> returns?
> How are they related to the above Qend, Rend, and Jend?
> How is the user supposed to use them?
>

I'm not sure. To play with your example it would be very helpful if you
could provide it as a Python script.

Ralf
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20121001/92c70e95/attachment.html 


More information about the SciPy-User mailing list