Very basic question on Ch. 9


driebel
02-24-2009, 07:10 AM
I am attempting to find the roots to three equation in three unknowns, using the globally convergent Newton method. My problem is simply with the syntax. roots_multidim.h is looking for a &vecfunc object, and I have no idea how to create one. The only line from the text that seems relevant is "Its declaration as a function is
VecDoub vecfunc(VecDoub_I x); "

I have tried _many_ variations on the above line, and pre-main declaration, I won't reproduce them all here. Most produce ~50 lines of error messages and won't compile. My latest attempt, which does not work, is:

VecDoub vecfunc(VecDoub_I x) {
vecfunc[0] = x[0] * x[1] * x[2] - 6;
vecfunc[1] = x[0] * x[0] * x[1] + x[1] * x[1] * x[2]+x[2] * x[2] * x[0]-23;
vecfunc[2] = exp(x[0]+x[1]+x[2])-403;
}
using namespace std;
int main() {
const int N=3;
bool check;
int i;
VecDoub x(N);
VecDoub vecfunc(VecDoub x);
x[0] = 0;
x[1] = 5;
x[2] = 7;
newt(x,check,vecfunc);

As you can tell, I have minimal experience with C++. Any help on the basic syntax of how to create the necessary vector of functions would be much appreciated. Thanks,

Dave

davekw7x
02-24-2009, 10:21 AM
...My problem is simply with the syntax....
The function takes a vector input and returns a vector output.

#include "../code/nr3.h"
#include "../code/qrdcmp.h"
#include "../code/ludcmp.h"
#include "../code/roots_multidim.h"

VecDoub vecfunc(VecDoub_I x)
{
VecDoub y(3);
y[0] = x[0] * x[1] * x[2] - 6;

y[1] = x[0] * x[0] * x[1] +
x[1] * x[1] * x[2] +
x[2] * x[2] * x[0] -
23;

y[2] = exp(x[0] + x[1] + x[2]) - 403;

return y;
}

int main()
{
bool check;

VecDoub x(3);
x[0] = 0;
x[1] = 5;
x[2] = 7;

cout << "Initial guess: " << endl;
for (int i = 0; i < x.size(); i++) {
cout << " x[" << i << "] = " << x[i] << endl;
}
cout << endl;

cout << "Calling newt . . . " << endl << endl;
newt(x, check, vecfunc);

if (check) {
cout << "Check is true (convergence to a local minimum)." << endl;
cout << "Try another initial guess." << endl;
}
else {
cout << "Check is false (a \"normal\" return)." << endl;
}
cout << endl;

cout << scientific;
cout << "Solution from newt:" << endl;
for (int i = 0; i < x.size(); i++) {
cout << " x[" << i << "] = " << x[i] << endl;
}
cout << endl;

VecDoub y;
y = vecfunc(x);
cout << "Value of the solution vector:" << endl;
for (int i = 0; i < y.size(); i++) {
cout << " y[" << i << "] = " << y[i] << endl;
}
cout << endl;


return 0;
}


Output (GNU g++ version 4.1.2)

Initial guess:
x[0] = 0
x[1] = 5
x[2] = 7

Calling newt . . .

Check is false (a "normal" return).

Solution from newt:
x[0] = 1.305823e+00
x[1] = 1.391819e+00
x[2] = 3.301295e+00

Value of the solution vector:
y[0] = -1.472073e-13
y[1] = -5.798816e-13
y[2] = -2.330580e-12


Regards,

Dave

driebel
02-24-2009, 11:43 AM
I could have sworn I'd tried that, but I'd obviously left something out. That worked perfectly. Thank you so much!


Dave