Using Newton Newton-Raphson iteration to solve a system on non-linnear equations

Hello,

I'm a student in a matlab class and I'm having trouble with a problem. I am not looking for someone to give me the solution, but rather to help me troubleshoot why my solution isn't working.

Here is the problem:

Consider the following system of equations:

Axy-(x^2 + y^2)^2 =0 y-Bx^3 = 0

where the constants A and B are parameters. Write a function that uses Newton-Raphson iteration to solve for x and y given values for the two parameters in the system. Your function should accept the following inputs (in order):

A scalar value for A. A scalar value for B. A column vector of initial guesses for x and y . A stopping criterion for the iteration. Your function should return the following outputs (in order):

A two-element column vector of the solutions for x and y. The Jacobian matrix evaluated at the initial guesses. (A 2x2 matrix). A two-element column vector of residuals associated with the solution. The number of iterations required for convergence (scalar). Set maximum to 50. Notes: Do not rearrange the equations before formulating the Jacobian if you want your code to pass the test suites. Also, the example NRsys.m will load with the test suite so you can use a function call to NRsys in your own solution if you like.

-------------------------------------------------------------------------------------

Here is my purposed function:

function [x_sol, J_initial, residuals, iter_count] = nr_nonlinear(A, B, x_init, es)

%Establish Max iteration count
max_it=50;

%Establish initial function
f = @(x) [A*x(1)*x(2)-(x(1)^2 +x(2)^2)^2;x(2)-B*x(1)^3];

f1 = @(x) A*x(1)*x(2)-(x(1)^2 +x(2)^2)^2;
f2 = @(x) x(2)-B*x(1)^3;

del = 1e-6;
df1x = @(x)( (f1(x(1)+del*x(1))) - (f1(x(1))) ) / (del*x(1));
df1y = @(x)( (f1(x(2)+del*x(2))) - (f1(x(2))) ) / (del*x(2));
df2x = @(x)( (f2(x(1)+del*x(1))) - (f2(x(1))) ) / (del*x(1));
df2y = @(x)( (f2(x(2)+del*x(2))) - (f2(x(2))) ) / (del*x(2));

%Build Jacobi Matrix
J = @(x) [df1x(x) , df1y(x); df2x(x) , df2y(x)];
J_initial = J(x_init);

%Apply NRsys
[x_sol,residuals, ea ,iter_count] = NRsys(f, J,x_init,es,max_it); 

And here is the list of errors it produces:

EDU>> Test
Attempted to access x(2); index out of bounds because numel(x)=1.

Error in nr_nonlinear>@(x)A*x(1)*x(2)-(x(1)^2+x(2)^2)^2 (line 22)
f1 = @(x) A*x(1)*x(2)-(x(1)^2 +x(2)^2)^2;

Error in nr_nonlinear>@(x)((f1(x(1)+del*x(1)))-(f1(x(1))))/(del*x(1)) (line 26)
df1x = @(x)( (f1(x(1)+del*x(1))) - (f1(x(1))) ) / (del*x(1));

Error in nr_nonlinear>@(x)[df1x(x),df1y(x);df2x(x),df2y(x)] (line 33)
J = @(x) [df1x(x) , df1y(x); df2x(x) , df2y(x)];

Error in nr_nonlinear (line 34)
J_initial = J(x_init);

Error in Test (line 1)
[x_sol, J_initial, residuals, iter_count] = nr_nonlinear(5, 1.5, [1;2], 1e-8)

Can anyone tell me what I'm doing wrong?

Answer 1

[Geoff Hayes]

Dylan - the error

 Attempted to access x(2); index out of bounds because numel(x)=1.

is telling you that your code is treating x as if it is an array with at least two elements when in fact it is a single element (scalar). Look at how you have defined your two functions

 f1 = @(x) A*x(1)*x(2)-(x(1)^2 +x(2)^2)^2;
 f2 = @(x) x(2)-B*x(1)^3;

Each of the above assumes that the input x is a 2 element array. Now look at how the code calls these functions. Given your input x_init which is (presumably) a two element vector, the Jacobian is built as

 J_initial = J(x_init);

where

 J = @(x) [df1x(x) , df1y(x); df2x(x) , df2y(x)];

So the two element array is passed to each of the partial derivatives for your two functions as

 df1x = @(x)( (f1(x(1)+del*x(1))) - (f1(x(1))) ) / (del*x(1));
 df1y = @(x)( (f1(x(2)+del*x(2))) - (f1(x(2))) ) / (del*x(2));
 df2x = @(x)( (f2(x(1)+del*x(1))) - (f2(x(1))) ) / (del*x(1));
 df2y = @(x)( (f2(x(2)+del*x(2))) - (f2(x(2))) ) / (del*x(2));

Each receives the two element array BUT only passes a single element into the f1 and f2 functions when these two functions are expecting a two element array. Hence the error.

It's been a while since I've looked at a problem like this, but I think that you want to modify each of your partial derivative functions so that only the variable being differentiated is affected by the small value del. For example, the first equation would become

 df1x = @(x)( (f1([x(1)+del x(2)])) - (f1(x)) ) / (del);

Note how the input the left f1 is an array of two elements where we add the small value (usually denoted by h) to just the first element x(1). We use the square brackets around

 [x(1)+del x(2)]

to indicate that we are concatenating the modified x(1) with x(2) so that f1 receives the two element array. For the right most call to f1, we just pass in x (the untouched two element array).

I've removed the multiplication with x(1) in both the numerator and denominator (you may need to put this back - please check your course work to be sure). Some similar changes would need to be made to your other three partial derivative functions. Try doing this and see what happens!