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

MATLABnonlinear optimizationnonlinear system

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:

\begin{equation} c[\alpha I+ k_f+k_d+k_ns+k_p(1-q)]-I \alpha =0 \end{equation}
\begin{equation} s[\lambda_b c P_C +\lambda_r (1-q)]- \lambda_b c P_C =0 \end{equation}
\begin{equation} q[\gamma +c k_p \frac{P_C}{P_Q}]- c k_p \frac{P_C}{P_Q}=0 \end{equation}

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

=>
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;
x0= [c;s;q];

iter= 0; %iterations
xnew =[100;100;100];
while norm(xnew -x0) > tol
    iter= iter + 1;
%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 = x0 - J\[f;g;h];
disp(sprintf('iter=%6.15f,  c=%6.15f,  s=%6.15f, q=%6.15f', iter,xnew)); 
end

can someone please check my code, there are some errors so, its not working. Thanks in advance.

Best Answer

Here is Newton's method in general in MATLAB code. Define a function J that takes a vector x as its argument and returns the Jacobian; likewise, define a function f that takes a vector x as its argument and represents the vector-valued function. For instance, if $$\mathbf{f}(x,y) = \begin{pmatrix} x^2 + y^2 \\ 2x-y \end{pmatrix}$$ then we would write:

f = @(x)[ x(1)^2+x(2)^2;...
          2*x(1)-x(2)];
J = @(x)[ 2*x(1), 2*x(2);...
          2, -1];

Next, use this code to solve the system. I wrote this blindly; there might be a trivial error in there somewhere, as I am not near a MATLAB machine at the moment.

tol = 1e-4; % Or some other tolerance
err = 1000; % Any value larger than tol
x = x_initial; % However this is defined.
iter = 1; max_iter = 30; % Or whatever.
while (err > tol)
    delta_x = J(x)\(-f(x)); % Compute x_{n+1}-x_n
    err = norm(delta_x);
    x = x + delta_x;
    iter = iter + 1;
    [iter x'] % This line simply outputs the current iteration and the solution. You can dress this up by using sprintf if you like.
    if (iter > max_iter)
        disp 'Failed to converge';
        break;
    end
end
Related Question