Newton-Raphson Method for Non-linear System of 3 variables in Matlab

I am trying to solve 3 non-linear system of 3 variables using the newton-raphson method in matlab. Here are the 3 non-linear equations:

 c[alpha I+ k_f+k_d+k_ns+k_p(1-q)]-I alpha =0 
 s[lambda_b c P_C +lambda_r (1-q)]- lambda_b c P_C =0 
 q[gamma +c k_p (P_C/P_Q)]- c k_p (P_C/P_Q)=0 

I need to find the values of c,s, and q  using the newton-raphson method.

=> can someone please check my code, there are no errors so, its converges after six iterations. can give the accurate values of c,s and q.

I am only concern about this line in my matlab code:

disp(sprintf('iter=%6.15f, c=%6.15f, s=%6.15f, q=%6.15f', iter,xnew));

I want to display xnew as new and accurate values for c,s and q. Is this code giving the correct and update values of c,s,q as xnew? Sorry if my question looks silly. Thanks in advance.

This is my matlab code :

    format long
    clear;

    %values of parameters
    I=1200;
    k_f= 6.7*10.^7;
    k_d= 6.03*10.^8; 
    k_n=2.92*10.^9; 
    k_p=4.94*10.^9;
    lambda_b= 0.0087;
    lambda_r =835; 
    gamma =2.74; 
    alpha =1.14437*10.^-3;
    P_C= 3 * 10.^(11);
    P_Q= 2.87 * 10.^(10);

    tol = 10.^-4;  %tol is a converge tolerance

    %initial guess or values
  c=1; 
  s=0.015;
  q=0.98; 
  iter= 0; %iterations
  xnew =[c;s;q];
  xold = zeros(size(xnew));
  while norm(xnew - xold) > tol
      iter= iter + 1;
      xold = xnew;

      % update c, s, and q
      c = xold(1);
      s = xold(2);
      q = xold(3);

      %Defining the functions for c,s and q.
      f = c * (alpha*I + k_f + k_d + k_n * s + k_p*(1-q))-I *alpha;
      g = s * (lambda_b * c* P_C + lambda_r *(1-q))- lambda_b* c * P_C; 
      h = q * ( gamma + c * k_p *(P_C / P_Q))- (c * k_p * (P_C / P_Q));

      %Partial derivatives in terms of c,s and q.
      dfdc = alpha*I + k_f + k_d + k_n * s + k_p*(1-q);
      dfds = k_n *c ;
      dfdq = - k_p *c;

      dgdc = lambda_b * P_C *(s-1);
      dgds = lambda_b * c* P_C + lambda_r *(1-q);
      dgdq = - lambda_r * s;

      dhdc = k_p *(P_C / P_Q)*(q-1);
      dhds = 0;
      dhdq = gamma + c * k_p *(P_C / P_Q);

      %Jacobian matrix 
      J = [dfdc dfds dfdq; dgdc dgds dgdq; dhdc dhds dhdq];    
      % Applying the Newton-Raphson method
      xnew = xold - J\[f;g;h];
      disp(sprintf('iter=%6.15f,  c=%6.15f,  s=%6.15f, q=%6.15f', iter,xnew)); 
  end

Answer 1

[Roger Stafford]

I do see an error in your code, though it does not account for why you get only one iteration. You are not renewing the value of x0 inside the while-loop and therefore it continues to be the original estimate values. Also you are not renewing the values of c, s, and q from your 'xnew' values. The way the code reads, if it repeats for a second iteration, the while-loop should run forever, since nothing is ever changed after that.

[dave]

Thank you very much for you comment. I see it now and just change it. now, it converges after six iterations. I am only concern about this line in my matlab code:

disp(sprintf('iter=%6.15f, c=%6.15f, s=%6.15f, q=%6.15f', iter,xnew));

I want to display xnew as new and accurate values for c,s and q. Is this code giving the correct and update values of c,s,q as xnew?

[Roger Stafford]

Your code correction looks good to me.

The 'disp' line should give you whatever values you have at this points in 'xnew'. Even though the c, s, and q variables are not yet updated at that point 'sprintf' is not influenced by the "c=", "s=", and "q=" strings it prints as far as the values that follow them. Those value come strictly from 'xnew'. You could have printed a crazy "inf=" and it would still give the values that are currently in 'xnew'.

By the way, when you write"%6.15f" the spacing in the result may not always be what you expect. The number 6 is understood to be the total length used to express the number's entire value (unless more are required,) not just the part to the left of the decimal point. If you want to allow for, say, a possible minus sign, five digits to the left of the decimal point, the decimal point itself, and 15 places to its right, you should use "%22.15f", because 22 = 1+5+1+15. However, possibly you want it to take only as many places as it needs in which case "%6.15f" is okay. You could also use "%.15f"

[dave]

Thank you very much for your help sir. Really appreciated :)

[Ankit Barthwal]

can you please send the code after corrections, I have to solve 3 non-linear equations as well.