6778
|
1 @c Copyright (C) 1996, 1997, 2007 John W. Eaton |
7018
|
2 @c |
|
3 @c This file is part of Octave. |
|
4 @c |
|
5 @c Octave is free software; you can redistribute it and/or modify it |
|
6 @c under the terms of the GNU General Public License as published by the |
|
7 @c Free Software Foundation; either version 3 of the License, or (at |
|
8 @c your option) any later version. |
|
9 @c |
|
10 @c Octave is distributed in the hope that it will be useful, but WITHOUT |
|
11 @c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
12 @c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
13 @c for more details. |
|
14 @c |
|
15 @c You should have received a copy of the GNU General Public License |
|
16 @c along with Octave; see the file COPYING. If not, see |
|
17 @c <http://www.gnu.org/licenses/>. |
3294
|
18 |
4167
|
19 @node Nonlinear Equations |
3294
|
20 @chapter Nonlinear Equations |
|
21 @cindex nonlinear equations |
|
22 @cindex equations, nonlinear |
|
23 |
|
24 Octave can solve sets of nonlinear equations of the form |
|
25 @iftex |
|
26 @tex |
|
27 $$ |
|
28 f (x) = 0 |
|
29 $$ |
|
30 @end tex |
|
31 @end iftex |
6850
|
32 @ifnottex |
3294
|
33 |
|
34 @example |
|
35 F (x) = 0 |
|
36 @end example |
6850
|
37 @end ifnottex |
3294
|
38 |
|
39 @noindent |
|
40 using the function @code{fsolve}, which is based on the @sc{Minpack} |
6850
|
41 subroutine @code{hybrd}. This is an iterative technique so a starting |
|
42 point will have to be provided. This also has the consequence that |
|
43 convergence is not guarantied even if a solution exists. |
3294
|
44 |
3368
|
45 @DOCSTRING(fsolve) |
3294
|
46 |
3368
|
47 @DOCSTRING(fsolve_options) |
3294
|
48 |
|
49 Here is a complete example. To solve the set of equations |
|
50 @iftex |
|
51 @tex |
|
52 $$ |
|
53 \eqalign{-2x^2 + 3xy + 4\sin(y) - 6 &= 0\cr |
|
54 3x^2 - 2xy^2 + 3\cos(x) + 4 &= 0} |
|
55 $$ |
|
56 @end tex |
|
57 @end iftex |
|
58 @ifinfo |
|
59 |
|
60 @example |
|
61 -2x^2 + 3xy + 4 sin(y) = 6 |
|
62 3x^2 - 2xy^2 + 3 cos(x) = -4 |
|
63 @end example |
|
64 @end ifinfo |
|
65 |
|
66 @noindent |
|
67 you first need to write a function to compute the value of the given |
|
68 function. For example: |
|
69 |
|
70 @example |
|
71 function y = f (x) |
|
72 y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6; |
|
73 y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4; |
|
74 endfunction |
|
75 @end example |
|
76 |
|
77 Then, call @code{fsolve} with a specified initial condition to find the |
|
78 roots of the system of equations. For example, given the function |
|
79 @code{f} defined above, |
|
80 |
|
81 @example |
6850
|
82 [x, info] = fsolve (@@f, [1; 2]) |
3294
|
83 @end example |
|
84 |
|
85 @noindent |
|
86 results in the solution |
|
87 |
|
88 @example |
|
89 x = |
|
90 |
|
91 0.57983 |
|
92 2.54621 |
|
93 |
|
94 info = 1 |
|
95 @end example |
|
96 |
6850
|
97 @noindent |
3294
|
98 A value of @code{info = 1} indicates that the solution has converged. |
|
99 |
|
100 The function @code{perror} may be used to print English messages |
|
101 corresponding to the numeric error codes. For example, |
|
102 |
|
103 @example |
|
104 @group |
|
105 perror ("fsolve", 1) |
|
106 @print{} solution converged to requested tolerance |
|
107 @end group |
|
108 @end example |
|
109 |
6850
|
110 When no Jacobian is supplied (as in the example above) it is approximated |
|
111 numerically. This requires more function evaluations, and hence is |
|
112 less efficient. In the example above we could compute the Jacobian |
|
113 analytically as |
3294
|
114 |
6850
|
115 @iftex |
|
116 @tex |
|
117 $$ |
|
118 \left[\matrix{ {\partial f_1 \over \partial x_1} & |
|
119 {\partial f_1 \over \partial x_2} \cr |
|
120 {\partial f_2 \over \partial x_1} & |
|
121 {\partial f_2 \over \partial x_2} \cr}\right] = |
|
122 \left[\matrix{ 3 x_2 - 4 x_1 & |
|
123 4 \cos(x_2) + 3 x_1 \cr |
|
124 -2 x_2^2 - 3 \sin(x_1) + 6 x_1 & |
|
125 -4 x_1 x_2 \cr }\right] |
|
126 $$ |
|
127 @end tex |
|
128 which is computed with the following Octave function |
|
129 @end iftex |
|
130 |
|
131 @example |
|
132 function J = jacobian(x) |
|
133 J(1,1) = 3*x(2) - 4*x(1); |
|
134 J(1,2) = 4*cos(x(2)) + 3*x(1); |
|
135 J(2,1) = -2*x(2)^2 - 3*sin(x(1)) + 6*x(1); |
|
136 J(2,2) = -4*x(1)*x(2); |
|
137 endfunction |
|
138 @end example |
|
139 |
|
140 @noindent |
|
141 Using this Jacobian is done with the following code |
|
142 |
|
143 @example |
|
144 [x, info] = fsolve (@{@@f, @@jacobian@}, [1; 2]); |
|
145 @end example |
|
146 |
|
147 @noindent |
|
148 which gives the same solution as before. |
|
149 |