0
$\begingroup$

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
$\endgroup$
5
  • 2
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Dec 1, 2018 at 3:23
  • 3
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Dec 1, 2018 at 3:30
  • $\begingroup$ 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. $\endgroup$
    – H. Alanzi
    Dec 1, 2018 at 4:15
  • $\begingroup$ Perhaps, MATLink may help use Matlab code in ·Mathematica more than translating scritps. see matlink.org. $\endgroup$
    – Jerry
    Dec 1, 2018 at 15:50
  • $\begingroup$ @H.Alanzi What should be the result? $\endgroup$ Dec 1, 2018 at 16:43

1 Answer 1

1
$\begingroup$

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}*)
$\endgroup$

Your Answer

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

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