# Translating the Newton-Raphson method from Matlab to Mathematica

I am trying to optimize the variables of two (or three depending on how you think about it) matrices using the Newton-Raphson Method. I made the bottom code in Matlab. The code is pretty simple it uses a while loop with the Newton-Raphson over a number of equations until I get a fixed point or value. I was wondering how I can write a similar code in Mathematica.

This is a slightly off-topic, However, generally, what are the best practices when dealing with matrices in Mathematica? As I am finding that Matrices are the hardest thing to wrap my head around when it comes to coding in Mathematica.

Note: in the code, inside the while loop, the equations just needs to be rearranged if one wants to copy and paste this sample code.

clear all

q = -500;
dx = 500;
dy = 500;
mu = 1;
Bo = 1;
phi = 0.25;
co = 10e-5;
dt = 30;
re = 0.14*sqrt(dx^2+dy^2);
rw = 0.3;

pi1 = 4000; %pi_n
pi2 = 4000; %p2_n
p2 = 2722; %p1_n+1
pwf1 = 2308;

%guessed values
p1 = 2722; %p1_n+1
kx = 70;
h = 250;

Delta = ones(3,1);
A = zeros(3,3);
R = zeros(3,1);

while max(abs(Delta)) > 10e-8
A(1,1) = (kx*1.127e-3 * dy)/(dx) * (1/(mu*Bo))*(p2-p1) -
((dx*dy*phi*co)/(Bo*5.615*dt))*(p1-pi1);
A(1,2) = ((dy*h)/dx) * (1/(mu*Bo)) * (p2-p1);
A(1,3) = -((kx*1.127e-3*dy*h)/(dx))*(1/(mu*Bo)) -
((h*dy*dx*phi*co)/(Bo*5.615*dt));

A(2,1) = ((kx*1.127e-3*dy)/(dx))* (1/(mu*Bo))*(p1-p2) -
((dx*dy*phi*co)/(Bo*5.615*dt))* (p2-pi2);
A(2,2) = ((dy*h)/(dx)* (1/mu*Bo)) * (p1-p2);
A(2,3) = (kx*1.127e-3*dy*h)/(dx)* (1/ (mu*Bo));

A(3,1) = ((2*pi*kx*1.127e-3)/(mu*Bo*log(re/rw)))*(p1-pwf1);
A(3,2) = ((2*pi*h)/(mu*Bo*log(re/rw)))*(p1-pwf1);
A(3,3) = (2*pi*h*kx*1.127e-3)/(mu*Bo*log(re/rw));

R(1,1) = ((kx*1.127e-3*dy*h)/(dx*mu*Bo)) * (p2-p1) + q -
((h*dy*dx*phi*co)/(Bo*5.615*dt))* (p1-pi1);
R(2,1) = ((kx*1.127e-3*dy*h)/(dx*mu*Bo)) * (p1-p2) -
((h*dy*dx*phi*co)/(Bo*5.615*dt))* (p2-pi2);
R(3,1) = q + ((2*pi*kx*1.127e-3*h)/(mu*Bo*log(re/rw)))* (p1 -
pwf1);

Delta = (A\(-R));

h = h + Delta(1,1);
kx = kx + Delta(2,1);
p1 = p1 + Delta(3,1);

end

• I can write a similar code in Mathematica either using Module, FindRoot, or a special Mathematica function if you do not want to translate the code literally as is, and can use special functions, then why not just use NSolve? It is one call. So it is not clear if you are looking to translate Matlab code, or asking how to solve nonlinear equations in Mathematica. – Nasser Dec 1 '18 at 3:23
• BTW, there are number of posts on this site that shows how to write NR in Mathematica. For example code-for-newtons-method and implementing-the-newton-raphson-method-for-finding-the-zeros-of-a-function and newton-raphson-method-in-mathematica may be the code in these is something you can use. – Nasser Dec 1 '18 at 3:30
• Thank you for sharing those links with me! To answer your concerns, I mainly just want to know how to how to translate this particular code as it is slightly more complicated to deal with than a single variable polynomial. – H. Alanzi Dec 1 '18 at 4:15
• Perhaps, MATLink may help use Matlab code in ·Mathematica more than translating scritps. see matlink.org. – Jerry Dec 1 '18 at 15:50
• @H.Alanzi What should be the result? – Alex Trounev Dec 1 '18 at 16:43

Not one to one but close. I used fast code with Do

 ClearAll["*"];

q = -500;
dx = 500;
dy = 500;
mu = 1;
Bo = 1;
phi = 0.25;
co = 10^-4;
dt = 30;
re = 0.14*Sqrt[dx^2 + dy^2];
rw = 0.3;

pi1 = 4000;(*Pi_n*)
pi2 = 4000;(*p2_n*)
p2 = 2722;(*p1_n+1*)
pwf1 = 2308;

(*guessed values*)
p1 = 2722;(*p1_n+1*)
kx = 70;
h = 250;

Delta = {1, 1, 1};
A = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
R = {0, 0, 0};

n = 10^6;
Do[A[[1, 1]] = (kx*1.127*10^-3*dy)/(dx)*(1/(mu*Bo))*(p2 -
p1) - ((dx*dy*phi*co)/(Bo*5.615*dt))*(p1 - pi1);
A[[1, 2]] = ((dy*h)/dx)*(1/(mu*Bo))*(p2 - p1);
A[[1, 3]] = -((kx*1.127*10^-3*dy*h)/(dx))*(1/(mu*Bo)) - ((h*dy*dx*
phi*co)/(Bo*5.615*dt));

A[[2, 1]] = ((kx*1.127*10^-3*dy)/(dx))*(1/(mu*Bo))*(p1 -
p2) - ((dx*dy*phi*co)/(Bo*5.615*dt))*(p2 - pi2);
A[[2, 2]] = ((dy*h)/(dx)*(1/mu*Bo))*(p1 - p2);
A[[2, 3]] = (kx*1.127*10^-3*dy*h)/(dx)*(1/(mu*Bo));

A[[3, 1]] = ((2*Pi*kx*1.127*10^-3)/(mu*Bo*Log[(re/rw)]))*(p1 - pwf1);
A[[3, 2]] = ((2*Pi*h)/(mu*Bo*Log[(re/rw)]))*(p1 - pwf1);
A[[3, 3]] = (2*Pi*h*kx*1.127*10^-3)/(mu*Bo*Log[(re/rw)]);

R[] = ((kx*1.127*10^-3*dy*h)/(dx*mu*Bo))*(p2 - p1) +
q - ((h*dy*dx*phi*co)/(Bo*5.615*dt))*(p1 - pi1);
R[] = ((kx*1.127*10^-3*dy*h)/(dx*mu*Bo))*(p1 -
p2) - ((h*dy*dx*phi*co)/(Bo*5.615*dt))*(p2 - pi2);
R[] =
q + ((2*Pi*kx*1.127*10^-3*h)/(mu*Bo*Log[re/rw]))*(p1 - pwf1);

Delta = LinearSolve[A, -R];

h = h + Delta[];
kx = kx + Delta[];
p1 = p1 + Delta[];
If[Norm[Delta] >= 10^-7, Continue[], Break[]];, {i, 1, n}]

Delta
(*{1.15441*10^-15, -9.99819*10^-8, 1.99263*10^-13}*)