How to write Newton-Raphson module and return a value?

Question

I am writing a code which solve nonlinear algebraic systems via Newton-Raphson algorithm. I want to make a Newton Module with another file and use the Module in my test file, but i didn't achieve it. Because, If the method is not converged, I break the code and return "not converged". Hence, iteration is not continuing and my code is not working properly. Here is my code:

F[y1_,y2_] = {
((y1 - 11)^2 / 64) - ((y2-7)^2/100) -1,
(y1-3)^2 + (y2-1)^2 -400
}; (*Define nonlinear algebraic system vector *)

J[y1_,y2_] = Outer[D,F[y1,y2],{y1,y2}]; (* Find symbolic jacobian *)

{y10,y20}= {20,-4}; (* Initial values for newton iterations *)

Tol   = 10^-12;
counter = 0; (*How many newton iterations*)
DeltaY = {0.1,0.1};
MaxIter = 50;
list  = Table[{0,0},{i,1,MaxIter}];
(* Begin loop for newton iterations *)
i=1;
While[Norm[ DeltaY, 2] > Tol,
    Result = Solve[J[y10,y20].{{dy1},{dy2}}==-F[y10,y20],{dy1,dy2}] //N //Flatten ;  (*Solve linear system for Delta y *)
    DeltaY = {dy1,dy2}/.Result;
    {y10,y20} = {y10,y20} + DeltaY; (*Find new Y values with computed delta y solutions *)
    list[[i]] = {y10,y20}; (* Append new computed Y elements into list array *)
    Print[list[[i]]," and error:",DeltaY]       
    Print[Norm[ DeltaY, Infinity]]
    If[i>MaxIter,Print["Not converged"] Return[{0,0}];]
counter++
i++;
] //AbsoluteTiming

How can I make n dimensional newton iteration module without breaking code. Also how to return a value with a module.

Note: I will use my newton module for some implicit numerical schemes.

Best regards.

Sorry, Newton.m is including just a function
Jac[f_, x_] := Outer[D, f, x]; — drxy, Feb 06 '15 at 07:51
Ok, then please change the code, and just define the function in the code. The idea is to allow one to run the code. — Nasser, Feb 06 '15 at 07:53
Ok @Nasser. I did it. Now you can run the code in your computer. — drxy, Feb 06 '15 at 07:58
Mean while, the API to the function should be very clear: What is the input to the function should be, and what is the output. how will you be calling it? What is should return? — Nasser, Feb 06 '15 at 08:22
I will give the parameters which are
Function,Jacobian,Max Iteration, Error Tolerance, Dimension. Our module will seems like NewtonMethod[Function,Jacobian,Max Iteration, Error Tolerance, Dimension] = Solution vector or false — drxy, Feb 06 '15 at 08:29
Your code certainly runs and produces results, although I have not verified their accuracy. For clarity, are you simply trying to convert the code into a Module, or something else? Also, are you certain that FindRoot will not meet your needs? — bbgodfrey, Feb 06 '15 at 11:30
My friend @bbgodfrey, I couldn't convert my above code into a Module. I have already ask how can i do it?. FindRoot will meet my needs but i want to control everything while solving a problem.
Briefly, I want to write an implicit Runge-Kutta package with Mathematica like this matlab script https://bitbucket.org/drreynolds/rklab. If I achieve this, I can publish a paper or book to help people about numerical computations on implicit methods via Mathematica.

Thanks in advance. — drxy, Feb 06 '15 at 11:52
Since you plan to use this and it's not just an exercise, why not use the built-in Newton's method function, FindRoot? — Michael E2, Feb 06 '15 at 13:47
Dear @MichaelE2, as I say above, I want to control everything in my iteration scheme. FindRoot is written by another person. If I want to add something to FindRoot method, I think I can not do it. So, I want to write my own code from zero to hero ;) — drxy, Feb 06 '15 at 14:04
I see. (BTW, FindRoot does everything you mention you want to do in the question. You can use Check to see if it complains about convergence. There are also the options EvaluationMonitor and StepMonitor, from which you can Throw exceptions. FindRoot does do some other things automatically or mysteriously, like check step size, so I can understand if you want to write your own. ) — Michael E2, Feb 06 '15 at 14:14

score 4 · Answer 1 · answered Feb 06 '15 at 21:41

Why not make use of Mathematica built-in functional capabilities? The following code can be modified to add whatever stopping condition is desired.

   jac[f_, vars_] := Outer[D, f[Sequence @@ vars], vars]
   jacobian[f_, vars_, pt_] := jac[f, vars] /. Thread[vars -> pt]
   newt[f_, vars_, pt_] := pt - Inverse[jacobian[f, vars, pt]].f[Sequence @@ pt]

   f[x_, y_] := {x y - 3, x^2 + y^2 - 9}
   NestList[newt[f, {x, y}, #] &, {3., 1.}, 5]

(*
    {{3.,1.},{2.8125,1.0625},{2.80256,1.07042},{2.80252,1.07047},{2.80252,1.07047},{2.80252,1.07047}}
*)

Dear @murray, your way is highly efficient. But I am lack of functional programming now. I will work Mathematica's functional ways if I find an appropriate time. Thank you very much. — drxy, Feb 07 '15 at 11:35

score 2 · Accepted Answer · answered Feb 06 '15 at 12:32

2

The simplest code to create your module is

nr[func_, jac_, y1_, y2_, MaxIter_, Tol_] := 
 Module[{y10 = y1, y20 = y2, counter = 0, DeltaY = {0.1, 0.1}, 
   list = Table[{0, 0}, {i, 1, MaxIter}], i = 1},
  While[Norm[DeltaY, 2] > Tol, 
   Result = Solve[jac[y10, y20].{{dy1}, {dy2}} == -func[y10, y20], {dy1, 
        dy2}] // N // Flatten;
   DeltaY = {dy1, dy2} /. Result;
   {y10, y20} = {y10, y20} + DeltaY;
   list[[i]] = {y10, y20};
   Print[list[[i]], " and error:", DeltaY, "  ", i, "  ", counter] ;
   Print[Norm[DeltaY, Infinity]]; 
   If[i > MaxIter, Print["Not converged"] Return[{0, 0}];] ;
   counter++ ;
   i++;]; {y10, y20}]

Call it with

nr[F, J, y10, y20, MaxIter, Tol]

and it will produce the Print statements produced by the code in your question, plus the answer, {22.519, -3.35977}. Of course, your final product should contain options, documentation, additional error checks, and the like. Hope this helps.

answered Feb 06 '15 at 12:32

bbgodfrey

61,439
17
89
156

It works now. But i will do some modifications and will publish the code here. Thank you ;) – drxy Feb 06 '15 at 12:44
@bbgodfrey: It would be great if it could do it for functions of general $n$ dimensional size. However, even being able to add 3D examples like $f(x,y,z)$ would be great as it currently only supports @D examples $f(x,y)$. – Amzoti Feb 07 '15 at 21:13
@drxy I am leaving soon on a trip and shall not be able to think about this for several days. Sorry. – bbgodfrey Feb 07 '15 at 22:28
@bbgodfrey friend its amzoti's comment. I will try to solve his problem. Thank you. Good luck at your job. – drxy Feb 07 '15 at 22:31

How to write Newton-Raphson module and return a value?

2 Answers2

Linked