1
$\begingroup$

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.

$\endgroup$
12
  • $\begingroup$ Sorry, Newton.m is including just a function Jac[f_, x_] := Outer[D, f, x]; $\endgroup$
    – drxy
    Feb 6, 2015 at 7:51
  • $\begingroup$ Ok, then please change the code, and just define the function in the code. The idea is to allow one to run the code. $\endgroup$
    – Nasser
    Feb 6, 2015 at 7:53
  • $\begingroup$ Ok @Nasser. I did it. Now you can run the code in your computer. $\endgroup$
    – drxy
    Feb 6, 2015 at 7:58
  • $\begingroup$ 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? $\endgroup$
    – Nasser
    Feb 6, 2015 at 8:22
  • $\begingroup$ 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 $\endgroup$
    – drxy
    Feb 6, 2015 at 8:29

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ It works now. But i will do some modifications and will publish the code here. Thank you ;) $\endgroup$
    – drxy
    Feb 6, 2015 at 12:44
  • $\begingroup$ @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)$. $\endgroup$
    – Amzoti
    Feb 7, 2015 at 21:13
  • $\begingroup$ @drxy I am leaving soon on a trip and shall not be able to think about this for several days. Sorry. $\endgroup$
    – bbgodfrey
    Feb 7, 2015 at 22:28
  • $\begingroup$ @bbgodfrey friend its amzoti's comment. I will try to solve his problem. Thank you. Good luck at your job. $\endgroup$
    – drxy
    Feb 7, 2015 at 22:31
4
$\begingroup$

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}}
*)
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – drxy
    Feb 7, 2015 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.