help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: fsolve for nonlinear system


From: Jaroslav Hajek
Subject: Re: fsolve for nonlinear system
Date: Mon, 7 Dec 2009 16:02:38 +0100

On Sun, Dec 6, 2009 at 1:50 PM, Bart Vandewoestyne
<address@hidden> wrote:
> Hello all,
>
> I'm trying to solve a nonlinear system using fsolve.  I'm using
> Octave 3.2.2 on Windows and I'm following the PDF docs for this
> version of Octave (the PDF that come installed together with
> Octave 3.2.2 on Windows).
>
> I define my function as:
>
>  function y = f_3d(x)
>
>  y(1) = 16*x(1)^4 + 16*x(2)^4 + x(3)^4 - 16;
>  y(2) =    x(1)^2 +    x(2)^2 + x(3)^2 - 3;
>  y(3) =    x(1)^3 -    x(2);
>
> then I also create my Jacobian as:
>
>  function J = jac_3d(x)
>
>  J = zeros(3,3);
>  J = [64*x(1)^3 64*x(2)^3 4*x(3)^3; ...
>        2*x(1)    2*x(2)   2*x(3);   ...
>        3*x(1)^2   -1      0];
>
> If i don't use the Jacobian, then all works well:
>
>  > x = [0.8; 0.6; 1.3];
>  octave-3.2.2.exe:9:C:\Octave\3.2.2_gcc-4.3.0\bin
>  > [x, fval, info] = fsolve(@f_3d, x0)
>  x =
>
>     0.87797
>     0.67676
>     1.33086
>
>  fval =
>
>    9.2213e-008  1.9237e-009  1.7637e-009
>
>  info =  1
>
> But if I try to use my own Jacobian i get errors:
>
>  > [x, fval, info] = fsolve(address@hidden, @jac_3d}, x0)
>  error: subscript indices must be either positive integers or logicals.
>  error: called from:
>  error:   
> C:\Octave\3.2.2_gcc-4.3.0\share\octave\3.2.2\m\optimization\fsolve.m at
>   line 177, column 6
>
>
> Can somebody tell me what I'm doing wrong here?  I can't see a difference
> between what I do and the manual...
>
> Kind regards,
> Bart
>

The style of supplying Jacobian is not supported in the 3.2 version of
fsolve. Instead, one should use the Matlab-compatible style:

 function [y,J] = f_3d(x)

 y(1) = 16*x(1)^4 + 16*x(2)^4 + x(3)^4 - 16;
 y(2) =    x(1)^2 +    x(2)^2 + x(3)^2 - 3;
 y(3) =    x(1)^3 -    x(2);

if (nargout > 1)
 J = zeros(3,3);
 J = [64*x(1)^3 64*x(2)^3 4*x(3)^3; ...
       2*x(1)    2*x(2)   2*x(3);   ...
       3*x(1)^2   -1      0];
endif


I think this style is preferable not only because of
Matlab-compatibility, but also because y and J are usually closely
related. It wouldn't be hard to support also the old version, but the
fact is that it is not there, because nobody asked for it.

Apparently the manual wasn't updated. But it has the docstring
included; so you should pay attention to that.

regards

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz



reply via email to

[Prev in Thread] Current Thread [Next in Thread]