Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Implementation of Newton Raphson and Modified Newton Raphson Method by Mathematica

This article discusses the Mathematica approaches to code Newton Raphson (NR) and Modified Newton Raphson (MNR) and the differences if any in computation. Let's assume that stress and strain of a rigid bar is defined by the following relationship:

𝜎 = 100000 (√𝜀 + (1/80) sin (𝜋𝜀/0.002))

Now, taking the first derivative of the above relationship, we may write,

𝜎' = 100000 ((1/2√𝜀) + 19.653 cos (1570.80𝜀))

The NR method implements an iterative approach by considering a tangent to a curve, for which the function under consideration becomes zero. If f(xⱼ) is any function and f′(xⱼ) is the derivative of that function at point xⱼ, then according to NR approach, we may write as,

xⱼ₊₁ = xⱼ - f(xⱼ) / f′(xⱼ)

The MNR works in a similar way, but is helpful for converging a solution efficiently for roots with multiplicity and may be represented as,

xⱼ₊₁ = xⱼ - ⟮ f(xⱼ)  f′(xⱼ)⟯ /⟮ f′(xⱼ)² -  f(xⱼ) f′′(xⱼ)⟯



The following Mathematica program is developed to plot the function and its derivatives.

(*Function Defined in the Problem*)
sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]);
y1 = 4000;
y2 = 6000;
sigmap = D[sigma, epsilon];
sigmap2 = D[sigma, {epsilon, 2}];

(*Plotting the Function, its Derivatives, Upper and Lower Bounds*)
A1 = Plot[{sigma, y1, y2}, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma"}, PlotStyle -> {Red, Green, Black, Thick}]
A2 = Plot[sigmap, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma'"}, PlotStyle -> {Red, Thick}]
A3 = Plot[sigmap2, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma''"}, PlotStyle -> {Red, Thick}]

(*Initial Guesses*)
sigma /. epsilon -> 7.928*10^-4
sigma /. epsilon -> 7.924*10^-4
sigma /. epsilon -> 3.882*10^-3
sigma /. epsilon -> 3.879*10^-3


plotting the function

plotting the function

plotting the function

function results














Mathematica Program for the NR Method:


(*Mathematica Codes for Newton-Raphson Approach*)
Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;
SigmaDerivative = D[Sigma, epsilon];
epsilontable = {0.0002};
Error = {1};
SetTolerance = 0.0005;
MaximumIteration = 100;
i = 1;

(*Newton-Raphson Algorithm Implementation*)
While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],
  epsilonnew = 
   epsilontable[[i]] - (Sigma /. 
       epsilon -> epsilontable[[i]])/(SigmaDerivative /. 
       epsilon -> epsilontable[[i]]); 
  epsilontable = Append[epsilontable, epsilonnew]; 
  Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 
  Error = Append[Error, Errornew];
  i++];

L = Length[epsilontable];
SolutionTable = 
  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, L}];
SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};
L = Prepend[SolutionTable, SolutionTable1];
Print["Showing Results for the Newton-Raphson Approach when Sigma is \
4,000 and Initial Guess is 0.0002"]
ScientificForm[L // MatrixForm, 4]



Mathematica Program for the MNR Method:

(*Mathematica Codes for the Modified Newton-Raphson Approach*)
Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;
epsilontable = {0.0002};
SigmaDerivative = D[Sigma, epsilon] /. epsilon -> epsilontable[[1]];
Error = {1};
SetTolerance = 0.0005;
MaximumIteration = 100;
i = 1;

(*Implementation of the Modified Newton-Raphson Algorithm*)
While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],
  epsilonnew = 
   epsilontable[[i]] - (Sigma /. epsilon -> epsilontable[[i]])/
     SigmaDerivative; epsilontable = Append[epsilontable, epsilonnew];
   Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 
  Error = Append[Error, Errornew];
  i++];

T = Length[epsilontable];
SolutionTable = 
  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};
T = Prepend[SolutionTable, SolutionTable1];
Print["Showing Results for the Modified Newton-Raphson Approach when \
Sigma is 4,000 and Initial Guess is 0.0002"]
ScientificForm[T // MatrixForm, 4]


Results:

The following table shows the comparison of results from NR and MNR methods. We see that at iteration number 17, the Modified Newton-Raphson approach converges, whereas the Newton-Raphson approach is still iterating. However, it is noticeable that the step is faster than the Modified approach. Nevertheless, Newton-Raphson approach fails to converge early while the Modified approach, although iterates slowly, but converges eventually.

Comparison between Newton-Raphson and Modified Newton-Raphson Approach























#NewtonRaphson #Mathematica #MatrixForm #Blog #Blogger

1 comment: