Newton Rhapson - this one

Posted on 2008-04-14 13:31 By bongiman (

24)

Categories:

Math Advice (Student)

I decided to repost this as the previous one went crazy.......

I'm trying to construct an iterative procedure NR1(f,x0,N) where f is a function, x0 is an initial estimate (which can be complex) and N the number of iterates.

Given: x_k+1=x_k - f(x_k)/f'(x_k), k=0,1,2....,

Now for example by defining f:= x-> x^2-2 and inputting x₀=1 and N=10 i should be able to check that my procedure calculates sqrt(2) correctly.

This is what i've got so far.....

NR1:= proc(f,x0,N)
> local x,k,f:
> x := x0:
> for k to N do:
> x := x-f/D(f)
> end do:
> evalf(x):
> end proc:f:= x-> x^2-2:
> NR1(f,1,10);

Could someone point out where i'm going wrong?

Thanks.

Thank you for your response.

Comment on 2008-04-19 17:06 from bongiman (

24)

Thank you for your response. I have a second query:

If i want to modify NR1 above to NR2(f,x0,N,eps) where eps is a real number so that it returns the first iterate x_kfor which |x[k]-x[k-1]|<eps and k<=N. However, if |x[k]-x[k-1]|>=eps for k=1,...,N, then an error message should be displayed.

If the procedure is correct then for f:= x-> x^5-1 i should get an answer when x₀=0.60+0.20i and the error message when x₀=0.60+0.16i.

The code i've written so far is incorrect but i'm not sure where i'm going wrong. Would you be able to give me some pointers?

NR2:=proc(f::mathfunc,x0::complex,N::posint,eps)
> local x,k:
> x := x0:
> for k to N do:
>   if abs(x[k]-x[k-1]) >= eps then
>   print("Convergence has not been achieved!"):
>   return :
> end if:
> x := x-f(x)/D(f)(x)
> end do:
> evalf(x):
> end proc:f:= x-> x^5-1:
> NR2(f,0.6+I*0.16,10,0.00001);
Error, (in NR2) cannot determine if this expression is true or false: .1e-4 <= abs((.6+.16*I)[1]-(.6+.16*I)[0])

> NR2(f,0.6+I*0.20,10,0.00001);
Error, (in NR2) cannot determine if this expression is true or false: .1e-4 < abs((.6+.20*I)[1]-(.6+.20*I)[0])

Thanks in advance.

suggestions

Comment on 2008-04-19 18:18 from acer ( Maple Rating 4

1377)

You need to make sure that x[k] and x[k-1] have been assigned a numeric value, for each time that they are compared and for the value of k at that moment. Make sure that you stick with either a scheme with indexed x[k],x[k-1],etc or a scheme with x,xnew,xold,etc. You were mixing x and indexed x[k], which wouldn't work.

Also, you indicated that you only wanted to print the message if convergence failed for all k from 1 to N, so put it after the loop and not inside the loop.

NR2:=proc(f::mathfunc,x0::complex,N::posint,eps)
> local x,k:
>   x[0] := x0:
>   for k to N do:
>     x[k] := evalf( x[k-1]-f(x[k-1])/D(f)(x[k-1]) );
>   end do;
>   if abs(x[N]-x[N-1]) >= eps then
>     printf("Convergence has not been achieved after %a iterations!\n",N);
>   else
>     return x[N];
>   end if;
> end proc:
>
> f:= x-> x^5-1:
>
> NR2(f,0.6+I*0.6,10,0.00001);
Convergence has not been achieved after 10 iterations!
> NR2(f,0.2+I*0.6,10,0.00001);
                         0.3090169944 + 0.9510565163 I

Side tip: maple's for-loop counters finish with a value incremented one step more than the last used value, when they have finished. For example, a for-loop counting k from 1 to 10 will have value 11 after it's finished. This matters, if you plan to refer to x[k] after it's finished. Notice that I referred to x[N] after the loop. I could also have referred to x[k-1]=x[10] but not to x[k]=x[11] which is unassigned.

Lastly, Robert's suggestion to use evalf was so that a large (potentially huge) symbolic expression did not accumulate via the iterative process. Using evalf can cure that, but only if it's done prior to assigning to x or x[k]. You had it done only as a separate task afterwards. I put it right in the iterative step.

acer

Brilliant

Comment on 2008-04-20 08:54 from bongiman (

24)

Thank you for your answer. Very well explained and extremely helpful. I've only been working on Maple for about 3 weeks now and it makes me wonder how on earth i'm actually going get my head around it all :).

minor point

Comment on 2008-04-20 17:33 from acer ( Maple Rating 4

1377)

The printed message about nonconvergence seems to make most sense outside and after the loop. That's why I moved it outside the loop, because of that purpose. But a test on abs(x[k]-x[k-1]) is not out of place inside the loop, for another purpose. It would also be quite sensible to put such a check inside the loop so that an early return could be made if convergence occurred while k was still less than N. You may not want it to continue with the full N iterations if the tolerance has already been met.

acer

Newton-Raphson

Comment on 2008-04-14 14:47 from Robert Israel ( Maple Rating 4

1460)

There are several problems with your procedure.

1) You have both a formal parameter f and a local variable f.
Get rid of the local variable.

2) You need f(x)/D(f)(x) instead of f/D(f).

3) It's probably a good idea to use floating point throughout on this.
I've seen students crash Maple in a Newton-Raphson loop, getting rational numbers
with huge numerators and denominators. So make it

x := evalf(x - f(x)/D(f)(x));

Newton's method again

Comment on 2008-05-24 15:49 from adnan06 (

13)

Hi all,

The discussion here is very interesting and useful. I still have a question regarding the Newton's method iteration or any other iteration. how can I implement the proc above for a function like : f(x) = sqrt(x-r) for x>= r and f(x) = - sqrt(r-x) for x<= r, where r is any ( choosen) real number. I know this function does not converge. and I want to see the result as output=plot.

Thanks to all

output = plot

Comment on 2008-05-25 02:17 from Robert Israel ( Maple Rating 4

1460)

The simplest way to define your function is

> f := x -> piecewise(x>=r, sqrt(x-r), -sqrt(r-x));

The Newton procedure for one iteration:

> Newt:= x -> evalf(x - f(x)/D(f)(x));

Now the way I like to visualize this is by drawing a vertical line from [x[i], 0] to [x[i], f(x[i])] and then the tangent line from there to [x[i+1],0]. The following procedure will do this.

> NewtStep:= x -> plot([[x,0],[x,f(x)],[Newt(x), 0]], colour=blue);

In your particular case, it turns out that the Newton function is an involution: Newt(Newt(x)) = x. Here's a plot:

> r:= 1; 
  plots[display]([plot(f(x),x=-1..3), NewtStep(2.5),
    NewtStep(Newt(2.5))]);