[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