[Tutor] newton's method for system of nonlinear equations

[Dave Angel]

kgotlelelok at galmail.co.za wrote:

>
> Hi,
> I am trying to write a program in python that solves a system of nonlinear
> equations using newton's method. I don't know what I am doing wrong.
> Please help
>
>   
First thing wrong is posting two messages in the mailing list in rapid 
succession (7 minutes) with same content and different titles.  Not the 
best way to get sympathy.
Next:  paste your error log, or if there is none, show the results and 
explain how they differ from what you wanted.  Also show versions of 
both Python and the scipy library you're importing.
Third:  you have lots of obvious things wrong with the code itself, but 
it's usually best to start fixing it where the error message points;  
I'm not familiar with scipi, but I'll try to give you a few pointers anyway.
Fourth:  you use the "import *" form, which is generally frowned upon 
for all but trivial scripts.  It pollutes the global namespace, and 
possibly masks built-in functions and variables.
Fifth:  your names aren't very mnemonic, so somebody not familiar with 
the code has to spend longer just guessing what things are for.

Now for specifics in the code:
> from scipy import*
>
> x = array([0.0,0.0,0.0])
> n=len(x)
> tol= 0.00001
> N=30
>
> k=1
> while k <= N:
>   
Nowhere in this loop do you modify k or N, so it'll never terminate.
>     def f(x):
>   
Defining a function inside a loop is seldom useful, and not a technique 
for beginners.  In fact, a beginner should expect all functions to be at 
module level (left margin), and all methods to be at class level (one 
indent in).  Besides, you never call it or the nested function J.
>         f= zeros((len(x)),float)
>   
It's bad form to re-use the function name as a local within the 
function.  It works, but probably contributes to your confusion below.
>         f[0][k]= x[0][k]**2 + x[1][k]-37
>         f[1][k]=x[0][k]- x[1][k]**2- 5
>         f[2][k]= x[0][k] + x[1][k]+ x[2][k]- 3
>         return f[k]
>         def J(x):
>   
This definition is never reached, so the code is dead.
>             J= zeros((n,n),float)
>             for i in range(n):
>                 ei=zeros(n,float)
>                 ei[i]=1.0
>                 J[:i]=(f(x[k]+tol*ei)-f(x[k]))/tol
>                 return J
>
>         y[k] = -(J.I)*f[k]
>   
Here, you're misusing J and f.  If you're trying to call those 
functions, your syntax is way-off.  But if you're trying to use the 
local variables by those names, realize that the local J is not in scope.
>         x[k+1]=x[k]+y[k]
>   
You'll never get to these next lines, because the previous while loop 
never terminates.  Perhaps an indentation problem?
> if sqrt(dot(f0,f0)/len(x)) < tol: print x
> else:
>     k=k+1
>   
Not much point in incrementing a loop variable outside the loop.
> print 'Too many iterations'
>
>
>   

There are other things to ask, such as whether you intended to quit 
after printing x.  But the fix for that depends on other things you'll 
fix first.