Translating the Newton-Raphson method from Matlab to Mathematica

Question

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

Answer 1

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[[1]] = ((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[[2]] = ((kx*1.127`*10^-3*dy*h)/(dx*mu*Bo))*(p1 - 
      p2) - ((h*dy*dx*phi*co)/(Bo*5.615*dt))*(p2 - pi2);
 R[[3]] = 
  q + ((2*Pi*kx*1.127`*10^-3*h)/(mu*Bo*Log[re/rw]))*(p1 - pwf1);

 Delta = LinearSolve[A, -R];

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


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