[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